File indexing completed on 2023-03-17 11:25:45
0001
0002
0003 #include "DataFormats/Common/interface/Ref.h"
0004
0005 #include "TrackAssociatorByHitsImpl.h"
0006
0007 #include "SimTracker/TrackerHitAssociation/interface/TrackerHitAssociator.h"
0008 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0009 #include "FWCore/Utilities/interface/Exception.h"
0010
0011
0012 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0013 #include "DataFormats/TrackReco/interface/Track.h"
0014 #include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
0015
0016 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
0017 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticleFwd.h"
0018 #include "SimGeneral/TrackingAnalysis/interface/SimHitTPAssociationProducer.h"
0019 #include "SimDataFormats/EncodedEventId/interface/EncodedEventId.h"
0020
0021 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0022 #include "DataFormats/DetId/interface/DetId.h"
0023 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0024 #include "DataFormats/SiStripDetId/interface/SiStripDetId.h"
0025 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0026 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0027
0028 using namespace reco;
0029 using namespace std;
0030
0031
0032
0033 TrackAssociatorByHitsImpl::TrackAssociatorByHitsImpl(edm::EDProductGetter const& productGetter,
0034 std::unique_ptr<TrackerHitAssociator> iAssociate,
0035 TrackerTopology const* iTopo,
0036 SimHitTPAssociationList const* iSimHitsTPAssoc,
0037 SimToRecoDenomType iSimToRecoDenominator,
0038 double iQuality_SimToReco,
0039 double iPurity_SimToReco,
0040 double iCut_RecoToSim,
0041 bool iUsePixels,
0042 bool iUseGrouped,
0043 bool iUseSplitting,
0044 bool iThreeHitTracksAreSpecial,
0045 bool iAbsoluteNumberOfHits)
0046 : productGetter_(&productGetter),
0047 associate(std::move(iAssociate)),
0048 tTopo(iTopo),
0049 simHitsTPAssoc(iSimHitsTPAssoc),
0050 SimToRecoDenominator(iSimToRecoDenominator),
0051 quality_SimToReco(iQuality_SimToReco),
0052 purity_SimToReco(iPurity_SimToReco),
0053 cut_RecoToSim(iCut_RecoToSim),
0054 UsePixels(iUsePixels),
0055 UseGrouped(iUseGrouped),
0056 UseSplitting(iUseSplitting),
0057 ThreeHitTracksAreSpecial(iThreeHitTracksAreSpecial),
0058 AbsoluteNumberOfHits(iAbsoluteNumberOfHits) {}
0059
0060
0061
0062
0063
0064
0065
0066
0067
0068
0069
0070
0071
0072
0073
0074
0075
0076
0077
0078
0079
0080
0081
0082
0083
0084
0085
0086
0087
0088
0089
0090
0091
0092 RecoToSimCollection TrackAssociatorByHitsImpl::associateRecoToSim(
0093 const edm::RefToBaseVector<reco::Track>& tC,
0094 const edm::RefVector<TrackingParticleCollection>& TPCollectionH) const {
0095
0096 int nshared = 0;
0097 float quality = 0;
0098 std::vector<SimHitIdpr> SimTrackIds;
0099 std::vector<SimHitIdpr> matchedIds;
0100 RecoToSimCollection outputCollection(productGetter_);
0101
0102
0103 std::vector<TrackingParticle const*> tPC;
0104 tPC.reserve(tPC.size());
0105 for (auto const& ref : TPCollectionH) {
0106 tPC.push_back(&(*ref));
0107 }
0108
0109
0110 int tindex = 0;
0111 for (edm::RefToBaseVector<reco::Track>::const_iterator track = tC.begin(); track != tC.end(); track++, tindex++) {
0112 matchedIds.clear();
0113 int ri = 0;
0114
0115 getMatchedIds<trackingRecHit_iterator>(
0116 matchedIds, SimTrackIds, ri, (*track)->recHitsBegin(), (*track)->recHitsEnd(), associate.get());
0117
0118
0119
0120
0121
0122
0123
0124
0125
0126 std::vector<SimHitIdpr> idcachev;
0127 if (!matchedIds.empty()) {
0128 int tpindex = 0;
0129 for (auto t = tPC.begin(); t != tPC.end(); ++t, ++tpindex) {
0130
0131
0132 idcachev.clear();
0133 nshared = getShared(matchedIds, idcachev, **t);
0134
0135
0136 if (abs((*t)->pdgId()) == 11 && ((*t)->g4Track_end() - (*t)->g4Track_begin()) > 1) {
0137 nshared -= getDoubleCount<trackingRecHit_iterator>(
0138 (*track)->recHitsBegin(), (*track)->recHitsEnd(), associate.get(), **t);
0139 }
0140
0141 if (AbsoluteNumberOfHits)
0142 quality = static_cast<double>(nshared);
0143 else if (ri != 0)
0144 quality = (static_cast<double>(nshared) / static_cast<double>(ri));
0145 else
0146 quality = 0;
0147
0148
0149 if (quality > cut_RecoToSim && !(ThreeHitTracksAreSpecial && ri == 3 && nshared < 3)) {
0150
0151 outputCollection.insert(tC[tindex], std::make_pair(TPCollectionH[tpindex], quality));
0152
0153
0154
0155
0156
0157 } else {
0158
0159
0160 }
0161 }
0162 }
0163 }
0164
0165 outputCollection.post_insert();
0166 return outputCollection;
0167 }
0168
0169 SimToRecoCollection TrackAssociatorByHitsImpl::associateSimToReco(
0170 const edm::RefToBaseVector<reco::Track>& tC,
0171 const edm::RefVector<TrackingParticleCollection>& TPCollectionH) const {
0172
0173
0174
0175
0176 float quality = 0;
0177 int nshared = 0;
0178 std::vector<SimHitIdpr> SimTrackIds;
0179 std::vector<SimHitIdpr> matchedIds;
0180 SimToRecoCollection outputCollection(productGetter_);
0181
0182
0183 std::vector<TrackingParticle const*> tPC;
0184 tPC.reserve(tPC.size());
0185 for (auto const& ref : TPCollectionH) {
0186 tPC.push_back(&(*ref));
0187 }
0188
0189
0190
0191
0192
0193
0194
0195
0196
0197
0198
0199
0200
0201 int tindex = 0;
0202 for (edm::RefToBaseVector<reco::Track>::const_iterator track = tC.begin(); track != tC.end(); track++, tindex++) {
0203
0204 int ri = 0;
0205 getMatchedIds<trackingRecHit_iterator>(
0206 matchedIds, SimTrackIds, ri, (*track)->recHitsBegin(), (*track)->recHitsEnd(), associate.get());
0207
0208
0209 std::vector<SimHitIdpr> idcachev;
0210 if (!matchedIds.empty()) {
0211 int tpindex = 0;
0212 for (auto t = tPC.begin(); t != tPC.end(); ++t, ++tpindex) {
0213 idcachev.clear();
0214 float totsimhit = 0;
0215
0216
0217 std::vector<PSimHit> tphits;
0218
0219
0220 nshared = getShared(matchedIds, idcachev, **t);
0221
0222
0223
0224
0225
0226
0227
0228
0229 if (nshared != 0) {
0230
0231
0232
0233
0234 std::pair<TrackingParticleRef, TrackPSimHitRef> clusterTPpairWithDummyTP(
0235 TPCollectionH[tpindex], TrackPSimHitRef());
0236
0237 auto range = std::equal_range(simHitsTPAssoc->begin(),
0238 simHitsTPAssoc->end(),
0239 clusterTPpairWithDummyTP,
0240 SimHitTPAssociationProducer::simHitTPAssociationListGreater);
0241 for (auto ip = range.first; ip != range.second; ++ip) {
0242 TrackPSimHitRef TPhit = ip->second;
0243 DetId dId = DetId(TPhit->detUnitId());
0244
0245 unsigned int subdetId = static_cast<unsigned int>(dId.subdetId());
0246 if (!UsePixels && (subdetId == PixelSubdetector::PixelBarrel || subdetId == PixelSubdetector::PixelEndcap))
0247 continue;
0248
0249
0250 SiStripDetId* stripDetId = nullptr;
0251 if (subdetId == SiStripDetId::TIB || subdetId == SiStripDetId::TOB || subdetId == SiStripDetId::TID ||
0252 subdetId == SiStripDetId::TEC)
0253 stripDetId = new SiStripDetId(dId);
0254
0255
0256
0257 bool newhit = true;
0258 for (std::vector<PSimHit>::const_iterator TPhitOK = tphits.begin(); TPhitOK != tphits.end(); TPhitOK++) {
0259 DetId dIdOK = DetId(TPhitOK->detUnitId());
0260
0261
0262
0263
0264
0265 if (!UseGrouped && !UseSplitting)
0266 if (tTopo->layer(dId) == tTopo->layer(dIdOK) && dId.subdetId() == dIdOK.subdetId())
0267 newhit = false;
0268
0269 if (!UseGrouped && UseSplitting)
0270 if (tTopo->layer(dId) == tTopo->layer(dIdOK) && dId.subdetId() == dIdOK.subdetId() &&
0271 (stripDetId == nullptr || stripDetId->partnerDetId() != dIdOK.rawId()))
0272 newhit = false;
0273
0274 if (UseGrouped && !UseSplitting)
0275 if (tTopo->layer(dId) == tTopo->layer(dIdOK) && dId.subdetId() == dIdOK.subdetId() &&
0276 stripDetId != nullptr && stripDetId->partnerDetId() == dIdOK.rawId())
0277 newhit = false;
0278
0279 if (UseGrouped && UseSplitting)
0280 newhit = true;
0281 }
0282 if (newhit) {
0283
0284 tphits.push_back(*TPhit);
0285 }
0286
0287 delete stripDetId;
0288 }
0289 totsimhit = tphits.size();
0290 }
0291
0292 if (AbsoluteNumberOfHits)
0293 quality = static_cast<double>(nshared);
0294 else if (SimToRecoDenominator == denomsim && totsimhit != 0)
0295 quality = ((double)nshared) / ((double)totsimhit);
0296 else if (SimToRecoDenominator == denomreco && ri != 0)
0297 quality = ((double)nshared) / ((double)ri);
0298 else
0299 quality = 0;
0300
0301
0302
0303 float purity = 1.0 * nshared / ri;
0304 if (quality > quality_SimToReco && !(ThreeHitTracksAreSpecial && totsimhit == 3 && nshared < 3) &&
0305 (AbsoluteNumberOfHits || (purity > purity_SimToReco))) {
0306
0307 outputCollection.insert(TPCollectionH[tpindex], std::make_pair(tC[tindex], quality));
0308
0309
0310
0311
0312 } else {
0313
0314
0315
0316 }
0317 }
0318 }
0319 }
0320
0321 outputCollection.post_insert();
0322 return outputCollection;
0323 }
0324
0325 RecoToSimCollectionSeed TrackAssociatorByHitsImpl::associateRecoToSim(
0326 const edm::Handle<edm::View<TrajectorySeed> >& seedCollectionH,
0327 const edm::Handle<TrackingParticleCollection>& TPCollectionH) const {
0328 edm::LogVerbatim("TrackAssociator") << "Starting TrackAssociatorByHitsImpl::associateRecoToSim - #seeds="
0329 << seedCollectionH->size() << " #TPs=" << TPCollectionH->size();
0330 int nshared = 0;
0331 float quality = 0;
0332 std::vector<SimHitIdpr> SimTrackIds;
0333 std::vector<SimHitIdpr> matchedIds;
0334 RecoToSimCollectionSeed outputCollection(productGetter_);
0335
0336 const TrackingParticleCollection& tPC = *(TPCollectionH.product());
0337
0338 const edm::View<TrajectorySeed> sC = *(seedCollectionH.product());
0339
0340
0341 int tindex = 0;
0342 for (edm::View<TrajectorySeed>::const_iterator seed = sC.begin(); seed != sC.end(); seed++, tindex++) {
0343 matchedIds.clear();
0344 int ri = 0;
0345 int nsimhit = seed->nHits();
0346 LogTrace("TrackAssociator") << "\nNEW SEED - seed number " << tindex << " # valid=" << nsimhit;
0347 getMatchedIds<edm::OwnVector<TrackingRecHit>::const_iterator>(
0348 matchedIds, SimTrackIds, ri, seed->recHits().begin(), seed->recHits().end(), associate.get());
0349
0350
0351 std::vector<SimHitIdpr> idcachev;
0352 if (!matchedIds.empty()) {
0353 int tpindex = 0;
0354 for (TrackingParticleCollection::const_iterator t = tPC.begin(); t != tPC.end(); ++t, ++tpindex) {
0355 LogTrace("TrackAssociator") << "TP number " << tpindex << " pdgId=" << t->pdgId()
0356 << " with number of PSimHits: " << nsimhit;
0357 idcachev.clear();
0358 nshared = getShared(matchedIds, idcachev, *t);
0359
0360
0361 if (abs(t->pdgId()) == 11 && (t->g4Track_end() - t->g4Track_begin()) > 1) {
0362 nshared -= getDoubleCount<edm::OwnVector<TrackingRecHit>::const_iterator>(
0363 seed->recHits().begin(), seed->recHits().end(), associate.get(), *t);
0364 }
0365
0366 if (AbsoluteNumberOfHits)
0367 quality = static_cast<double>(nshared);
0368 else if (ri != 0)
0369 quality = (static_cast<double>(nshared) / static_cast<double>(ri));
0370 else
0371 quality = 0;
0372
0373 if (quality > cut_RecoToSim && !(ThreeHitTracksAreSpecial && ri == 3 && nshared < 3)) {
0374
0375 outputCollection.insert(
0376 edm::RefToBase<TrajectorySeed>(seedCollectionH, tindex),
0377 std::make_pair(edm::Ref<TrackingParticleCollection>(TPCollectionH, tpindex), quality));
0378 LogTrace("TrackAssociator") << "Seed number " << tindex << " with #hits=" << ri
0379 << "associated to TP (pdgId, nb segments, p) = " << (*t).pdgId() << " "
0380 << (*t).g4Tracks().size() << " " << (*t).momentum() << " number " << tpindex
0381 << " with quality =" << quality;
0382 } else {
0383 LogTrace("TrackAssociator") << "Seed number " << tindex << " with #hits=" << ri
0384 << " NOT associated with quality =" << quality;
0385 }
0386 }
0387 }
0388 }
0389 LogTrace("TrackAssociator") << "% of Assoc Seeds="
0390 << ((double)outputCollection.size()) / ((double)seedCollectionH->size());
0391
0392 outputCollection.post_insert();
0393 return outputCollection;
0394 }
0395
0396 SimToRecoCollectionSeed TrackAssociatorByHitsImpl::associateSimToReco(
0397 const edm::Handle<edm::View<TrajectorySeed> >& seedCollectionH,
0398 const edm::Handle<TrackingParticleCollection>& TPCollectionH) const {
0399 edm::LogVerbatim("TrackAssociator") << "Starting TrackAssociatorByHitsImpl::associateSimToReco - #seeds="
0400 << seedCollectionH->size() << " #TPs=" << TPCollectionH->size();
0401 float quality = 0;
0402 int nshared = 0;
0403 std::vector<SimHitIdpr> SimTrackIds;
0404 std::vector<SimHitIdpr> matchedIds;
0405 SimToRecoCollectionSeed outputCollection(productGetter_);
0406
0407 const TrackingParticleCollection& tPC = *TPCollectionH.product();
0408
0409 const edm::View<TrajectorySeed> sC = *(seedCollectionH.product());
0410
0411
0412 int tindex = 0;
0413 for (edm::View<TrajectorySeed>::const_iterator seed = sC.begin(); seed != sC.end(); seed++, tindex++) {
0414 int ri = 0;
0415 LogTrace("TrackAssociator") << "\nNEW SEED - seed number " << tindex << " # valid=" << seed->nHits();
0416 getMatchedIds<edm::OwnVector<TrackingRecHit>::const_iterator>(
0417 matchedIds, SimTrackIds, ri, seed->recHits().begin(), seed->recHits().end(), associate.get());
0418
0419
0420 std::vector<SimHitIdpr> idcachev;
0421 if (!matchedIds.empty()) {
0422 int tpindex = 0;
0423 for (TrackingParticleCollection::const_iterator t = tPC.begin(); t != tPC.end(); ++t, ++tpindex) {
0424 idcachev.clear();
0425 int nsimhit = t->numberOfTrackerHits();
0426 LogTrace("TrackAssociator") << "TP number " << tpindex << " pdgId=" << t->pdgId()
0427 << " with number of PSimHits: " << nsimhit;
0428 nshared = getShared(matchedIds, idcachev, *t);
0429
0430 if (AbsoluteNumberOfHits)
0431 quality = static_cast<double>(nshared);
0432 else if (ri != 0)
0433 quality = ((double)nshared) / ((double)ri);
0434 else
0435 quality = 0;
0436
0437
0438
0439 if (quality > quality_SimToReco && !(ThreeHitTracksAreSpecial && ri == 3 && nshared < 3)) {
0440 outputCollection.insert(edm::Ref<TrackingParticleCollection>(TPCollectionH, tpindex),
0441 std::make_pair(edm::RefToBase<TrajectorySeed>(seedCollectionH, tindex), quality));
0442 LogTrace("TrackAssociator") << "TrackingParticle number " << tpindex << " with #hits=" << nsimhit
0443 << " associated to seed number " << tindex << " with #hits=" << ri
0444 << " with hit quality =" << quality;
0445 } else {
0446 LogTrace("TrackAssociator") << "TrackingParticle number " << tpindex << " with #hits=" << nsimhit
0447 << " NOT associated with quality =" << quality;
0448 }
0449 }
0450 }
0451 }
0452 LogTrace("TrackAssociator") << "% of Assoc TPs="
0453 << ((double)outputCollection.size()) / ((double)TPCollectionH->size());
0454 outputCollection.post_insert();
0455 return outputCollection;
0456 }
0457
0458 template <typename iter>
0459 void TrackAssociatorByHitsImpl::getMatchedIds(std::vector<SimHitIdpr>& matchedIds,
0460 std::vector<SimHitIdpr>& SimTrackIds,
0461 int& ri,
0462 iter begin,
0463 iter end,
0464 TrackerHitAssociator* associate) const {
0465 matchedIds.clear();
0466 ri = 0;
0467 for (iter it = begin; it != end; it++) {
0468 const TrackingRecHit* hit = getHitPtr(it);
0469 if (hit->isValid()) {
0470 ri++;
0471 uint32_t t_detID = hit->geographicalId().rawId();
0472 SimTrackIds.clear();
0473 associate->associateHitId(*hit, SimTrackIds);
0474
0475 if (!SimTrackIds.empty()) {
0476 for (size_t j = 0; j < SimTrackIds.size(); j++) {
0477 LogTrace("TrackAssociator") << " hit # " << ri << " valid=" << hit->isValid() << " det id = " << t_detID
0478 << " SimId " << SimTrackIds[j].first << " evt=" << SimTrackIds[j].second.event()
0479 << " bc=" << SimTrackIds[j].second.bunchCrossing();
0480 matchedIds.push_back(SimTrackIds[j]);
0481 }
0482 }
0483
0484
0485
0486
0487
0488
0489
0490
0491
0492
0493
0494
0495
0496
0497
0498
0499 } else {
0500 LogTrace("TrackAssociator") << "\t\t Invalid Hit On " << hit->geographicalId().rawId();
0501 }
0502 }
0503 }
0504
0505 int TrackAssociatorByHitsImpl::getShared(std::vector<SimHitIdpr>& matchedIds,
0506 std::vector<SimHitIdpr>& idcachev,
0507 TrackingParticle const& t) const {
0508 int nshared = 0;
0509 if (t.numberOfHits() == 0)
0510 return nshared;
0511
0512 for (size_t j = 0; j < matchedIds.size(); j++) {
0513
0514 if (find(idcachev.begin(), idcachev.end(), matchedIds[j]) == idcachev.end()) {
0515
0516 idcachev.push_back(matchedIds[j]);
0517
0518 for (TrackingParticle::g4t_iterator g4T = t.g4Track_begin(); g4T != t.g4Track_end(); ++g4T) {
0519
0520
0521
0522
0523
0524
0525
0526 if ((*g4T).trackId() == matchedIds[j].first && t.eventId() == matchedIds[j].second) {
0527 int countedhits = std::count(matchedIds.begin(), matchedIds.end(), matchedIds[j]);
0528 nshared += countedhits;
0529
0530
0531
0532 }
0533 }
0534 }
0535 }
0536 return nshared;
0537 }
0538
0539 template <typename iter>
0540 int TrackAssociatorByHitsImpl::getDoubleCount(iter begin,
0541 iter end,
0542 TrackerHitAssociator* associate,
0543 TrackingParticle const& t) const {
0544 int doublecount = 0;
0545 std::vector<SimHitIdpr> SimTrackIdsDC;
0546
0547 for (iter it = begin; it != end; it++) {
0548 int idcount = 0;
0549 SimTrackIdsDC.clear();
0550 associate->associateHitId(*getHitPtr(it), SimTrackIdsDC);
0551
0552 if (SimTrackIdsDC.size() > 1) {
0553
0554 for (TrackingParticle::g4t_iterator g4T = t.g4Track_begin(); g4T != t.g4Track_end(); ++g4T) {
0555 if (find(SimTrackIdsDC.begin(),
0556 SimTrackIdsDC.end(),
0557 SimHitIdpr((*g4T).trackId(), SimTrackIdsDC.begin()->second)) != SimTrackIdsDC.end()) {
0558 idcount++;
0559 }
0560 }
0561 }
0562 if (idcount > 1)
0563 doublecount += (idcount - 1);
0564 }
0565 return doublecount;
0566 }