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