Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2021-02-14 13:34:17 

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             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         //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());
0193         //std::cout << " trk1 trk2 nhits1 nhits2 nover " << i << " " << j << " " << track->numberOfValidHits() << " "  << track2->numberOfValidHits() << " " << noverlap << " " << fi << " " << fj  <<std::endl;
0194         if ((noverlap - firstoverlap) > (std::min(nhit1, nhit2) - firstoverlap) * shareFrac &&
0195             (!checkCharge || track->track()->charge() * track2->track()->charge() > 0)) {
0196           //printf("overlapping tracks\n");
0197           //printf ("ndof1 = %5f, chisq1 = %5f, ndof2 = %5f, chisq2 = %5f\n",track->track()->ndof(),track->track()->chi2(),track2->track()->ndof(),track2->track()->chi2());
0198 
0199           double probFirst = ChiSquaredProbability(track->track()->chi2(), track->track()->ndof());
0200           double probSecond = ChiSquaredProbability(track2->track()->chi2(), track2->track()->ndof());
0201 
0202           //arbitrate by number of hits and reduced chisq
0203           bool keepFirst = (nhit1 > nhit2 ||
0204                             (nhit1 == nhit2 && track->track()->normalizedChi2() < track2->track()->normalizedChi2()));
0205 
0206           //override decision in case one track is radically worse quality than the other
0207           keepFirst |= (probFirst > minProb && probSecond <= minProb);
0208           keepFirst &= !(probFirst <= minProb && probSecond > minProb);
0209 
0210           bool keepSecond = !keepFirst;
0211 
0212           //set flags based on arbitration decision and precedence settings
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         }  //end got a duplicate
0251       }    //end track2 loop
0252     }      //end track loop
0253   }        //end more than 1 track
0254 
0255   //
0256   //  output selected tracks - if any
0257   //
0258 
0259   if (!tC1.empty()) {
0260     i = 0;
0261     for (reco::ConversionTrackCollection::const_iterator track = tC1.begin(); track != tC1.end(); ++track, ++i) {
0262       //don't store tracks rejected as duplicates
0263       if (!selected1[i]) {
0264         continue;
0265       }
0266       //fill the TrackCollection
0267       outputTrks->push_back(*track);
0268     }  //end faux loop over tracks
0269   }    //end more than 0 track
0270 
0271   //Fill the trajectories, etc. for 1st collection
0272 
0273   if (!tC2.empty()) {
0274     i = 0;
0275     for (reco::ConversionTrackCollection::const_iterator track = tC2.begin(); track != tC2.end(); ++track, ++i) {
0276       //don't store tracks rejected as duplicates
0277       if (!selected2[i]) {
0278         continue;
0279       }
0280       //fill the TrackCollection
0281       outputTrks->push_back(*track);
0282     }  //end faux loop over tracks
0283   }    //end more than 0 track
0284 
0285   e.put(std::move(outputTrks));
0286   return;
0287 
0288 }  //end produce
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!