Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2025-01-08 03:36:32

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