EgammaPhotonProducers/src/ConversionTrackMerger.cc

0001 //
0002 // Package:         RecoTracker/FinalTrackSelectors
0003 // Class:           ConversionTrackMerger
0004 //
0005 // Description:     Merger for ConversionTracks, adapted from SimpleTrackListMerger
0006 //
0007 // Original Author: J.Bendavid
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   // retrieve producer name of input TrackCollection(s)
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 // Virtual destructor needed.
0061 ConversionTrackMerger::~ConversionTrackMerger() {}
0062
0063 // Functions that gets called by framework every event
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   // get Inputs
0078   // if 1 input list doesn't exist, make an empty list, issue a warning, and continue
0079   // this allows ConversionTrackMerger to be used as a cleaner only if handed just one list
0080   // if both input lists don't exist, will issue 2 warnings and generate an empty output collection
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     //std::cout << "1st collection " << trackProducer1 << " has "<< TC1->size() << " tracks" << std::endl ;
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     //std::cout << "2nd collection " << trackProducer2 << " has "<< TC2->size() << " tracks" << std::endl ;
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   // Step B: create empty output collection
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       //clear flags if preferCollection was set to 0
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         //clear flags if preferCollection was set to 0
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             for (unsigned jh = 0; jh < nh2; ++jh) {
0177               const TrackingRecHit* jt = jHits[jh];
0178               if (jt->isValid()) {
0179                 if (it->sharesInput(jt, TrackingRecHit::some)) {
0180                   noverlap++;
0181                   if (allowFirstHitShare && (ih == 0) && (jh == 0))
0182                     firstoverlap = 1;
0183                 }
0184               }
0185             }
0186           }
0187         }
0188         int nhit1 = track->track()->numberOfValidHits();
0189         int nhit2 = track2->track()->numberOfValidHits();
0190         //if (noverlap>0) printf("noverlap = %i, firstoverlap = %i, nhit1 = %i, nhit2 = %i, algo1 = %i, algo2 = %i, q1 = %i, q2 = %i\n",noverlap,firstoverlap,nhit1,nhit2,track->track()->algo(),track2->track()->algo(),track->track()->charge(),track2->track()->charge());
0191         //std::cout << " trk1 trk2 nhits1 nhits2 nover " << i << " " << j << " " << track->numberOfValidHits() << " "  << track2->numberOfValidHits() << " " << noverlap << " " << fi << " " << fj  <<std::endl;
0192         if ((noverlap - firstoverlap) > (std::min(nhit1, nhit2) - firstoverlap) * shareFrac &&
0193             (!checkCharge || track->track()->charge() * track2->track()->charge() > 0)) {
0194           //printf("overlapping tracks\n");
0195           //printf ("ndof1 = %5f, chisq1 = %5f, ndof2 = %5f, chisq2 = %5f\n",track->track()->ndof(),track->track()->chi2(),track2->track()->ndof(),track2->track()->chi2());
0196
0197           double probFirst = ChiSquaredProbability(track->track()->chi2(), track->track()->ndof());
0198           double probSecond = ChiSquaredProbability(track2->track()->chi2(), track2->track()->ndof());
0199
0200           //arbitrate by number of hits and reduced chisq
0201           bool keepFirst = (nhit1 > nhit2 ||
0202                             (nhit1 == nhit2 && track->track()->normalizedChi2() < track2->track()->normalizedChi2()));
0203
0204           //override decision in case one track is radically worse quality than the other
0205           keepFirst |= (probFirst > minProb && probSecond <= minProb);
0206           keepFirst &= !(probFirst <= minProb && probSecond > minProb);
0207
0208           bool keepSecond = !keepFirst;
0209
0210           //set flags based on arbitration decision and precedence settings
0211
0212           selected1[i] = selected1[i] && ((keepFirst && outputPreferCollection == 3) || outputPreferCollection == -1 ||
0213                                           outputPreferCollection == 1);
0214           track->setIsTrackerOnly(track->isTrackerOnly() &&
0215                                   ((keepFirst && trackerOnlyPreferCollection == 3) ||
0216                                    trackerOnlyPreferCollection == -1 || trackerOnlyPreferCollection == 1));
0217           track->setIsArbitratedEcalSeeded(track->isArbitratedEcalSeeded() &&
0218                                            ((keepFirst && arbitratedEcalSeededPreferCollection == 3) ||
0219                                             arbitratedEcalSeededPreferCollection == -1 ||
0220                                             arbitratedEcalSeededPreferCollection == 1));
0221           track->setIsArbitratedMerged(track->isArbitratedMerged() &&
0222                                        ((keepFirst && arbitratedMergedPreferCollection == 3) ||
0223                                         arbitratedMergedPreferCollection == -1 ||
0224                                         arbitratedMergedPreferCollection == 1));
0225           track->setIsArbitratedMergedEcalGeneral(track->isArbitratedMergedEcalGeneral() &&
0226                                                   ((keepFirst && arbitratedMergedEcalGeneralPreferCollection == 3) ||
0227                                                    arbitratedMergedEcalGeneralPreferCollection == -1 ||
0228                                                    arbitratedMergedEcalGeneralPreferCollection == 1));
0229
0230           selected2[j] = selected2[j] && ((keepSecond && outputPreferCollection == 3) || outputPreferCollection == -1 ||
0231                                           outputPreferCollection == 2);
0232           track2->setIsTrackerOnly(track2->isTrackerOnly() &&
0233                                    ((keepSecond && trackerOnlyPreferCollection == 3) ||
0234                                     trackerOnlyPreferCollection == -1 || trackerOnlyPreferCollection == 2));
0235           track2->setIsArbitratedEcalSeeded(track2->isArbitratedEcalSeeded() &&
0236                                             ((keepSecond && arbitratedEcalSeededPreferCollection == 3) ||
0237                                              arbitratedEcalSeededPreferCollection == -1 ||
0238                                              arbitratedEcalSeededPreferCollection == 2));
0239           track2->setIsArbitratedMerged(track2->isArbitratedMerged() &&
0240                                         ((keepSecond && arbitratedMergedPreferCollection == 3) ||
0241                                          arbitratedMergedPreferCollection == -1 ||
0242                                          arbitratedMergedPreferCollection == 2));
0243           track2->setIsArbitratedMergedEcalGeneral(track2->isArbitratedMergedEcalGeneral() &&
0244                                                    ((keepSecond && arbitratedMergedEcalGeneralPreferCollection == 3) ||
0245                                                     arbitratedMergedEcalGeneralPreferCollection == -1 ||
0246                                                     arbitratedMergedEcalGeneralPreferCollection == 2));
0247
0248         }  //end got a duplicate
0249       }    //end track2 loop
0250     }      //end track loop
0251   }        //end more than 1 track
0252
0253   //
0254   //  output selected tracks - if any
0255   //
0256
0257   if (!tC1.empty()) {
0258     i = 0;
0259     for (reco::ConversionTrackCollection::const_iterator track = tC1.begin(); track != tC1.end(); ++track, ++i) {
0260       //don't store tracks rejected as duplicates
0261       if (!selected1[i]) {
0262         continue;
0263       }
0264       //fill the TrackCollection
0265       outputTrks->push_back(*track);
0266     }  //end faux loop over tracks
0267   }    //end more than 0 track
0268
0269   //Fill the trajectories, etc. for 1st collection
0270
0271   if (!tC2.empty()) {
0272     i = 0;
0273     for (reco::ConversionTrackCollection::const_iterator track = tC2.begin(); track != tC2.end(); ++track, ++i) {
0274       //don't store tracks rejected as duplicates
0275       if (!selected2[i]) {
0276         continue;
0277       }
0278       //fill the TrackCollection
0279       outputTrks->push_back(*track);
0280     }  //end faux loop over tracks
0281   }    //end more than 0 track
0282
0283   e.put(std::move(outputTrks));
0284   return;
0285
0286 }  //end produce