File indexing completed on 2025-01-08 03:36:32
0001
0002
0003
0004
0005
0006
0007
0008
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
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;
0164 unsigned long long timeOv;
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 }
0205
0206 namespace {
0207 edm::ProductID clusterProductB(const TrackingRecHit* hit) {
0208 return reinterpret_cast<const BaseTrackerRecHit*>(hit)->firstClusterRef().id();
0209 }
0210 }
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
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
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
0295
0296
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
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
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
0349 void TrackListMerger::produce(edm::Event& e, const edm::EventSetup& es) {
0350
0351
0352
0353
0354
0355
0356
0357 edm::ESHandle<TrackAlgoPriorityOrder> priorityH = es.getHandle(priorityToken);
0358 auto const& trackAlgoPriorityOrder = *priorityH;
0359
0360
0361
0362
0363
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
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
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
0399 this->returnEmptyCollections(e);
0400 return;
0401 }
0402
0403 statCount.begin(rSize);
0404
0405
0406
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;
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
0482 indexG[i] = ngood++;
0483
0484 }
0485 }
0486 }
0487
0488 statCount.pre(ngood);
0489
0490
0491 typedef std::pair<unsigned int, const TrackingRecHit*> IHit;
0492 std::vector<std::vector<IHit>> rh1(ngood);
0493
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);
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
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
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
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;
0554 float score1 = score[k1];
0555
0556
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
0564 if (notActive[collNum2])
0565 continue;
0566
0567 int k2 = indexG[j];
0568
0569 int newQualityMask = -9;
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);
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
0589 int noverlap = 0;
0590 int firstoverlap = 0;
0591
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
0600 unsigned int jh = 0;
0601 unsigned int ih = 0;
0602 while (ih != nh1 && jh != nh2) {
0603
0604
0605
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
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 }
0630
0631 }
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;
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];
0646 };
0647
0648 if (dupfound) {
0649 float score2 = score[k2];
0650 constexpr float almostSame =
0651 0.01f;
0652 if (score1 - score2 > almostSame) {
0653 seti(i, j);
0654 } else if (score2 - score1 > almostSame) {
0655 seti(j, i);
0656 } else {
0657
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
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 }
0677 statCount.overlap();
0678
0679
0680
0681
0682
0683
0684 }
0685 else {
0686 statCount.noOverlap();
0687 }
0688
0689 if (selected[i] == 0)
0690 break;
0691 }
0692 }
0693 }
0694
0695 auto vmMVA = std::make_unique<edm::ValueMap<float>>();
0696 edm::ValueMap<float>::Filler fillerMVA(*vmMVA);
0697
0698
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);
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
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
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
0801
0802
0803 if (copyExtras_) {
0804 edm::RefToBase<TrajectorySeed> origSeedRef = track->seedRef();
0805
0806 if (makeReKeyedSeeds_) {
0807 bool doRekeyOnThisSeed = false;
0808
0809 edm::InputTag clusterRemovalInfos("");
0810
0811 if (origSeedRef->nHits() != 0) {
0812 TrackingRecHit const& hit = *origSeedRef->recHits().begin();
0813 if (hit.isValid()) {
0814 edm::ProductID pID = clusterProductB(&hit);
0815
0816
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 }
0822 }
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
0836 else {
0837
0838 outputSeeds->push_back(TrajectorySeed(*origSeedRef));
0839 }
0840 edm::Ref<TrajectorySeedCollection> pureRef(refTrajSeeds, outputSeeds->size() - 1);
0841 origSeedRef = edm::RefToBase<TrajectorySeed>(pureRef);
0842 }
0843
0844
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
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 }
0874
0875
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
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 }
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);