Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:28:08 

0001 //
0002 // Package:         RecoTracker/FinalTrackSelectors
0003 // Class:           TrackListMerger
0004 //
0005 // Description:     Hit Dumper
0006 //
0007 // Original Author: Steve Wagner, stevew@pizero.colorado.edu
0008 // Created:         Sat Jan 14 22:00:00 UTC 2006
0009 //
0010 //
0011 
0012 #include "DataFormats/Common/interface/Handle.h"
0013 #include "DataFormats/Common/interface/ValueMap.h"
0014 #include "DataFormats/TrackReco/interface/Track.h"
0015 #include "DataFormats/TrackReco/interface/TrackBase.h"
0016 #include "DataFormats/TrackReco/interface/TrackExtra.h"
0017 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0018 #include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
0019 #include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h"
0020 #include "FWCore/Framework/interface/Event.h"
0021 #include "FWCore/Framework/interface/EventSetup.h"
0022 #include "FWCore/Framework/interface/stream/EDProducer.h"
0023 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0024 #include "FWCore/Utilities/interface/InputTag.h"
0025 #include "FWCore/Utilities/interface/thread_safety_macros.h"
0026 #include "FWCore/Utilities/interface/ESGetToken.h"
0027 #include "RecoTracker/FinalTrackSelectors/interface/TrackAlgoPriorityOrder.h"
0028 #include "RecoTracker/Record/interface/CkfComponentsRecord.h"
0029 #include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h"
0030 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0031 
0032 class dso_hidden TrackListMerger : public edm::stream::EDProducer<> {
0033 public:
0034   explicit TrackListMerger(const edm::ParameterSet& conf);
0035 
0036   ~TrackListMerger() override = default;
0037 
0038   void produce(edm::Event& e, const edm::EventSetup& c) override;
0039 
0040 private:
0041   void returnEmptyCollections(edm::Event& e);
0042 
0043   using MVACollection = std::vector<float>;
0044   using QualityMaskCollection = std::vector<unsigned char>;
0045 
0046   std::unique_ptr<reco::TrackCollection> outputTrks;
0047   std::unique_ptr<reco::TrackExtraCollection> outputTrkExtras;
0048   std::unique_ptr<TrackingRecHitCollection> outputTrkHits;
0049   std::unique_ptr<std::vector<Trajectory>> outputTrajs;
0050   std::unique_ptr<TrajTrackAssociationCollection> outputTTAss;
0051   std::unique_ptr<TrajectorySeedCollection> outputSeeds;
0052 
0053   reco::TrackRefProd refTrks;
0054   reco::TrackExtraRefProd refTrkExtras;
0055   TrackingRecHitRefProd refTrkHits;
0056   edm::RefProd<std::vector<Trajectory>> refTrajs;
0057   edm::RefProd<TrajectorySeedCollection> refTrajSeeds;
0058 
0059   bool copyExtras_;
0060   bool makeReKeyedSeeds_;
0061 
0062   edm::ESGetToken<TrackAlgoPriorityOrder, CkfComponentsRecord> priorityToken;
0063 
0064   struct TkEDGetTokenss {
0065     edm::InputTag tag;
0066     edm::EDGetTokenT<reco::TrackCollection> tk;
0067     edm::EDGetTokenT<std::vector<Trajectory>> traj;
0068     edm::EDGetTokenT<TrajTrackAssociationCollection> tass;
0069     edm::EDGetTokenT<edm::ValueMap<int>> tsel;
0070     edm::EDGetTokenT<edm::ValueMap<float>> tmva;
0071     TkEDGetTokenss() {}
0072     TkEDGetTokenss(const edm::InputTag& tag_,
0073                    edm::EDGetTokenT<reco::TrackCollection>&& tk_,
0074                    edm::EDGetTokenT<std::vector<Trajectory>>&& traj_,
0075                    edm::EDGetTokenT<TrajTrackAssociationCollection>&& tass_,
0076                    edm::EDGetTokenT<edm::ValueMap<int>>&& tsel_,
0077                    edm::EDGetTokenT<edm::ValueMap<float>>&& tmva_)
0078         : tag(tag_), tk(tk_), traj(traj_), tass(tass_), tsel(tsel_), tmva(tmva_) {}
0079   };
0080   TkEDGetTokenss edTokens(const edm::InputTag& tag, const edm::InputTag& seltag, const edm::InputTag& mvatag) {
0081     return TkEDGetTokenss(tag,
0082                           consumes<reco::TrackCollection>(tag),
0083                           consumes<std::vector<Trajectory>>(tag),
0084                           consumes<TrajTrackAssociationCollection>(tag),
0085                           consumes<edm::ValueMap<int>>(seltag),
0086                           consumes<edm::ValueMap<float>>(mvatag));
0087   }
0088   TkEDGetTokenss edTokens(const edm::InputTag& tag, const edm::InputTag& mvatag) {
0089     return TkEDGetTokenss(tag,
0090                           consumes<reco::TrackCollection>(tag),
0091                           consumes<std::vector<Trajectory>>(tag),
0092                           consumes<TrajTrackAssociationCollection>(tag),
0093                           edm::EDGetTokenT<edm::ValueMap<int>>(),
0094                           consumes<edm::ValueMap<float>>(mvatag));
0095   }
0096   std::vector<TkEDGetTokenss> trackProducers_;
0097 
0098   std::string priorityName_;
0099 
0100   double maxNormalizedChisq_;
0101   double minPT_;
0102   unsigned int minFound_;
0103   float epsilon_;
0104   float shareFrac_;
0105   float foundHitBonus_;
0106   float lostHitPenalty_;
0107   std::vector<double> indivShareFrac_;
0108 
0109   std::vector<std::vector<int>> listsToMerge_;
0110   std::vector<bool> promoteQuality_;
0111   std::vector<int> hasSelector_;
0112   bool copyMVA_;
0113 
0114   bool allowFirstHitShare_;
0115   reco::TrackBase::TrackQuality qualityToSet_;
0116   bool use_sharesInput_;
0117   bool trkQualMod_;
0118 };
0119 
0120 #include <memory>
0121 #include <string>
0122 #include <iostream>
0123 #include <cmath>
0124 #include <vector>
0125 
0126 #include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2DCollection.h"
0127 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2DCollection.h"
0128 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit1DCollection.h"
0129 #include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h"
0130 #include "DataFormats/TrackCandidate/interface/TrackCandidateCollection.h"
0131 #include "DataFormats/Common/interface/ValueMap.h"
0132 #include "FWCore/Framework/interface/ESHandle.h"
0133 
0134 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0135 
0136 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0137 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0138 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0139 
0140 //#include "DataFormats/TrackReco/src/classes.h"
0141 
0142 #include "TrackingTools/PatternTools/interface/ClusterRemovalRefSetter.h"
0143 
0144 #ifdef STAT_TSB
0145 #include <x86intrin.h>
0146 #endif
0147 
0148 namespace {
0149 #ifdef STAT_TSB
0150   inline volatile unsigned long long rdtsc() { return __rdtsc(); }
0151 
0152   struct StatCount {
0153     float maxDP = 0.;
0154     float maxDE = 0.;
0155     unsigned long long st;
0156     long long totBegin = 0;
0157     long long totPre = 0;
0158     long long totEnd = 0;
0159     unsigned long long timeNo;  // no-overlap
0160     unsigned long long timeOv;  // overlap
0161     void begin(int tt) { totBegin += tt; }
0162     void start() { st = rdtsc(); }
0163     void noOverlap() { timeNo += (rdtsc() - st); }
0164     void overlap() { timeOv += (rdtsc() - st); }
0165     void pre(int tt) { totPre += tt; }
0166     void end(int tt) { totEnd += tt; }
0167     void de(float d) {
0168       if (d > maxDE)
0169         maxDE = d;
0170     }
0171     void dp(float d) {
0172       if (d > maxDP)
0173         maxDP = d;
0174     }
0175 
0176     void print() const {
0177       std::cout << "TrackListMerger stat\nBegin/Pre/End/maxDPhi/maxDEta/Overlap/NoOverlap " << totBegin << '/' << totPre
0178                 << '/' << totEnd << '/' << maxDP << '/' << maxDE << '/' << timeOv / 1000 << '/' << timeNo / 1000
0179                 << std::endl;
0180     }
0181     StatCount() {}
0182     ~StatCount() { print(); }
0183   };
0184   StatCount statCount;
0185 
0186 #else
0187   struct StatCount {
0188     void begin(int) {}
0189     void pre(int) {}
0190     void end(int) {}
0191     void start() {}
0192     void noOverlap() {}
0193     void overlap() {}
0194     void de(float) {}
0195     void dp(float) {}
0196   };
0197   CMS_THREAD_SAFE StatCount statCount;
0198 #endif
0199 
0200 }  // namespace
0201 
0202 namespace {
0203   edm::ProductID clusterProductB(const TrackingRecHit* hit) {
0204     return reinterpret_cast<const BaseTrackerRecHit*>(hit)->firstClusterRef().id();
0205   }
0206 }  // namespace
0207 
0208 TrackListMerger::TrackListMerger(edm::ParameterSet const& conf) {
0209   copyExtras_ = conf.getUntrackedParameter<bool>("copyExtras", true);
0210   priorityName_ = conf.getParameter<std::string>("trackAlgoPriorityOrder");
0211 
0212   std::vector<edm::InputTag> trackProducerTags(conf.getParameter<std::vector<edm::InputTag>>("TrackProducers"));
0213   //which of these do I need to turn into vectors?
0214   maxNormalizedChisq_ = conf.getParameter<double>("MaxNormalizedChisq");
0215   minPT_ = conf.getParameter<double>("MinPT");
0216   minFound_ = (unsigned int)conf.getParameter<int>("MinFound");
0217   epsilon_ = conf.getParameter<double>("Epsilon");
0218   shareFrac_ = conf.getParameter<double>("ShareFrac");
0219   allowFirstHitShare_ = conf.getParameter<bool>("allowFirstHitShare");
0220   foundHitBonus_ = conf.getParameter<double>("FoundHitBonus");
0221   lostHitPenalty_ = conf.getParameter<double>("LostHitPenalty");
0222   indivShareFrac_ = conf.getParameter<std::vector<double>>("indivShareFrac");
0223   std::string qualityStr = conf.getParameter<std::string>("newQuality");
0224   priorityToken = esConsumes<TrackAlgoPriorityOrder, CkfComponentsRecord>(edm::ESInputTag("", priorityName_));
0225 
0226   if (!qualityStr.empty()) {
0227     qualityToSet_ = reco::TrackBase::qualityByName(conf.getParameter<std::string>("newQuality"));
0228   } else
0229     qualityToSet_ = reco::TrackBase::undefQuality;
0230 
0231   use_sharesInput_ = true;
0232   if (epsilon_ > 0.0)
0233     use_sharesInput_ = false;
0234 
0235   edm::VParameterSet setsToMerge = conf.getParameter<edm::VParameterSet>("setsToMerge");
0236 
0237   for (unsigned int i = 0; i < setsToMerge.size(); i++) {
0238     listsToMerge_.push_back(setsToMerge[i].getParameter<std::vector<int>>("tLists"));
0239     promoteQuality_.push_back(setsToMerge[i].getParameter<bool>("pQual"));
0240   }
0241   hasSelector_ = conf.getParameter<std::vector<int>>("hasSelector");
0242   copyMVA_ = conf.getParameter<bool>("copyMVA");
0243 
0244   std::vector<edm::InputTag> selectors(conf.getParameter<std::vector<edm::InputTag>>("selectedTrackQuals"));
0245   std::vector<edm::InputTag> mvaStores;
0246   if (conf.exists("mvaValueTags")) {
0247     mvaStores = conf.getParameter<std::vector<edm::InputTag>>("mvaValueTags");
0248   } else {
0249     for (int i = 0; i < (int)selectors.size(); i++) {
0250       edm::InputTag ntag(selectors[i].label(), "MVAVals");
0251       mvaStores.push_back(ntag);
0252     }
0253   }
0254   unsigned int numTrkColl = trackProducerTags.size();
0255   if (numTrkColl != hasSelector_.size() || numTrkColl != selectors.size()) {
0256     throw cms::Exception("Inconsistent size") << "need same number of track collections and selectors";
0257   }
0258   if (numTrkColl != hasSelector_.size() || numTrkColl != mvaStores.size()) {
0259     throw cms::Exception("Inconsistent size") << "need same number of track collections and MVA stores";
0260   }
0261   for (unsigned int i = indivShareFrac_.size(); i < numTrkColl; i++) {
0262     //      edm::LogWarning("TrackListMerger") << "No indivShareFrac for " << trackProducersTags <<". Using default value of 1";
0263     indivShareFrac_.push_back(1.0);
0264   }
0265 
0266   trkQualMod_ = conf.getParameter<bool>("writeOnlyTrkQuals");
0267   if (trkQualMod_) {
0268     bool ok = true;
0269     for (unsigned int i = 1; i < numTrkColl; i++) {
0270       if (!(trackProducerTags[i] == trackProducerTags[0]))
0271         ok = false;
0272     }
0273     if (!ok) {
0274       throw cms::Exception("Bad input") << "to use writeOnlyTrkQuals=True all input InputTags must be the same";
0275     }
0276     produces<edm::ValueMap<int>>();
0277     produces<QualityMaskCollection>("QualityMasks");
0278   } else {
0279     makeReKeyedSeeds_ = conf.getUntrackedParameter<bool>("makeReKeyedSeeds", false);
0280     if (makeReKeyedSeeds_) {
0281       copyExtras_ = true;
0282       produces<TrajectorySeedCollection>();
0283     }
0284 
0285     if (copyExtras_) {
0286       produces<reco::TrackExtraCollection>();
0287       produces<TrackingRecHitCollection>();
0288     }
0289 
0290     // TrackCollection refers to TrackingRechit and TrackExtra
0291     // collections, need to declare its production after them to work
0292     // around a rare race condition in framework scheduling
0293     produces<reco::TrackCollection>();
0294 
0295     produces<std::vector<Trajectory>>();
0296     produces<TrajTrackAssociationCollection>();
0297   }
0298   produces<edm::ValueMap<float>>("MVAVals");
0299   produces<MVACollection>("MVAValues");
0300 
0301   // Do all the consumes
0302   trackProducers_.resize(numTrkColl);
0303   for (unsigned int i = 0; i < numTrkColl; ++i) {
0304     trackProducers_[i] = hasSelector_[i] > 0 ? edTokens(trackProducerTags[i], selectors[i], mvaStores[i])
0305                                              : edTokens(trackProducerTags[i], mvaStores[i]);
0306   }
0307 }
0308 
0309 // Functions that gets called by framework every event
0310 void TrackListMerger::produce(edm::Event& e, const edm::EventSetup& es) {
0311   // extract tracker geometry
0312   //
0313   //edm::ESHandle<TrackerGeometry> theG;
0314   //es.get<TrackerDigiGeometryRecord>().get(theG);
0315 
0316   //    using namespace reco;
0317 
0318   edm::ESHandle<TrackAlgoPriorityOrder> priorityH = es.getHandle(priorityToken);
0319   auto const& trackAlgoPriorityOrder = *priorityH;
0320 
0321   // get Inputs
0322   // if 1 input list doesn't exist, make an empty list, issue a warning, and continue
0323   // this allows TrackListMerger to be used as a cleaner only if handed just one list
0324   // if both input lists don't exist, will issue 2 warnings and generate an empty output collection
0325   //
0326   static const reco::TrackCollection s_empty;
0327 
0328   std::vector<const reco::TrackCollection*> trackColls;
0329   std::vector<edm::Handle<reco::TrackCollection>> trackHandles(trackProducers_.size());
0330   for (unsigned int i = 0; i < trackProducers_.size(); i++) {
0331     trackColls.push_back(nullptr);
0332     //edm::Handle<reco::TrackCollection> trackColl;
0333     e.getByToken(trackProducers_[i].tk, trackHandles[i]);
0334     if (trackHandles[i].isValid()) {
0335       trackColls[i] = trackHandles[i].product();
0336     } else {
0337       edm::LogWarning("TrackListMerger") << "TrackCollection " << trackProducers_[i].tag << " not found";
0338       trackColls[i] = &s_empty;
0339     }
0340   }
0341 
0342   if (trackColls.empty()) {
0343     // output empty collections and early return
0344     this->returnEmptyCollections(e);
0345     return;
0346   }
0347 
0348   unsigned int collsSize = trackColls.size();
0349   unsigned int rSize = 0;
0350   unsigned int trackCollSizes[collsSize];
0351   unsigned int trackCollFirsts[collsSize];
0352   for (unsigned int i = 0; i != collsSize; i++) {
0353     trackCollSizes[i] = trackColls[i]->size();
0354     trackCollFirsts[i] = rSize;
0355     rSize += trackCollSizes[i];
0356   }
0357 
0358   if (rSize == 0) {
0359     // output empty collections and early return
0360     this->returnEmptyCollections(e);
0361     return;
0362   }
0363 
0364   statCount.begin(rSize);
0365 
0366   //
0367   //  quality cuts first
0368   //
0369   int i = -1;
0370 
0371   int selected[rSize];
0372   int indexG[rSize];
0373   bool trkUpdated[rSize];
0374   int trackCollNum[rSize];
0375   int trackQuals[rSize];
0376   float trackMVAs[rSize];
0377   reco::TrackBase::TrackAlgorithm oriAlgo[rSize];
0378   std::vector<reco::TrackBase::AlgoMask> algoMask(rSize);
0379   for (unsigned int j = 0; j < rSize; j++) {
0380     indexG[j] = -1;
0381     selected[j] = 1;
0382     trkUpdated[j] = false;
0383     trackCollNum[j] = 0;
0384     trackQuals[j] = 0;
0385     trackMVAs[j] = -998.0;
0386     oriAlgo[j] = reco::TrackBase::undefAlgorithm;
0387   }
0388 
0389   int ngood = 0;
0390   for (unsigned int j = 0; j != collsSize; j++) {
0391     const reco::TrackCollection* tC1 = trackColls[j];
0392 
0393     edm::Handle<edm::ValueMap<int>> trackSelColl;
0394     edm::Handle<edm::ValueMap<float>> trackMVAStore;
0395     if (copyMVA_)
0396       e.getByToken(trackProducers_[j].tmva, trackMVAStore);
0397     if (hasSelector_[j] > 0) {
0398       e.getByToken(trackProducers_[j].tsel, trackSelColl);
0399     }
0400 
0401     if (!tC1->empty()) {
0402       unsigned int iC = 0;
0403       for (reco::TrackCollection::const_iterator track = tC1->begin(); track != tC1->end(); track++) {
0404         i++;
0405         trackCollNum[i] = j;
0406         trackQuals[i] = track->qualityMask();
0407         oriAlgo[i] = track->originalAlgo();
0408         algoMask[i] = track->algoMask();
0409 
0410         reco::TrackRef trkRef = reco::TrackRef(trackHandles[j], iC);
0411         if (copyMVA_)
0412           if ((*trackMVAStore).contains(trkRef.id()))
0413             trackMVAs[i] = (*trackMVAStore)[trkRef];
0414         if (hasSelector_[j] > 0) {
0415           int qual = (*trackSelColl)[trkRef];
0416           if (qual < 0) {
0417             selected[i] = 0;
0418             iC++;
0419             continue;
0420           } else {
0421             trackQuals[i] = qual;
0422           }
0423         }
0424         iC++;
0425         selected[i] = trackQuals[i] + 10;  //10 is magic number used throughout...
0426         if ((short unsigned)track->ndof() < 1) {
0427           selected[i] = 0;
0428           continue;
0429         }
0430         if (track->normalizedChi2() > maxNormalizedChisq_) {
0431           selected[i] = 0;
0432           continue;
0433         }
0434         if (track->found() < minFound_) {
0435           selected[i] = 0;
0436           continue;
0437         }
0438         if (track->pt() < minPT_) {
0439           selected[i] = 0;
0440           continue;
0441         }
0442         // good!
0443         indexG[i] = ngood++;
0444         //if ( beVerb) std::cout << "inverb " << track->pt() << " " << selected[i] << std::endl;
0445       }  //end loop over tracks
0446     }    //end more than 0 track
0447   }      // loop over trackcolls
0448 
0449   statCount.pre(ngood);
0450 
0451   //cache the id and rechits of valid hits
0452   typedef std::pair<unsigned int, const TrackingRecHit*> IHit;
0453   std::vector<std::vector<IHit>> rh1(ngood);  // "not an array" of vectors!
0454   //const TrackingRecHit*  fh1[ngood];  // first hit...
0455   reco::TrackBase::TrackAlgorithm algo[std::max(1, ngood)];
0456   float score[std::max(1, ngood)];
0457 
0458   for (unsigned int j = 0; j < rSize; j++) {
0459     if (selected[j] == 0)
0460       continue;
0461     int i = indexG[j];
0462     assert(i >= 0);
0463     unsigned int collNum = trackCollNum[j];
0464     unsigned int trackNum = j - trackCollFirsts[collNum];
0465     const reco::Track* track = &((trackColls[collNum])->at(trackNum));
0466 
0467     algo[i] = track->algo();
0468     int validHits = track->numberOfValidHits();
0469     int validPixelHits = track->hitPattern().numberOfValidPixelHits();
0470     int lostHits = track->numberOfLostHits();
0471     score[i] =
0472         foundHitBonus_ * validPixelHits + foundHitBonus_ * validHits - lostHitPenalty_ * lostHits - track->chi2();
0473 
0474     rh1[i].reserve(validHits);
0475     auto compById = [](IHit const& h1, IHit const& h2) { return h1.first < h2.first; };
0476     for (trackingRecHit_iterator it = track->recHitsBegin(); it != track->recHitsEnd(); ++it) {
0477       const TrackingRecHit* hit = (*it);
0478       unsigned int id = hit->rawId();
0479       if (hit->geographicalId().subdetId() > 2)
0480         id &= (~3);  // mask mono/stereo in strips...
0481       if LIKELY (hit->isValid()) {
0482         rh1[i].emplace_back(id, hit);
0483         std::push_heap(rh1[i].begin(), rh1[i].end(), compById);
0484       }
0485     }
0486     std::sort_heap(rh1[i].begin(), rh1[i].end(), compById);
0487   }
0488 
0489   //DL here
0490   if LIKELY (ngood > 1 && collsSize > 1)
0491     for (unsigned int ltm = 0; ltm < listsToMerge_.size(); ltm++) {
0492       int saveSelected[rSize];
0493       bool notActive[collsSize];
0494       for (unsigned int cn = 0; cn != collsSize; ++cn)
0495         notActive[cn] = find(listsToMerge_[ltm].begin(), listsToMerge_[ltm].end(), cn) == listsToMerge_[ltm].end();
0496 
0497       for (unsigned int i = 0; i < rSize; i++)
0498         saveSelected[i] = selected[i];
0499 
0500       //DL protect against 0 tracks?
0501       for (unsigned int i = 0; i < rSize - 1; i++) {
0502         if (selected[i] == 0)
0503           continue;
0504         unsigned int collNum = trackCollNum[i];
0505 
0506         //check that this track is in one of the lists for this iteration
0507         if (notActive[collNum])
0508           continue;
0509 
0510         int k1 = indexG[i];
0511         unsigned int nh1 = rh1[k1].size();
0512         int qualityMaskT1 = trackQuals[i];
0513 
0514         int nhit1 = nh1;  // validHits[k1];
0515         float score1 = score[k1];
0516 
0517         // start at next collection
0518         for (unsigned int j = i + 1; j < rSize; j++) {
0519           if (selected[j] == 0)
0520             continue;
0521           unsigned int collNum2 = trackCollNum[j];
0522           if ((collNum == collNum2) && indivShareFrac_[collNum] > 0.99)
0523             continue;
0524           //check that this track is in one of the lists for this iteration
0525           if (notActive[collNum2])
0526             continue;
0527 
0528           int k2 = indexG[j];
0529 
0530           int newQualityMask = -9;  //avoid resetting quality mask if not desired 10+ -9 =1
0531           if (promoteQuality_[ltm]) {
0532             int maskT1 = saveSelected[i] > 1 ? saveSelected[i] - 10 : qualityMaskT1;
0533             int maskT2 = saveSelected[j] > 1 ? saveSelected[j] - 10 : trackQuals[j];
0534             newQualityMask = (maskT1 | maskT2);  // take OR of trackQuality
0535           }
0536           unsigned int nh2 = rh1[k2].size();
0537           int nhit2 = nh2;
0538 
0539           auto share = use_sharesInput_ ? [](const TrackingRecHit* it, const TrackingRecHit* jt, float) -> bool {
0540             return it->sharesInput(jt, TrackingRecHit::some);
0541           }
0542           : [](const TrackingRecHit* it, const TrackingRecHit* jt, float eps) -> bool {
0543               float delta = std::abs(it->localPosition().x() - jt->localPosition().x());
0544               return (it->geographicalId() == jt->geographicalId()) && (delta < eps);
0545             };
0546 
0547           statCount.start();
0548 
0549           //loop over rechits
0550           int noverlap = 0;
0551           int firstoverlap = 0;
0552           // check first hit  (should use REAL first hit?)
0553           if UNLIKELY (allowFirstHitShare_ && rh1[k1][0].first == rh1[k2][0].first) {
0554             const TrackingRecHit* it = rh1[k1][0].second;
0555             const TrackingRecHit* jt = rh1[k2][0].second;
0556             if (share(it, jt, epsilon_))
0557               firstoverlap = 1;
0558           }
0559 
0560           // exploit sorting
0561           unsigned int jh = 0;
0562           unsigned int ih = 0;
0563           while (ih != nh1 && jh != nh2) {
0564             // break if not enough to go...
0565             // if ( nprecut-noverlap+firstoverlap > int(nh1-ih)) break;
0566             // if ( nprecut-noverlap+firstoverlap > int(nh2-jh)) break;
0567             auto const id1 = rh1[k1][ih].first;
0568             auto const id2 = rh1[k2][jh].first;
0569             if (id1 < id2)
0570               ++ih;
0571             else if (id2 < id1)
0572               ++jh;
0573             else {
0574               // in case of split-hit do full conbinatorics
0575               auto li = ih;
0576               while ((++li) != nh1 && id1 == rh1[k1][li].first) {
0577               }
0578               auto lj = jh;
0579               while ((++lj) != nh2 && id2 == rh1[k2][lj].first) {
0580               }
0581               for (auto ii = ih; ii != li; ++ii)
0582                 for (auto jj = jh; jj != lj; ++jj) {
0583                   const TrackingRecHit* it = rh1[k1][ii].second;
0584                   const TrackingRecHit* jt = rh1[k2][jj].second;
0585                   if (share(it, jt, epsilon_))
0586                     noverlap++;
0587                 }
0588               jh = lj;
0589               ih = li;
0590             }  // equal ids
0591 
0592           }  //loop over ih & jh
0593 
0594           bool dupfound =
0595               (collNum != collNum2)
0596                   ? (noverlap - firstoverlap) > (std::min(nhit1, nhit2) - firstoverlap) * shareFrac_
0597                   : (noverlap - firstoverlap) > (std::min(nhit1, nhit2) - firstoverlap) * indivShareFrac_[collNum];
0598 
0599           auto seti = [&](unsigned int ii, unsigned int jj) {
0600             selected[jj] = 0;
0601             selected[ii] = 10 + newQualityMask;  // add 10 to avoid the case where mask = 1
0602             trkUpdated[ii] = true;
0603             if (trackAlgoPriorityOrder.priority(oriAlgo[jj]) < trackAlgoPriorityOrder.priority(oriAlgo[ii]))
0604               oriAlgo[ii] = oriAlgo[jj];
0605             algoMask[ii] |= algoMask[jj];
0606             algoMask[jj] = algoMask[ii];  // in case we keep discarded
0607           };
0608 
0609           if (dupfound) {
0610             float score2 = score[k2];
0611             constexpr float almostSame =
0612                 0.01f;  // difference rather than ratio due to possible negative values for score
0613             if (score1 - score2 > almostSame) {
0614               seti(i, j);
0615             } else if (score2 - score1 > almostSame) {
0616               seti(j, i);
0617             } else {
0618               // If tracks from both iterations are virtually identical, choose the one with the best quality or with lower algo
0619               if ((trackQuals[j] &
0620                    (1 << reco::TrackBase::loose | 1 << reco::TrackBase::tight | 1 << reco::TrackBase::highPurity)) ==
0621                   (trackQuals[i] &
0622                    (1 << reco::TrackBase::loose | 1 << reco::TrackBase::tight | 1 << reco::TrackBase::highPurity))) {
0623                 //same quality, pick earlier algo
0624                 if (trackAlgoPriorityOrder.priority(algo[k1]) <= trackAlgoPriorityOrder.priority(algo[k2])) {
0625                   seti(i, j);
0626                 } else {
0627                   seti(j, i);
0628                 }
0629               } else if ((trackQuals[j] & (1 << reco::TrackBase::loose | 1 << reco::TrackBase::tight |
0630                                            1 << reco::TrackBase::highPurity)) <
0631                          (trackQuals[i] & (1 << reco::TrackBase::loose | 1 << reco::TrackBase::tight |
0632                                            1 << reco::TrackBase::highPurity))) {
0633                 seti(i, j);
0634               } else {
0635                 seti(j, i);
0636               }
0637             }  //end fi < fj
0638             statCount.overlap();
0639             /*
0640         if (at0[k1]&&at0[k2]) {
0641           statCount.dp(dphi);
0642           if (dz<1.f) statCount.de(deta);
0643         }
0644         */
0645           }  //end got a duplicate
0646           else {
0647             statCount.noOverlap();
0648           }
0649           //stop if the ith track is now unselected
0650           if (selected[i] == 0)
0651             break;
0652         }  //end track2 loop
0653       }    //end track loop
0654     }      //end loop over track list sets
0655 
0656   auto vmMVA = std::make_unique<edm::ValueMap<float>>();
0657   edm::ValueMap<float>::Filler fillerMVA(*vmMVA);
0658 
0659   // special case - if just doing the trkquals
0660   if (trkQualMod_) {
0661     unsigned int tSize = trackColls[0]->size();
0662     auto vm = std::make_unique<edm::ValueMap<int>>();
0663     edm::ValueMap<int>::Filler filler(*vm);
0664 
0665     std::vector<int> finalQuals(tSize, -1);  //default is unselected
0666     for (unsigned int i = 0; i < rSize; i++) {
0667       unsigned int tNum = i % tSize;
0668 
0669       if (selected[i] > 1) {
0670         finalQuals[tNum] = selected[i] - 10;
0671         if (trkUpdated[i])
0672           finalQuals[tNum] = (finalQuals[tNum] | (1 << qualityToSet_));
0673       }
0674       if (selected[i] == 1)
0675         finalQuals[tNum] = trackQuals[i];
0676     }
0677 
0678     filler.insert(trackHandles[0], finalQuals.begin(), finalQuals.end());
0679     filler.fill();
0680 
0681     e.put(std::move(vm));
0682     for (auto& q : finalQuals)
0683       q = std::max(q, 0);
0684     auto quals = std::make_unique<QualityMaskCollection>(finalQuals.begin(), finalQuals.end());
0685     e.put(std::move(quals), "QualityMasks");
0686 
0687     std::vector<float> mvaVec(tSize, -99);
0688 
0689     for (unsigned int i = 0; i < rSize; i++) {
0690       unsigned int tNum = i % tSize;
0691       mvaVec[tNum] = trackMVAs[tNum];
0692     }
0693 
0694     fillerMVA.insert(trackHandles[0], mvaVec.begin(), mvaVec.end());
0695     fillerMVA.fill();
0696     if (copyMVA_) {
0697       e.put(std::move(vmMVA), "MVAVals");
0698       auto mvas = std::make_unique<MVACollection>(mvaVec.begin(), mvaVec.end());
0699       e.put(std::move(mvas), "MVAValues");
0700     }
0701     return;
0702   }
0703 
0704   //
0705   //  output selected tracks - if any
0706   //
0707 
0708   std::vector<reco::TrackRef> trackRefs(rSize);
0709   std::vector<edm::RefToBase<TrajectorySeed>> seedsRefs(rSize);
0710 
0711   unsigned int nToWrite = 0;
0712   for (unsigned int i = 0; i < rSize; i++)
0713     if (selected[i] != 0)
0714       nToWrite++;
0715 
0716   std::vector<float> mvaVec;
0717 
0718   outputTrks = std::make_unique<reco::TrackCollection>();
0719   outputTrks->reserve(nToWrite);
0720   refTrks = e.getRefBeforePut<reco::TrackCollection>();
0721 
0722   if (copyExtras_) {
0723     outputTrkExtras = std::make_unique<reco::TrackExtraCollection>();
0724     outputTrkExtras->reserve(nToWrite);
0725     refTrkExtras = e.getRefBeforePut<reco::TrackExtraCollection>();
0726     outputTrkHits = std::make_unique<TrackingRecHitCollection>();
0727     outputTrkHits->reserve(nToWrite * 25);
0728     refTrkHits = e.getRefBeforePut<TrackingRecHitCollection>();
0729     if (makeReKeyedSeeds_) {
0730       outputSeeds = std::make_unique<TrajectorySeedCollection>();
0731       outputSeeds->reserve(nToWrite);
0732       refTrajSeeds = e.getRefBeforePut<TrajectorySeedCollection>();
0733     }
0734   }
0735 
0736   outputTrajs = std::make_unique<std::vector<Trajectory>>();
0737   outputTrajs->reserve(rSize);
0738 
0739   for (unsigned int i = 0; i < rSize; i++) {
0740     if (selected[i] == 0) {
0741       trackRefs[i] = reco::TrackRef();
0742       continue;
0743     }
0744 
0745     unsigned int collNum = trackCollNum[i];
0746     unsigned int trackNum = i - trackCollFirsts[collNum];
0747     const reco::Track* track = &((trackColls[collNum])->at(trackNum));
0748     outputTrks->push_back(reco::Track(*track));
0749     mvaVec.push_back(trackMVAs[i]);
0750     if (selected[i] > 1) {
0751       outputTrks->back().setQualityMask(selected[i] - 10);
0752       if (trkUpdated[i])
0753         outputTrks->back().setQuality(qualityToSet_);
0754     }
0755     //might duplicate things, but doesnt hurt
0756     if (selected[i] == 1)
0757       outputTrks->back().setQualityMask(trackQuals[i]);
0758     outputTrks->back().setOriginalAlgorithm(oriAlgo[i]);
0759     outputTrks->back().setAlgoMask(algoMask[i]);
0760 
0761     // if ( beVerb ) std::cout << "selected " << outputTrks->back().pt() << " " << outputTrks->back().qualityMask() << " " << selected[i] << std::endl;
0762 
0763     //fill the TrackCollection
0764     if (copyExtras_) {
0765       edm::RefToBase<TrajectorySeed> origSeedRef = track->seedRef();
0766       //creating a seed with rekeyed clusters if required
0767       if (makeReKeyedSeeds_) {
0768         bool doRekeyOnThisSeed = false;
0769 
0770         edm::InputTag clusterRemovalInfos("");
0771         //grab on of the hits of the seed
0772         if (origSeedRef->nHits() != 0) {
0773           TrackingRecHit const& hit = *origSeedRef->recHits().begin();
0774           if (hit.isValid()) {
0775             edm::ProductID pID = clusterProductB(&hit);
0776             // the cluster collection either produced a removalInfo or mot
0777             //get the clusterremoval info from the provenance: will rekey if this is found
0778             edm::Handle<reco::ClusterRemovalInfo> CRIh;
0779             edm::StableProvenance const& prov = e.getStableProvenance(pID);
0780             clusterRemovalInfos = edm::InputTag(prov.moduleLabel(), prov.productInstanceName(), prov.processName());
0781             doRekeyOnThisSeed = e.getByLabel(clusterRemovalInfos, CRIh);
0782           }  //valid hit
0783         }    //nhit!=0
0784 
0785         if (doRekeyOnThisSeed && !(clusterRemovalInfos == edm::InputTag(""))) {
0786           ClusterRemovalRefSetter refSetter(e, clusterRemovalInfos);
0787           TrajectorySeed::RecHitContainer newRecHitContainer;
0788           newRecHitContainer.reserve(origSeedRef->nHits());
0789           for (auto const& recHit : origSeedRef->recHits()) {
0790             newRecHitContainer.push_back(recHit);
0791             refSetter.reKey(&newRecHitContainer.back());
0792           }
0793           outputSeeds->push_back(
0794               TrajectorySeed(origSeedRef->startingState(), newRecHitContainer, origSeedRef->direction()));
0795         }
0796         //doRekeyOnThisSeed=true
0797         else {
0798           //just copy the one we had before
0799           outputSeeds->push_back(TrajectorySeed(*origSeedRef));
0800         }
0801         edm::Ref<TrajectorySeedCollection> pureRef(refTrajSeeds, outputSeeds->size() - 1);
0802         origSeedRef = edm::RefToBase<TrajectorySeed>(pureRef);
0803       }  //creating a new seed and rekeying it rechit clusters.
0804 
0805       // Fill TrackExtra collection
0806       outputTrkExtras->push_back(reco::TrackExtra(track->outerPosition(),
0807                                                   track->outerMomentum(),
0808                                                   track->outerOk(),
0809                                                   track->innerPosition(),
0810                                                   track->innerMomentum(),
0811                                                   track->innerOk(),
0812                                                   track->outerStateCovariance(),
0813                                                   track->outerDetId(),
0814                                                   track->innerStateCovariance(),
0815                                                   track->innerDetId(),
0816                                                   track->seedDirection(),
0817                                                   origSeedRef));
0818       seedsRefs[i] = origSeedRef;
0819       outputTrks->back().setExtra(reco::TrackExtraRef(refTrkExtras, outputTrkExtras->size() - 1));
0820       reco::TrackExtra& tx = outputTrkExtras->back();
0821       tx.setResiduals(track->residuals());
0822 
0823       // fill TrackingRecHits
0824       unsigned nh1 = track->recHitsSize();
0825       tx.setHits(refTrkHits, outputTrkHits->size(), nh1);
0826       tx.setTrajParams(track->extra()->trajParams(), track->extra()->chi2sX5());
0827       assert(tx.trajParams().size() == tx.recHitsSize());
0828       for (auto hh = track->recHitsBegin(), eh = track->recHitsEnd(); hh != eh; ++hh) {
0829         outputTrkHits->push_back((*hh)->clone());
0830       }
0831     }
0832     trackRefs[i] = reco::TrackRef(refTrks, outputTrks->size() - 1);
0833 
0834   }  //end faux loop over tracks
0835 
0836   //Fill the trajectories, etc. for 1st collection
0837   refTrajs = e.getRefBeforePut<std::vector<Trajectory>>();
0838 
0839   outputTTAss = std::make_unique<TrajTrackAssociationCollection>(refTrajs, refTrks);
0840 
0841   for (unsigned int ti = 0; ti < trackColls.size(); ti++) {
0842     edm::Handle<std::vector<Trajectory>> hTraj1;
0843     edm::Handle<TrajTrackAssociationCollection> hTTAss1;
0844     e.getByToken(trackProducers_[ti].traj, hTraj1);
0845     e.getByToken(trackProducers_[ti].tass, hTTAss1);
0846 
0847     if (hTraj1.failedToGet() || hTTAss1.failedToGet())
0848       continue;
0849 
0850     for (size_t i = 0, n = hTraj1->size(); i < n; ++i) {
0851       edm::Ref<std::vector<Trajectory>> trajRef(hTraj1, i);
0852       TrajTrackAssociationCollection::const_iterator match = hTTAss1->find(trajRef);
0853       if (match != hTTAss1->end()) {
0854         const edm::Ref<reco::TrackCollection>& trkRef = match->val;
0855         uint32_t oldKey = trackCollFirsts[ti] + static_cast<uint32_t>(trkRef.key());
0856         if (trackRefs[oldKey].isNonnull()) {
0857           outputTrajs->push_back(*trajRef);
0858           //if making extras and the seeds at the same time, change the seed ref on the trajectory
0859           if (copyExtras_ && makeReKeyedSeeds_)
0860             outputTrajs->back().setSeedRef(seedsRefs[oldKey]);
0861           outputTTAss->insert(edm::Ref<std::vector<Trajectory>>(refTrajs, outputTrajs->size() - 1), trackRefs[oldKey]);
0862         }
0863       }
0864     }
0865   }
0866 
0867   statCount.end(outputTrks->size());
0868 
0869   edm::ProductID nPID = refTrks.id();
0870   edm::TestHandle<reco::TrackCollection> outHandle(outputTrks.get(), nPID);
0871   fillerMVA.insert(outHandle, mvaVec.begin(), mvaVec.end());
0872   fillerMVA.fill();
0873 
0874   e.put(std::move(outputTrks));
0875   if (copyMVA_) {
0876     e.put(std::move(vmMVA), "MVAVals");
0877     auto mvas = std::make_unique<MVACollection>(mvaVec.begin(), mvaVec.end());
0878     e.put(std::move(mvas), "MVAValues");
0879   }
0880   if (copyExtras_) {
0881     e.put(std::move(outputTrkExtras));
0882     e.put(std::move(outputTrkHits));
0883     if (makeReKeyedSeeds_)
0884       e.put(std::move(outputSeeds));
0885   }
0886   e.put(std::move(outputTrajs));
0887   e.put(std::move(outputTTAss));
0888   return;
0889 
0890 }  //end produce
0891 
0892 void TrackListMerger::returnEmptyCollections(edm::Event& e) {
0893   if (trkQualMod_) {
0894     auto vm = std::make_unique<edm::ValueMap<int>>();
0895     e.put(std::move(vm));
0896     auto quals = std::make_unique<QualityMaskCollection>();
0897     e.put(std::move(quals), "QualityMasks");
0898   } else {
0899     auto outputTrks = std::make_unique<reco::TrackCollection>();
0900     e.put(std::move(outputTrks));
0901 
0902     if (makeReKeyedSeeds_) {
0903       auto outputSeeds = std::make_unique<TrajectorySeedCollection>();
0904       e.put(std::move(outputSeeds));
0905     }
0906 
0907     if (copyExtras_) {
0908       auto outputTrkExtras = std::make_unique<reco::TrackExtraCollection>();
0909       auto outputTrkHits = std::make_unique<TrackingRecHitCollection>();
0910       e.put(std::move(outputTrkExtras));
0911       e.put(std::move(outputTrkHits));
0912     }
0913 
0914     auto outputTrajs = std::make_unique<std::vector<Trajectory>>();
0915     outputTTAss = std::make_unique<TrajTrackAssociationCollection>();
0916     e.put(std::move(outputTrajs));
0917     e.put(std::move(outputTTAss));
0918   }
0919   auto vmMVA = std::make_unique<edm::ValueMap<float>>();
0920   e.put(std::move(vmMVA), "MVAVals");
0921   auto mvas = std::make_unique<MVACollection>();
0922   e.put(std::move(mvas), "MVAValues");
0923   return;
0924 }
0925 
0926 #include "FWCore/PluginManager/interface/ModuleDef.h"
0927 #include "FWCore/Framework/interface/MakerMacros.h"
0928 
0929 DEFINE_FWK_MODULE(TrackListMerger);
[ 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!