TrackerHitAssociation/src/TrackerHitAssociator.cc

0001 // File: TrackerHitAssociator.cc
0002
0003 #include <memory>
0004 #include <string>
0005 #include <vector>
0006
0007 #include "SimTracker/TrackerHitAssociation/interface/TrackerHitAssociator.h"
0008
0009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0010
0011 //--- for Geometry:
0012 #include "DataFormats/Common/interface/Ref.h"
0013 #include "DataFormats/DetId/interface/DetId.h"
0014 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
0015 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0016
0017 #include "SimDataFormats/CrossingFrame/interface/CrossingFrame.h"
0018
0019 //for accumulate
0020 #include <numeric>
0021 #include <iostream>
0022
0023 using namespace std;
0024 using namespace edm;
0025
0026 //
0027 // Constructor for Config helper class, using default parameters
0028 //
0029 TrackerHitAssociator::Config::Config(edm::ConsumesCollector&& iC)
0030     : doPixel_(true), doStrip_(true), useOTph2_(false), doTrackAssoc_(false), assocHitbySimTrack_(false) {
0031   if (doStrip_) {
0032     if (useOTph2_)
0033       ph2OTrToken_ = iC.consumes<edm::DetSetVector<PixelDigiSimLink>>(edm::InputTag("simSiPixelDigis", "Tracker"));
0034     else
0035       stripToken_ = iC.consumes<edm::DetSetVector<StripDigiSimLink>>(edm::InputTag("simSiStripDigis"));
0036   }
0037   if (doPixel_) {
0038     if (useOTph2_)
0039       pixelToken_ = iC.consumes<edm::DetSetVector<PixelDigiSimLink>>(edm::InputTag("simSiPixelDigis", "Pixel"));
0040     else
0041       pixelToken_ = iC.consumes<edm::DetSetVector<PixelDigiSimLink>>(edm::InputTag("simSiPixelDigis"));
0042   }
0043   if (!doTrackAssoc_) {
0044     std::vector<std::string> trackerContainers;
0045     trackerContainers.reserve(12);
0046     trackerContainers.emplace_back("g4SimHitsTrackerHitsTIBLowTof");
0047     trackerContainers.emplace_back("g4SimHitsTrackerHitsTIBHighTof");
0048     trackerContainers.emplace_back("g4SimHitsTrackerHitsTIDLowTof");
0049     trackerContainers.emplace_back("g4SimHitsTrackerHitsTIDHighTof");
0050     trackerContainers.emplace_back("g4SimHitsTrackerHitsTOBLowTof");
0051     trackerContainers.emplace_back("g4SimHitsTrackerHitsTOBHighTof");
0052     trackerContainers.emplace_back("g4SimHitsTrackerHitsTECLowTof");
0053     trackerContainers.emplace_back("g4SimHitsTrackerHitsTECHighTof");
0054     trackerContainers.emplace_back("g4SimHitsTrackerHitsPixelBarrelLowTof");
0055     trackerContainers.emplace_back("g4SimHitsTrackerHitsPixelBarrelHighTof");
0056     trackerContainers.emplace_back("g4SimHitsTrackerHitsPixelEndcapLowTof");
0057     trackerContainers.emplace_back("g4SimHitsTrackerHitsPixelEndcapHighTof");
0058     cfTokens_.reserve(trackerContainers.size());
0059     simHitTokens_.reserve(trackerContainers.size());
0060     for (auto const& trackerContainer : trackerContainers) {
0061       cfTokens_.push_back(iC.consumes<CrossingFrame<PSimHit>>(edm::InputTag("mix", trackerContainer)));
0062       simHitTokens_.push_back(iC.consumes<std::vector<PSimHit>>(edm::InputTag("g4SimHits", trackerContainer)));
0063     }
0064   }
0065 }
0066
0067 //
0068 // Constructor for Config helper class, using configured parameters
0069 //
0070 TrackerHitAssociator::Config::Config(const edm::ParameterSet& conf, edm::ConsumesCollector&& iC)
0071     : doPixel_(conf.getParameter<bool>("associatePixel")),
0072       doStrip_(conf.getParameter<bool>("associateStrip")),
0073       useOTph2_(conf.existsAs<bool>("usePhase2Tracker") ? conf.getParameter<bool>("usePhase2Tracker") : false),
0074       //
0075       doTrackAssoc_(conf.getParameter<bool>("associateRecoTracks")),
0076       assocHitbySimTrack_(
0077           conf.existsAs<bool>("associateHitbySimTrack") ? conf.getParameter<bool>("associateHitbySimTrack") : false) {
0078   if (doStrip_) {
0079     if (useOTph2_)
0080       ph2OTrToken_ =
0081           iC.consumes<edm::DetSetVector<PixelDigiSimLink>>(conf.getParameter<edm::InputTag>("phase2TrackerSimLinkSrc"));
0082     else
0083       stripToken_ =
0084           iC.consumes<edm::DetSetVector<StripDigiSimLink>>(conf.getParameter<edm::InputTag>("stripSimLinkSrc"));
0085   }
0086   if (doPixel_)
0087     pixelToken_ = iC.consumes<edm::DetSetVector<PixelDigiSimLink>>(conf.getParameter<edm::InputTag>("pixelSimLinkSrc"));
0088   if (!doTrackAssoc_) {
0089     std::vector<std::string> trackerContainers(conf.getParameter<std::vector<std::string>>("ROUList"));
0090     cfTokens_.reserve(trackerContainers.size());
0091     simHitTokens_.reserve(trackerContainers.size());
0092     for (auto const& trackerContainer : trackerContainers) {
0093       cfTokens_.push_back(iC.consumes<CrossingFrame<PSimHit>>(edm::InputTag("mix", trackerContainer)));
0094       simHitTokens_.push_back(iC.consumes<std::vector<PSimHit>>(edm::InputTag("g4SimHits", trackerContainer)));
0095     }
0096   }
0097 }
0098
0099 void TrackerHitAssociator::fillPSetDescription(edm::ParameterSetDescription& desc) {
0100   desc.setComment("auxilliary class to store information about recHit/simHit association");
0101   desc.add<bool>("associatePixel", false);
0102   desc.add<bool>("associateStrip", false);
0103   desc.add<bool>("usePhase2Tracker", false);
0104   desc.add<bool>("associateRecoTracks", false);
0105   desc.add<bool>("associateHitbySimTrack", false);
0106   desc.add<edm::InputTag>("phase2TrackerSimLinkSrc", edm::InputTag("simSiPixelDigis", "Tracker"));
0107   desc.add<edm::InputTag>("stripSimLinkSrc", edm::InputTag("simSiStripDigis"));
0108   desc.add<edm::InputTag>("pixelSimLinkSrc", edm::InputTag("simSiPixelDigis"));
0109   desc.add<std::vector<std::string>>(
0110       "ROUList", {"TrackerHitsTIBLowTof", "TrackerHitsTIBHighTof", "TrackerHitsTOBLowTof", "TrackerHitsTOBHighTof"});
0111 }
0112
0113 //
0114 // Constructor supporting consumes interface
0115 //
0116 TrackerHitAssociator::TrackerHitAssociator(const edm::Event& e, const TrackerHitAssociator::Config& config)
0117     : doPixel_(config.doPixel_),
0118       doStrip_(config.doStrip_),
0119       useOTph2_(config.useOTph2_),
0120       doTrackAssoc_(config.doTrackAssoc_),
0121       assocHitbySimTrack_(config.assocHitbySimTrack_) {
0122   //if track association there is no need to access the input collections
0123   if (!doTrackAssoc_) {
0124     makeMaps(e, config);
0125   }
0126
0127   if (doStrip_) {
0128     if (useOTph2_)
0129       e.getByToken(config.ph2OTrToken_, ph2trackerdigisimlink);
0130     else
0131       e.getByToken(config.stripToken_, stripdigisimlink);
0132   }
0133   if (doPixel_)
0134     e.getByToken(config.pixelToken_, pixeldigisimlink);
0135 }
0136
0137 void TrackerHitAssociator::makeMaps(const edm::Event& theEvent, const TrackerHitAssociator::Config& config) {
0138   // Step A: Get Inputs
0139   //  The collections are specified via ROUList in the configuration, and can
0140   //  be either crossing frames (e.g., mix/g4SimHitsTrackerHitsTIBLowTof)
0141   //  or just PSimHits (e.g., g4SimHits/TrackerHitsTIBLowTof)
0142   if (assocHitbySimTrack_) {
0143     for (auto const& cfToken : config.cfTokens_) {
0144       edm::Handle<CrossingFrame<PSimHit>> cf_simhit;
0145       int Nhits = 0;
0146       if (theEvent.getByToken(cfToken, cf_simhit)) {
0147         std::unique_ptr<MixCollection<PSimHit>> thisContainerHits(new MixCollection<PSimHit>(cf_simhit.product()));
0148         for (auto const& isim : *thisContainerHits) {
0149           DetId theDet(isim.detUnitId());
0150           SimHitMap[theDet].push_back(isim);
0151           ++Nhits;
0152         }
0153         LogDebug("TrkHitAssocTrace") << "simHits from crossing frames; map size = " << SimHitMap.size()
0154                                      << ", Hit count = " << Nhits << std::endl;
0155       }
0156     }
0157     for (auto const& simHitToken : config.simHitTokens_) {
0158       edm::Handle<std::vector<PSimHit>> simHits;
0159       int Nhits = 0;
0160       if (theEvent.getByToken(simHitToken, simHits)) {
0161         for (auto const& isim : *simHits) {
0162           DetId theDet(isim.detUnitId());
0163           SimHitMap[theDet].push_back(isim);
0164           ++Nhits;
0165         }
0166         LogDebug("TrkHitAssocTrace") << "simHits from prompt collections; map size = " << SimHitMap.size()
0167                                      << ", Hit count = " << Nhits << std::endl;
0168       }
0169     }
0170   } else {  // !assocHitbySimTrack_
0171     const char* const highTag = "HighTof";
0172     unsigned int tofBin;
0173     edm::EDConsumerBase::Labels labels;
0174     subDetTofBin theSubDetTofBin;
0175     unsigned int collectionIndex = 0;
0176     for (auto const& cfToken : config.cfTokens_) {
0177       collectionIndex++;
0178       edm::Handle<CrossingFrame<PSimHit>> cf_simhit;
0179       int Nhits = 0;
0180       if (theEvent.getByToken(cfToken, cf_simhit)) {
0181         std::unique_ptr<MixCollection<PSimHit>> thisContainerHits(new MixCollection<PSimHit>(cf_simhit.product()));
0182         theEvent.labelsForToken(cfToken, labels);
0183         if (std::strstr(labels.productInstance, highTag) != nullptr) {
0184           tofBin = StripDigiSimLink::HighTof;
0185         } else {
0186           tofBin = StripDigiSimLink::LowTof;
0187         }
0188         for (auto const& isim : *thisContainerHits) {
0189           DetId theDet(isim.detUnitId());
0190           theSubDetTofBin = std::make_pair(theDet.subdetId(), tofBin);
0191           SimHitCollMap[theSubDetTofBin] = collectionIndex;
0192           SimHitMap[SimHitCollMap[theSubDetTofBin]].push_back(isim);
0193           ++Nhits;
0194         }
0195         LogDebug("TrkHitAssocTrace") << "simHits from crossing frames " << collectionIndex << ":  " << Nhits
0196                                      << std::endl;
0197       }
0198     }
0199     collectionIndex = 0;
0200     for (auto const& simHitToken : config.simHitTokens_) {
0201       collectionIndex++;
0202       edm::Handle<std::vector<PSimHit>> simHits;
0203       int Nhits = 0;
0204       if (theEvent.getByToken(simHitToken, simHits)) {
0205         theEvent.labelsForToken(simHitToken, labels);
0206         if (std::strstr(labels.productInstance, highTag) != nullptr) {
0207           tofBin = StripDigiSimLink::HighTof;
0208         } else {
0209           tofBin = StripDigiSimLink::LowTof;
0210         }
0211         for (auto const& isim : *simHits) {
0212           DetId theDet(isim.detUnitId());
0213           theSubDetTofBin = std::make_pair(theDet.subdetId(), tofBin);
0214           SimHitCollMap[theSubDetTofBin] = collectionIndex;
0215           SimHitMap[SimHitCollMap[theSubDetTofBin]].push_back(isim);
0216           ++Nhits;
0217         }
0218         LogDebug("TrkHitAssocTrace") << "simHits from prompt collection " << collectionIndex << ":  " << Nhits
0219                                      << std::endl;
0220       }
0221     }
0222   }
0223 }
0224
0225 std::vector<PSimHit> TrackerHitAssociator::associateHit(const TrackingRecHit& thit) const {
0226   if (const SiTrackerMultiRecHit* rechit = dynamic_cast<const SiTrackerMultiRecHit*>(&thit)) {
0227     return associateMultiRecHit(rechit);
0228   }
0229
0230   //vector with the matched SimHit
0231   std::vector<PSimHit> result;
0232
0233   if (doTrackAssoc_)
0234     return result;  // We don't want the SimHits for this RecHit
0235
0236   // Vectors to contain lists of matched simTracks, simHits
0237   std::vector<SimHitIdpr> simtrackid;
0238   std::vector<simhitAddr> simhitCFPos;
0239
0240   //get the Detector type of the rechit
0241   DetId detid = thit.geographicalId();
0242   uint32_t detID = detid.rawId();
0243
0244   // Get the vectors of simtrackIDs and simHit addresses associated with this rechit
0245   associateHitId(thit, simtrackid, &simhitCFPos);
0246   LogDebug("TrkHitAssocTrace") << printDetBnchEvtTrk(detid, detID, simtrackid);
0247
0248   // Get the vector of simHits associated with this rechit
0249   if (!assocHitbySimTrack_ && !simhitCFPos.empty()) {
0250     // We use the indices to the simHit collections taken
0251     //  from the DigiSimLinks and returned in simhitCFPos.
0252     //  simhitCFPos[i] contains the full address of the ith simhit:
0253     //   <collection index, simhit index>
0254
0255     //check if the recHit is a SiStripMatchedRecHit2D
0256     if (dynamic_cast<const SiStripMatchedRecHit2D*>(&thit)) {
0257       for (auto const& theSimHitAddr : simhitCFPos) {
0258         simHitCollectionID theSimHitCollID = theSimHitAddr.first;
0259         auto it = SimHitMap.find(theSimHitCollID);
0260         if (it != SimHitMap.end()) {
0261           unsigned int theSimHitIndex = theSimHitAddr.second;
0262           if (theSimHitIndex < (it->second).size()) {
0263             const PSimHit& theSimHit = (it->second)[theSimHitIndex];
0264             // Try to remove ghosts by requiring a match to the simTrack also
0265             unsigned int simHitid = theSimHit.trackId();
0266             EncodedEventId simHiteid = theSimHit.eventId();
0267             for (auto const& id : simtrackid) {
0268               if (simHitid == id.first && simHiteid == id.second) {
0269                 result.push_back(theSimHit);
0270               }
0271             }
0272             LogDebug("TrkHitAssocTrace") << "by CFpos, simHit detId =  " << theSimHit.detUnitId() << " address = ("
0273                                          << theSimHitAddr.first << ", " << theSimHitIndex
0274                                          << "), process = " << theSimHit.processType() << " ("
0275                                          << theSimHit.eventId().bunchCrossing() << ", " << theSimHit.eventId().event()
0276                                          << ", " << theSimHit.trackId() << ")" << std::endl;
0277           }
0278         }
0279       }
0280     } else {  // Not a SiStripMatchedRecHit2D
0281       for (auto const& theSimHitAddr : simhitCFPos) {
0282         simHitCollectionID theSimHitCollID = theSimHitAddr.first;
0283         auto it = SimHitMap.find(theSimHitCollID);
0284         if (it != SimHitMap.end()) {
0285           unsigned int theSimHitIndex = theSimHitAddr.second;
0286           if (theSimHitIndex < (it->second).size()) {
0287             result.push_back((it->second)[theSimHitIndex]);
0288             LogDebug("TrkHitAssocTrace") << "by CFpos, simHit detId =  " << (it->second)[theSimHitIndex].detUnitId()
0289                                          << " address = (" << theSimHitCollID << ", " << theSimHitIndex
0290                                          << "), process = " << (it->second)[theSimHitIndex].processType() << " ("
0291                                          << (it->second)[theSimHitIndex].eventId().bunchCrossing() << ", "
0292                                          << (it->second)[theSimHitIndex].eventId().event() << ", "
0293                                          << (it->second)[theSimHitIndex].trackId() << ")" << std::endl;
0294           }
0295         }
0296       }
0297     }
0298     return result;
0299   }  // if !assocHitbySimTrack
0300
0301   // Get the SimHit from the trackid instead
0302   auto it = SimHitMap.find(detID);
0303   if (it != SimHitMap.end()) {
0304     for (auto const& ihit : it->second) {
0305       unsigned int simHitid = ihit.trackId();
0306       EncodedEventId simHiteid = ihit.eventId();
0307       for (auto id : simtrackid) {
0308         if (simHitid == id.first && simHiteid == id.second) {
0309           result.push_back(ihit);
0310           LogDebug("TrkHitAssocTrace") << "by TrackID, simHit detId =  " << ihit.detUnitId()
0311                                        << ", process = " << ihit.processType() << " (" << ihit.eventId().bunchCrossing()
0312                                        << ", " << ihit.eventId().event() << ", " << ihit.trackId() << ")" << std::endl;
0313           break;
0314         }
0315       }
0316     }
0317
0318   } else {
0319     /// Check if it's the gluedDet.
0320     auto itrphi = SimHitMap.find(detID + 2);  //iterator to the simhit in the rphi module
0321     auto itster = SimHitMap.find(detID + 1);  //iterator to the simhit in the stereo module
0322     if (itrphi != SimHitMap.end() && itster != SimHitMap.end()) {
0323       std::vector<PSimHit> simHitVector = itrphi->second;
0324       simHitVector.insert(simHitVector.end(), (itster->second).begin(), (itster->second).end());
0325       for (auto const& ihit : simHitVector) {
0326         unsigned int simHitid = ihit.trackId();
0327         EncodedEventId simHiteid = ihit.eventId();
0328         for (auto const& id : simtrackid) {
0329           if (simHitid == id.first && simHiteid == id.second) {
0330             result.push_back(ihit);
0331             LogDebug("TrkHitAssocTrace") << "by TrackID, simHit detId =  " << ihit.detUnitId()
0332                                          << ", process = " << ihit.processType() << " ("
0333                                          << ihit.eventId().bunchCrossing() << ", " << ihit.eventId().event() << ", "
0334                                          << ihit.trackId() << ")" << std::endl;
0335             break;
0336           }
0337         }
0338       }
0339     }
0340   }
0341
0342   return result;
0343 }
0344
0345 std::vector<SimHitIdpr> TrackerHitAssociator::associateHitId(const TrackingRecHit& thit) const {
0346   std::vector<SimHitIdpr> simhitid;
0347   associateHitId(thit, simhitid);
0348   return simhitid;
0349 }
0350
0351 void TrackerHitAssociator::associateHitId(const TrackingRecHit& thit,
0352                                           std::vector<SimHitIdpr>& simtkid,
0353                                           std::vector<simhitAddr>* simhitCFPos) const {
0354   simtkid.clear();
0355
0356   if (const SiTrackerMultiRecHit* rechit = dynamic_cast<const SiTrackerMultiRecHit*>(&thit))
0357     simtkid = associateMultiRecHitId(rechit, simhitCFPos);
0358
0359   //check if it is a simple SiStripRecHit2D
0360   if (const SiStripRecHit2D* rechit = dynamic_cast<const SiStripRecHit2D*>(&thit))
0361     associateSiStripRecHit(rechit, simtkid, simhitCFPos);
0362
0363   //check if it is a SiStripRecHit1D
0364   else if (const SiStripRecHit1D* rechit = dynamic_cast<const SiStripRecHit1D*>(&thit))
0365     associateSiStripRecHit(rechit, simtkid, simhitCFPos);
0366
0367   //check if it is a SiStripMatchedRecHit2D
0368   else if (const SiStripMatchedRecHit2D* rechit = dynamic_cast<const SiStripMatchedRecHit2D*>(&thit))
0369     simtkid = associateMatchedRecHit(rechit, simhitCFPos);
0370
0371   //check if it is a  ProjectedSiStripRecHit2D
0372   else if (const ProjectedSiStripRecHit2D* rechit = dynamic_cast<const ProjectedSiStripRecHit2D*>(&thit))
0373     simtkid = associateProjectedRecHit(rechit, simhitCFPos);
0374
0375   //check if it is a Phase2TrackerRecHit1D
0376   else if (const Phase2TrackerRecHit1D* rechit = dynamic_cast<const Phase2TrackerRecHit1D*>(&thit))
0377     associatePhase2TrackerRecHit(rechit, simtkid, simhitCFPos);
0378
0379   //check if it is a SiPixelRecHit
0380   else if (const SiPixelRecHit* rechit = dynamic_cast<const SiPixelRecHit*>(&thit))
0381     associatePixelRecHit(rechit, simtkid, simhitCFPos);
0382
0383   //check if these are GSRecHits (from FastSim)
0384   if (trackerHitRTTI::isFast(thit))
0385     simtkid = associateFastRecHit(static_cast<const FastTrackerRecHit*>(&thit));
0386 }
0387
0388 template <typename T>
0389 inline void TrackerHitAssociator::associateSiStripRecHit(const T* simplerechit,
0390                                                          std::vector<SimHitIdpr>& simtrackid,
0391                                                          std::vector<simhitAddr>* simhitCFPos) const {
0392   const SiStripCluster* clust = &(*simplerechit->cluster());
0393   associateSimpleRecHitCluster(clust, simplerechit->geographicalId(), simtrackid, simhitCFPos);
0394 }
0395
0396 //
0397 //  Method for obtaining simTracks and simHits from a cluster
0398 //
0399 void TrackerHitAssociator::associateCluster(const SiStripCluster* clust,
0400                                             const DetId& detid,
0401                                             std::vector<SimHitIdpr>& simtrackid,
0402                                             std::vector<PSimHit>& simhit) const {
0403   std::vector<simhitAddr> simhitCFPos;
0404   associateSimpleRecHitCluster(clust, detid, simtrackid, &simhitCFPos);
0405
0406   for (auto const& theSimHitAddr : simhitCFPos) {
0407     simHitCollectionID theSimHitCollID = theSimHitAddr.first;
0408     auto it = SimHitMap.find(theSimHitCollID);
0409
0410     if (it != SimHitMap.end()) {
0411       unsigned int theSimHitIndex = theSimHitAddr.second;
0412       if (theSimHitIndex < (it->second).size())
0413         simhit.push_back((it->second)[theSimHitIndex]);
0414       LogDebug("TrkHitAssocTrace") << "For cluster, simHit detId =  " << (it->second)[theSimHitIndex].detUnitId()
0415                                    << " address = (" << theSimHitCollID << ", " << theSimHitIndex
0416                                    << "), process = " << (it->second)[theSimHitIndex].processType()
0417                                    << " (bnch, evt, trk) = (" << (it->second)[theSimHitIndex].eventId().bunchCrossing()
0418                                    << ", " << (it->second)[theSimHitIndex].eventId().event() << ", "
0419                                    << (it->second)[theSimHitIndex].trackId() << ")" << std::endl;
0420     }
0421   }
0422 }
0423
0424 void TrackerHitAssociator::associateSimpleRecHitCluster(const SiStripCluster* clust,
0425                                                         const DetId& detid,
0426                                                         std::vector<SimHitIdpr>& simtrackid,
0427                                                         std::vector<simhitAddr>* simhitCFPos) const {
0428   uint32_t detID = detid.rawId();
0429   auto isearch = stripdigisimlink->find(detID);
0430   if (isearch != stripdigisimlink->end()) {  //if it is not empty
0431     auto link_detset = (*isearch);
0432
0433     if (clust != nullptr) {  //the cluster is valid
0434       int clusiz = clust->amplitudes().size();
0435       int first = clust->firstStrip();
0436       int last = first + clusiz;
0437
0438       LogDebug("TrkHitAssocDbg") << "Cluster size " << clusiz << " first strip = " << first
0439                                  << " last strip = " << last - 1 << std::endl
0440                                  << " detID = " << detID << " DETSET size = " << link_detset.data.size() << std::endl;
0441       int channel;
0442       for (const auto& linkiter : link_detset.data) {
0443         channel = (int)(linkiter.channel());
0444         if (channel >= first && channel < last) {
0445           LogDebug("TrkHitAssocDbg") << "Channel = " << std::setw(4) << linkiter.channel()
0446                                      << ", TrackID = " << std::setw(8) << linkiter.SimTrackId()
0447                                      << ", tofBin = " << std::setw(3) << linkiter.TofBin()
0448                                      << ", fraction = " << std::setw(8) << linkiter.fraction()
0449                                      << ", Position = " << linkiter.CFposition() << std::endl;
0450           SimHitIdpr currentId(linkiter.SimTrackId(), linkiter.eventId());
0451
0452           //create a vector with the list of SimTrack ID's of the tracks that contributed to the RecHit
0453           //write the id only once in the vector
0454
0455           if (find(simtrackid.begin(), simtrackid.end(), currentId) == simtrackid.end()) {
0456             LogDebug("TrkHitAssocDbg") << " Adding track id  = " << currentId.first
0457                                        << " Event id = " << currentId.second.event()
0458                                        << " Bunch Xing = " << currentId.second.bunchCrossing() << std::endl;
0459             simtrackid.push_back(currentId);
0460           }
0461
0462           if (simhitCFPos != nullptr) {
0463             //create a vector that contains all the positions (in the MixCollection) of the SimHits that contributed to the RecHit
0464             //write position only once
0465             unsigned int currentCFPos = linkiter.CFposition();
0466             unsigned int tofBin = linkiter.TofBin();
0467             subDetTofBin theSubDetTofBin = std::make_pair(detid.subdetId(), tofBin);
0468             auto it = SimHitCollMap.find(theSubDetTofBin);
0469             if (it != SimHitCollMap.end()) {
0470               simhitAddr currentAddr = std::make_pair(it->second, currentCFPos);
0471               if (find(simhitCFPos->begin(), simhitCFPos->end(), currentAddr) == simhitCFPos->end()) {
0472                 simhitCFPos->push_back(currentAddr);
0473               }
0474             }
0475           }
0476         }
0477       }
0478     } else {
0479       edm::LogError("TrackerHitAssociator") << "no cluster reference attached";
0480     }
0481   }
0482 }
0483
0484 std::vector<SimHitIdpr> TrackerHitAssociator::associateMatchedRecHit(const SiStripMatchedRecHit2D* matchedrechit,
0485                                                                      std::vector<simhitAddr>* simhitCFPos) const {
0486   std::vector<SimHitIdpr> matched_mono;
0487   std::vector<SimHitIdpr> matched_st;
0488
0489   const SiStripRecHit2D mono = matchedrechit->monoHit();
0490   const SiStripRecHit2D st = matchedrechit->stereoHit();
0491   //associate the two simple hits separately
0492   associateSiStripRecHit(&mono, matched_mono, simhitCFPos);
0493   associateSiStripRecHit(&st, matched_st, simhitCFPos);
0494
0495   //save in a vector all the simtrack-id's that are common to mono and stereo hits
0496   std::vector<SimHitIdpr> simtrackid;
0497   if (!(matched_mono.empty() || matched_st.empty())) {
0498     for (auto const& mhit : matched_mono) {
0499       //save only once the ID
0500       if (find(simtrackid.begin(), simtrackid.end(), mhit) == simtrackid.end()) {
0501         //save if the stereoID matched the monoID
0502         if (find(matched_st.begin(), matched_st.end(), mhit) != matched_st.end()) {
0503           simtrackid.push_back(mhit);
0504         }
0505       }
0506     }
0507   }
0508   return simtrackid;
0509 }
0510
0511 std::vector<SimHitIdpr> TrackerHitAssociator::associateProjectedRecHit(const ProjectedSiStripRecHit2D* projectedrechit,
0512                                                                        std::vector<simhitAddr>* simhitCFPos) const {
0513   //projectedRecHit is a "matched" rechit with only one component
0514
0515   std::vector<SimHitIdpr> matched_mono;
0516
0517   const SiStripRecHit2D mono = projectedrechit->originalHit();
0518   associateSiStripRecHit(&mono, matched_mono, simhitCFPos);
0519   return matched_mono;
0520 }
0521
0522 void TrackerHitAssociator::associatePhase2TrackerRecHit(const Phase2TrackerRecHit1D* rechit,
0523                                                         std::vector<SimHitIdpr>& simtrackid,
0524                                                         std::vector<simhitAddr>* simhitCFPos) const {
0525   //
0526   // Phase 2 outer tracker associator
0527   //
0528   DetId detid = rechit->geographicalId();
0529   uint32_t detID = detid.rawId();
0530
0531   auto isearch = ph2trackerdigisimlink->find(detID);
0532   if (isearch != ph2trackerdigisimlink->end()) {  //if it is not empty
0533     auto link_detset = (*isearch);
0534     Phase2TrackerRecHit1D::ClusterRef const& cluster = rechit->cluster();
0535
0536     //check the reference is valid
0537     if (!(cluster.isNull())) {  //if the cluster is valid
0538       int minRow = (*cluster).firstStrip();
0539       int maxRow = (*cluster).firstStrip() + (*cluster).size();
0540       int Col = (*cluster).column();
0541       LogDebug("TrkHitAssocDbg") << "    Cluster minRow " << minRow << " maxRow " << maxRow << " column " << Col
0542                                  << std::endl;
0543       int dsl = 0;
0544       for (auto const& linkiter : link_detset.data) {
0545         ++dsl;
0546         std::pair<int, int> coord = Phase2TrackerDigi::channelToPixel(linkiter.channel());
0547         LogDebug("TrkHitAssocDbg") << "    " << dsl << ") Digi link: row " << coord.first << " col " << coord.second
0548                                    << std::endl;
0549         if (coord.first <= maxRow && coord.first >= minRow && coord.second == Col) {
0550           LogDebug("TrkHitAssocDbg") << "      !-> trackid   " << linkiter.SimTrackId() << endl
0551                                      << "          fraction  " << linkiter.fraction() << endl;
0552           SimHitIdpr currentId(linkiter.SimTrackId(), linkiter.eventId());
0553           if (find(simtrackid.begin(), simtrackid.end(), currentId) == simtrackid.end()) {
0554             simtrackid.push_back(currentId);
0555           }
0556
0557           if (simhitCFPos != nullptr) {
0558             //create a vector that contains all the positions (in the MixCollection) of the SimHits that contributed to the RecHit
0559             //write position only once
0560             unsigned int currentCFPos = linkiter.CFposition();
0561             unsigned int tofBin = linkiter.TofBin();
0562             subDetTofBin theSubDetTofBin = std::make_pair(detid.subdetId(), tofBin);
0563             auto it = SimHitCollMap.find(theSubDetTofBin);
0564             if (it != SimHitCollMap.end()) {
0565               simhitAddr currentAddr = std::make_pair(it->second, currentCFPos);
0566               if (find(simhitCFPos->begin(), simhitCFPos->end(), currentAddr) == simhitCFPos->end()) {
0567                 simhitCFPos->push_back(currentAddr);
0568               }
0569             }
0570           }
0571         }
0572       }  // end of simlink loop
0573     } else {
0574       edm::LogError("TrackerHitAssociator") << "no Phase2 outer tracker cluster reference attached";
0575     }
0576   }
0577 }
0578
0579 void TrackerHitAssociator::associatePixelRecHit(const SiPixelRecHit* pixelrechit,
0580                                                 std::vector<SimHitIdpr>& simtrackid,
0581                                                 std::vector<simhitAddr>* simhitCFPos) const {
0582   //
0583   // Pixel associator
0584   //
0585   DetId detid = pixelrechit->geographicalId();
0586   uint32_t detID = detid.rawId();
0587
0588   auto isearch = pixeldigisimlink->find(detID);
0589   if (isearch != pixeldigisimlink->end()) {  //if it is not empty
0590     auto link_detset = (*isearch);
0591     SiPixelRecHit::ClusterRef const& cluster = pixelrechit->cluster();
0592
0593     //check the reference is valid
0594
0595     if (!(cluster.isNull())) {  //if the cluster is valid
0596
0597       int minPixelRow = (*cluster).minPixelRow();
0598       int maxPixelRow = (*cluster).maxPixelRow();
0599       int minPixelCol = (*cluster).minPixelCol();
0600       int maxPixelCol = (*cluster).maxPixelCol();
0601       LogDebug("TrkHitAssocDbg") << "    Cluster minRow " << minPixelRow << " maxRow " << maxPixelRow << std::endl
0602                                  << "    Cluster minCol " << minPixelCol << " maxCol " << maxPixelCol << std::endl;
0603       int dsl = 0;
0604       for (auto const& linkiter : link_detset.data) {
0605         ++dsl;
0606         std::pair<int, int> pixel_coord = PixelDigi::channelToPixel(linkiter.channel());
0607         LogDebug("TrkHitAssocDbg") << "    " << dsl << ") Digi link: row " << pixel_coord.first << " col "
0608                                    << pixel_coord.second << std::endl;
0609         if (pixel_coord.first <= maxPixelRow && pixel_coord.first >= minPixelRow && pixel_coord.second <= maxPixelCol &&
0610             pixel_coord.second >= minPixelCol) {
0611           LogDebug("TrkHitAssocDbg") << "      !-> trackid   " << linkiter.SimTrackId() << endl
0612                                      << "          fraction  " << linkiter.fraction() << endl;
0613           SimHitIdpr currentId(linkiter.SimTrackId(), linkiter.eventId());
0614           if (find(simtrackid.begin(), simtrackid.end(), currentId) == simtrackid.end()) {
0615             simtrackid.push_back(currentId);
0616           }
0617
0618           if (simhitCFPos != nullptr) {
0619             //create a vector that contains all the positions (in the MixCollection) of the SimHits that contributed to the RecHit
0620             //write position only once
0621             unsigned int currentCFPos = linkiter.CFposition();
0622             unsigned int tofBin = linkiter.TofBin();
0623             subDetTofBin theSubDetTofBin = std::make_pair(detid.subdetId(), tofBin);
0624             auto it = SimHitCollMap.find(theSubDetTofBin);
0625             if (it != SimHitCollMap.end()) {
0626               simhitAddr currentAddr = std::make_pair(it->second, currentCFPos);
0627               if (find(simhitCFPos->begin(), simhitCFPos->end(), currentAddr) == simhitCFPos->end()) {
0628                 simhitCFPos->push_back(currentAddr);
0629               }
0630             }
0631           }
0632         }
0633       }
0634     } else {
0635       edm::LogError("TrackerHitAssociator") << "no Pixel cluster reference attached";
0636     }
0637   }
0638 }
0639
0640 std::vector<PSimHit> TrackerHitAssociator::associateMultiRecHit(const SiTrackerMultiRecHit* multirechit) const {
0641   std::vector<const TrackingRecHit*> componenthits = multirechit->recHits();
0642   //        std::vector<PSimHit> assimhits;
0643   int size = multirechit->weights().size(), idmostprobable = 0;
0644
0645   for (int i = 0; i < size; ++i) {
0646     if (multirechit->weight(i) > multirechit->weight(idmostprobable))
0647       idmostprobable = i;
0648   }
0649
0650   return associateHit(*componenthits[idmostprobable]);
0651 }
0652
0653 std::vector<SimHitIdpr> TrackerHitAssociator::associateMultiRecHitId(const SiTrackerMultiRecHit* multirechit,
0654                                                                      std::vector<simhitAddr>* simhitCFPos) const {
0655   std::vector<const TrackingRecHit*> componenthits = multirechit->recHits();
0656   int size = multirechit->weights().size(), idmostprobable = 0;
0657
0658   for (int i = 0; i < size; ++i) {
0659     if (multirechit->weight(i) > multirechit->weight(idmostprobable))
0660       idmostprobable = i;
0661   }
0662
0663   std::vector<SimHitIdpr> simhitid;
0664   associateHitId(*componenthits[idmostprobable], simhitid, simhitCFPos);
0665   return simhitid;
0666 }
0667
0668 // fastsim
0669 std::vector<SimHitIdpr> TrackerHitAssociator::associateFastRecHit(const FastTrackerRecHit* rechit) const {
0670   vector<SimHitIdpr> simtrackid;
0671   simtrackid.clear();
0672   for (size_t index = 0, indexEnd = rechit->nSimTrackIds(); index < indexEnd; ++index) {
0673     SimHitIdpr currentId(rechit->simTrackId(index), EncodedEventId(rechit->simTrackEventId(index)));
0674     simtrackid.push_back(currentId);
0675   }
0676   return simtrackid;
0677 }
0678
0679 inline std::string TrackerHitAssociator::printDetBnchEvtTrk(const DetId& detid,
0680                                                             const uint32_t& detID,
0681                                                             std::vector<SimHitIdpr>& simtrackid) const {
0682   std::stringstream message;
0683   message << "recHit subdet, detID = " << detid.subdetId() << ", " << detID << ", (bnch, evt, trk) = ";
0684   for (size_t i = 0; i < simtrackid.size(); ++i)
0685     message << ", (" << simtrackid[i].second.bunchCrossing() << ", " << simtrackid[i].second.event() << ", "
0686             << simtrackid[i].first << ")";
0687   // message << std::endl;
0688   return message.str();
0689 }