File indexing completed on 2021-02-14 13:34:17
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011 #include "CommonTools/Statistics/interface/ChiSquaredProbability.h"
0012 #include "DataFormats/Common/interface/Handle.h"
0013 #include "DataFormats/EgammaTrackReco/interface/ConversionTrack.h"
0014 #include "DataFormats/EgammaTrackReco/interface/ConversionTrackFwd.h"
0015 #include "DataFormats/TrackCandidate/interface/TrackCandidateCollection.h"
0016 #include "DataFormats/TrackReco/interface/Track.h"
0017 #include "DataFormats/TrackReco/interface/TrackBase.h"
0018 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0019 #include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
0020 #include "FWCore/Framework/interface/ESHandle.h"
0021 #include "FWCore/Framework/interface/Event.h"
0022 #include "FWCore/Framework/interface/EventSetup.h"
0023 #include "FWCore/Framework/interface/stream/EDProducer.h"
0024 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0025 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0026 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0027 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0028 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0029
0030 #include <vector>
0031
0032 class ConversionTrackMerger : public edm::stream::EDProducer<> {
0033 public:
0034 explicit ConversionTrackMerger(const edm::ParameterSet& conf);
0035
0036 ~ConversionTrackMerger() override;
0037
0038 void produce(edm::Event& e, const edm::EventSetup& c) override;
0039
0040 private:
0041 edm::ParameterSet conf_;
0042
0043 edm::EDGetTokenT<reco::ConversionTrackCollection> trackProducer1;
0044 edm::EDGetTokenT<reco::ConversionTrackCollection> trackProducer2;
0045
0046 std::unique_ptr<reco::ConversionTrackCollection> outputTrks;
0047 };
0048
0049 #include "FWCore/Framework/interface/MakerMacros.h"
0050 DEFINE_FWK_MODULE(ConversionTrackMerger);
0051
0052 ConversionTrackMerger::ConversionTrackMerger(edm::ParameterSet const& conf) : conf_(conf) {
0053
0054 trackProducer1 = consumes<reco::ConversionTrackCollection>(conf_.getParameter<edm::InputTag>("TrackProducer1"));
0055 trackProducer2 = consumes<reco::ConversionTrackCollection>(conf_.getParameter<edm::InputTag>("TrackProducer2"));
0056
0057 produces<reco::ConversionTrackCollection>();
0058 }
0059
0060
0061 ConversionTrackMerger::~ConversionTrackMerger() {}
0062
0063
0064 void ConversionTrackMerger::produce(edm::Event& e, const edm::EventSetup& es) {
0065 double shareFrac = conf_.getParameter<double>("ShareFrac");
0066 bool allowFirstHitShare = conf_.getParameter<bool>("allowFirstHitShare");
0067 bool checkCharge = conf_.getParameter<bool>("checkCharge");
0068 double minProb = conf_.getParameter<double>("minProb");
0069
0070 int outputPreferCollection = conf_.getParameter<int>("outputPreferCollection");
0071 int trackerOnlyPreferCollection = conf_.getParameter<int>("trackerOnlyPreferCollection");
0072 int arbitratedEcalSeededPreferCollection = conf_.getParameter<int>("arbitratedEcalSeededPreferCollection");
0073 int arbitratedMergedPreferCollection = conf_.getParameter<int>("arbitratedMergedPreferCollection");
0074 int arbitratedMergedEcalGeneralPreferCollection =
0075 conf_.getParameter<int>("arbitratedMergedEcalGeneralPreferCollection");
0076
0077
0078
0079
0080
0081
0082 const reco::ConversionTrackCollection* TC1 = nullptr;
0083 static const reco::ConversionTrackCollection s_empty1, s_empty2;
0084 edm::Handle<reco::ConversionTrackCollection> trackCollection1;
0085 e.getByToken(trackProducer1, trackCollection1);
0086 if (trackCollection1.isValid()) {
0087 TC1 = trackCollection1.product();
0088
0089 } else {
0090 TC1 = &s_empty1;
0091 edm::LogWarning("ConversionTrackMerger") << "1st TrackCollection not found;"
0092 << " will only clean 2nd TrackCollection ";
0093 }
0094 reco::ConversionTrackCollection tC1 = *TC1;
0095
0096 const reco::ConversionTrackCollection* TC2 = nullptr;
0097 edm::Handle<reco::ConversionTrackCollection> trackCollection2;
0098 e.getByToken(trackProducer2, trackCollection2);
0099 if (trackCollection2.isValid()) {
0100 TC2 = trackCollection2.product();
0101
0102 } else {
0103 TC2 = &s_empty2;
0104 edm::LogWarning("ConversionTrackMerger") << "2nd TrackCollection not found;"
0105 << " will only clean 1st TrackCollection ";
0106 }
0107 reco::ConversionTrackCollection tC2 = *TC2;
0108
0109
0110 outputTrks = std::make_unique<reco::ConversionTrackCollection>();
0111 int i;
0112
0113 std::vector<int> selected1;
0114 for (unsigned int i = 0; i < tC1.size(); ++i) {
0115 selected1.push_back(1);
0116 }
0117 std::vector<int> selected2;
0118 for (unsigned int i = 0; i < tC2.size(); ++i) {
0119 selected2.push_back(1);
0120 }
0121
0122 std::map<reco::ConversionTrackCollection::const_iterator, std::vector<const TrackingRecHit*> > rh1;
0123 std::map<reco::ConversionTrackCollection::const_iterator, std::vector<const TrackingRecHit*> > rh2;
0124 for (reco::ConversionTrackCollection::const_iterator track = tC1.begin(); track != tC1.end(); ++track) {
0125 trackingRecHit_iterator itB = track->track()->recHitsBegin();
0126 trackingRecHit_iterator itE = track->track()->recHitsEnd();
0127 for (trackingRecHit_iterator it = itB; it != itE; ++it) {
0128 const TrackingRecHit* hit = &(**it);
0129 rh1[track].push_back(hit);
0130 }
0131 }
0132 for (reco::ConversionTrackCollection::const_iterator track = tC2.begin(); track != tC2.end(); ++track) {
0133 trackingRecHit_iterator jtB = track->track()->recHitsBegin();
0134 trackingRecHit_iterator jtE = track->track()->recHitsEnd();
0135 for (trackingRecHit_iterator jt = jtB; jt != jtE; ++jt) {
0136 const TrackingRecHit* hit = &(**jt);
0137 rh2[track].push_back(hit);
0138 }
0139 }
0140
0141 if ((!tC1.empty()) && (!tC2.empty())) {
0142 i = -1;
0143 for (reco::ConversionTrackCollection::iterator track = tC1.begin(); track != tC1.end(); ++track) {
0144 i++;
0145
0146
0147 selected1[i] = selected1[i] && outputPreferCollection != 0;
0148 track->setIsTrackerOnly(track->isTrackerOnly() && trackerOnlyPreferCollection != 0);
0149 track->setIsArbitratedEcalSeeded(track->isArbitratedEcalSeeded() && arbitratedEcalSeededPreferCollection != 0);
0150 track->setIsArbitratedMerged(track->isArbitratedMerged() && arbitratedMergedPreferCollection != 0);
0151 track->setIsArbitratedMergedEcalGeneral(track->isArbitratedMergedEcalGeneral() &&
0152 arbitratedMergedEcalGeneralPreferCollection != 0);
0153
0154 std::vector<const TrackingRecHit*>& iHits = rh1[track];
0155 unsigned nh1 = iHits.size();
0156 int j = -1;
0157 for (reco::ConversionTrackCollection::iterator track2 = tC2.begin(); track2 != tC2.end(); ++track2) {
0158 j++;
0159
0160
0161 selected2[j] = selected2[j] && outputPreferCollection != 0;
0162 track2->setIsTrackerOnly(track2->isTrackerOnly() && trackerOnlyPreferCollection != 0);
0163 track2->setIsArbitratedEcalSeeded(track2->isArbitratedEcalSeeded() &&
0164 arbitratedEcalSeededPreferCollection != 0);
0165 track2->setIsArbitratedMerged(track2->isArbitratedMerged() && arbitratedMergedPreferCollection != 0);
0166 track2->setIsArbitratedMergedEcalGeneral(track2->isArbitratedMergedEcalGeneral() &&
0167 arbitratedMergedEcalGeneralPreferCollection != 0);
0168
0169 std::vector<const TrackingRecHit*>& jHits = rh2[track2];
0170 unsigned nh2 = jHits.size();
0171 int noverlap = 0;
0172 int firstoverlap = 0;
0173 for (unsigned ih = 0; ih < nh1; ++ih) {
0174 const TrackingRecHit* it = iHits[ih];
0175 if (it->isValid()) {
0176 int jj = -1;
0177 for (unsigned jh = 0; jh < nh2; ++jh) {
0178 const TrackingRecHit* jt = jHits[jh];
0179 jj++;
0180 if (jt->isValid()) {
0181 if (it->sharesInput(jt, TrackingRecHit::some)) {
0182 noverlap++;
0183 if (allowFirstHitShare && (ih == 0) && (jh == 0))
0184 firstoverlap = 1;
0185 }
0186 }
0187 }
0188 }
0189 }
0190 int nhit1 = track->track()->numberOfValidHits();
0191 int nhit2 = track2->track()->numberOfValidHits();
0192
0193
0194 if ((noverlap - firstoverlap) > (std::min(nhit1, nhit2) - firstoverlap) * shareFrac &&
0195 (!checkCharge || track->track()->charge() * track2->track()->charge() > 0)) {
0196
0197
0198
0199 double probFirst = ChiSquaredProbability(track->track()->chi2(), track->track()->ndof());
0200 double probSecond = ChiSquaredProbability(track2->track()->chi2(), track2->track()->ndof());
0201
0202
0203 bool keepFirst = (nhit1 > nhit2 ||
0204 (nhit1 == nhit2 && track->track()->normalizedChi2() < track2->track()->normalizedChi2()));
0205
0206
0207 keepFirst |= (probFirst > minProb && probSecond <= minProb);
0208 keepFirst &= !(probFirst <= minProb && probSecond > minProb);
0209
0210 bool keepSecond = !keepFirst;
0211
0212
0213
0214 selected1[i] = selected1[i] && ((keepFirst && outputPreferCollection == 3) || outputPreferCollection == -1 ||
0215 outputPreferCollection == 1);
0216 track->setIsTrackerOnly(track->isTrackerOnly() &&
0217 ((keepFirst && trackerOnlyPreferCollection == 3) ||
0218 trackerOnlyPreferCollection == -1 || trackerOnlyPreferCollection == 1));
0219 track->setIsArbitratedEcalSeeded(track->isArbitratedEcalSeeded() &&
0220 ((keepFirst && arbitratedEcalSeededPreferCollection == 3) ||
0221 arbitratedEcalSeededPreferCollection == -1 ||
0222 arbitratedEcalSeededPreferCollection == 1));
0223 track->setIsArbitratedMerged(track->isArbitratedMerged() &&
0224 ((keepFirst && arbitratedMergedPreferCollection == 3) ||
0225 arbitratedMergedPreferCollection == -1 ||
0226 arbitratedMergedPreferCollection == 1));
0227 track->setIsArbitratedMergedEcalGeneral(track->isArbitratedMergedEcalGeneral() &&
0228 ((keepFirst && arbitratedMergedEcalGeneralPreferCollection == 3) ||
0229 arbitratedMergedEcalGeneralPreferCollection == -1 ||
0230 arbitratedMergedEcalGeneralPreferCollection == 1));
0231
0232 selected2[j] = selected2[j] && ((keepSecond && outputPreferCollection == 3) || outputPreferCollection == -1 ||
0233 outputPreferCollection == 2);
0234 track2->setIsTrackerOnly(track2->isTrackerOnly() &&
0235 ((keepSecond && trackerOnlyPreferCollection == 3) ||
0236 trackerOnlyPreferCollection == -1 || trackerOnlyPreferCollection == 2));
0237 track2->setIsArbitratedEcalSeeded(track2->isArbitratedEcalSeeded() &&
0238 ((keepSecond && arbitratedEcalSeededPreferCollection == 3) ||
0239 arbitratedEcalSeededPreferCollection == -1 ||
0240 arbitratedEcalSeededPreferCollection == 2));
0241 track2->setIsArbitratedMerged(track2->isArbitratedMerged() &&
0242 ((keepSecond && arbitratedMergedPreferCollection == 3) ||
0243 arbitratedMergedPreferCollection == -1 ||
0244 arbitratedMergedPreferCollection == 2));
0245 track2->setIsArbitratedMergedEcalGeneral(track2->isArbitratedMergedEcalGeneral() &&
0246 ((keepSecond && arbitratedMergedEcalGeneralPreferCollection == 3) ||
0247 arbitratedMergedEcalGeneralPreferCollection == -1 ||
0248 arbitratedMergedEcalGeneralPreferCollection == 2));
0249
0250 }
0251 }
0252 }
0253 }
0254
0255
0256
0257
0258
0259 if (!tC1.empty()) {
0260 i = 0;
0261 for (reco::ConversionTrackCollection::const_iterator track = tC1.begin(); track != tC1.end(); ++track, ++i) {
0262
0263 if (!selected1[i]) {
0264 continue;
0265 }
0266
0267 outputTrks->push_back(*track);
0268 }
0269 }
0270
0271
0272
0273 if (!tC2.empty()) {
0274 i = 0;
0275 for (reco::ConversionTrackCollection::const_iterator track = tC2.begin(); track != tC2.end(); ++track, ++i) {
0276
0277 if (!selected2[i]) {
0278 continue;
0279 }
0280
0281 outputTrks->push_back(*track);
0282 }
0283 }
0284
0285 e.put(std::move(outputTrks));
0286 return;
0287
0288 }