TrackingAnalysis/src/TrackingParticleNumberOfLayers.cc

0001 #include "SimGeneral/TrackingAnalysis/interface/TrackingParticleNumberOfLayers.h"
0002
0003 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0004 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
0005 #include "FWCore/Utilities/interface/isFinite.h"
0006
0007 namespace {
0008   bool trackIdHitPairLess(const std::pair<unsigned int, const PSimHit *> &a,
0009                           const std::pair<unsigned int, const PSimHit *> &b) {
0010     return a.first < b.first;
0011   }
0012
0013   bool trackIdHitPairLessSort(const std::pair<unsigned int, const PSimHit *> &a,
0014                               const std::pair<unsigned int, const PSimHit *> &b) {
0015     if (a.first == b.first) {
0016       const auto atof = edm::isFinite(a.second->timeOfFlight())
0017                             ? a.second->timeOfFlight()
0018                             : std::numeric_limits<decltype(a.second->timeOfFlight())>::max();
0019       const auto btof = edm::isFinite(b.second->timeOfFlight())
0020                             ? b.second->timeOfFlight()
0021                             : std::numeric_limits<decltype(b.second->timeOfFlight())>::max();
0022       return atof < btof;
0023     }
0024     return a.first < b.first;
0025   }
0026 }  // namespace
0027
0028 TrackingParticleNumberOfLayers::TrackingParticleNumberOfLayers(
0029     const edm::Event &iEvent, const std::vector<edm::EDGetTokenT<std::vector<PSimHit>>> &simHitTokens) {
0030   // A multimap linking SimTrack::trackId() to a pointer to PSimHit
0031   // Similar to TrackingTruthAccumulator
0032   for (const auto &simHitToken : simHitTokens) {
0033     edm::Handle<std::vector<PSimHit>> hsimhits;
0034     iEvent.getByToken(simHitToken, hsimhits);
0035     trackIdToHitPtr_.reserve(trackIdToHitPtr_.size() + hsimhits->size());
0036     for (const auto &simHit : *hsimhits) {
0037       trackIdToHitPtr_.emplace_back(simHit.trackId(), &simHit);
0038     }
0039   }
0040   std::stable_sort(trackIdToHitPtr_.begin(), trackIdToHitPtr_.end(), trackIdHitPairLessSort);
0041 }
0042
0043 std::tuple<std::unique_ptr<edm::ValueMap<unsigned int>>,
0044            std::unique_ptr<edm::ValueMap<unsigned int>>,
0045            std::unique_ptr<edm::ValueMap<unsigned int>>>
0046 TrackingParticleNumberOfLayers::calculate(const edm::Handle<TrackingParticleCollection> &htps,
0047                                           const TrackerTopology &tTopo) const {
0048   const TrackingParticleCollection &tps = *htps;
0049   std::vector<unsigned int> valuesLayers(tps.size(), 0);
0050   std::vector<unsigned int> valuesPixelLayers(tps.size(), 0);
0051   std::vector<unsigned int> valuesStripMonoAndStereoLayers(tps.size(), 0);
0052   for (size_t iTP = 0; iTP < tps.size(); ++iTP) {
0053     const TrackingParticle &tp = tps[iTP];
0054     const auto pdgId = tp.pdgId();
0055
0056     // I would prefer a better way...
0057     constexpr unsigned int maxSubdet = static_cast<unsigned>(StripSubdetector::TEC) + 1;
0058     constexpr unsigned int maxLayer = 0xF + 1;  // as in HitPattern.h
0059     bool hasHit[maxSubdet][maxLayer];
0060     bool hasPixel[maxSubdet][maxLayer];
0061     bool hasMono[maxSubdet][maxLayer];
0062     bool hasStereo[maxSubdet][maxLayer];
0063     memset(hasHit, 0, sizeof(hasHit));
0064     memset(hasPixel, 0, sizeof(hasPixel));
0065     memset(hasMono, 0, sizeof(hasMono));
0066     memset(hasStereo, 0, sizeof(hasStereo));
0067
0068     for (const SimTrack &simTrack : tp.g4Tracks()) {
0069       // Logic is from TrackingTruthAccumulator
0070       auto range = std::equal_range(trackIdToHitPtr_.begin(),
0071                                     trackIdToHitPtr_.end(),
0072                                     std::pair<unsigned int, const PSimHit *>(simTrack.trackId(), nullptr),
0073                                     trackIdHitPairLess);
0074       if (range.first == range.second)
0075         continue;
0076
0077       auto iHitPtr = range.first;
0078       for (; iHitPtr != range.second; ++iHitPtr) {
0079         // prevent simhits with particleType != pdgId from being the "first hit"
0080         if (iHitPtr->second->particleType() == pdgId)
0081           break;
0082       }
0083       if (iHitPtr == range.second)  // no simhits with particleType == pdgId
0084         continue;
0085       int processType = iHitPtr->second->processType();
0086       int particleType = iHitPtr->second->particleType();
0087
0088       for (; iHitPtr != range.second; ++iHitPtr) {
0089         const PSimHit &simHit = *(iHitPtr->second);
0090         if (simHit.eventId() != tp.eventId())
0091           continue;
0092         DetId newDetector = DetId(simHit.detUnitId());
0093
0094         // Check for delta and interaction products discards
0095         if (processType == simHit.processType() && particleType == simHit.particleType() &&
0096             pdgId == simHit.particleType()) {
0097           // The logic of this piece follows HitPattern
0098           bool isPixel = false;
0099           bool isStripStereo = false;
0100           bool isStripMono = false;
0101
0102           switch (newDetector.subdetId()) {
0103             case PixelSubdetector::PixelBarrel:
0104             case PixelSubdetector::PixelEndcap:
0105               isPixel = true;
0106               break;
0107             case StripSubdetector::TIB:
0108               isStripMono = tTopo.tibIsRPhi(newDetector);
0109               isStripStereo = tTopo.tibIsStereo(newDetector);
0110               break;
0111             case StripSubdetector::TID:
0112               isStripMono = tTopo.tidIsRPhi(newDetector);
0113               isStripStereo = tTopo.tidIsStereo(newDetector);
0114               break;
0115             case StripSubdetector::TOB:
0116               isStripMono = tTopo.tobIsRPhi(newDetector);
0117               isStripStereo = tTopo.tobIsStereo(newDetector);
0118               break;
0119             case StripSubdetector::TEC:
0120               isStripMono = tTopo.tecIsRPhi(newDetector);
0121               isStripStereo = tTopo.tecIsStereo(newDetector);
0122               break;
0123           }
0124
0125           const auto subdet = newDetector.subdetId();
0126           const auto layer = tTopo.layer(newDetector);
0127
0128           hasHit[subdet][layer] = true;
0129           if (isPixel)
0130             hasPixel[subdet][layer] = isPixel;
0131           else if (isStripMono)
0132             hasMono[subdet][layer] = isStripMono;
0133           else if (isStripStereo)
0134             hasStereo[subdet][layer] = isStripStereo;
0135         }
0136       }
0137     }
0138
0139     unsigned int nLayers = 0;
0140     unsigned int nPixelLayers = 0;
0141     unsigned int nStripMonoAndStereoLayers = 0;
0142     for (unsigned int i = 0; i < maxSubdet; ++i) {
0143       for (unsigned int j = 0; j < maxLayer; ++j) {
0144         nLayers += hasHit[i][j];
0145         nPixelLayers += hasPixel[i][j];
0146         nStripMonoAndStereoLayers += (hasMono[i][j] && hasStereo[i][j]);
0147       }
0148     }
0149
0150     valuesLayers[iTP] = nLayers;
0151     valuesPixelLayers[iTP] = nPixelLayers;
0152     valuesStripMonoAndStereoLayers[iTP] = nStripMonoAndStereoLayers;
0153   }
0154
0155   auto ret0 = std::make_unique<edm::ValueMap<unsigned int>>();
0156   {
0157     edm::ValueMap<unsigned int>::Filler filler(*ret0);
0158     filler.insert(htps, valuesLayers.begin(), valuesLayers.end());
0159     filler.fill();
0160   }
0161   auto ret1 = std::make_unique<edm::ValueMap<unsigned int>>();
0162   {
0163     edm::ValueMap<unsigned int>::Filler filler(*ret1);
0164     filler.insert(htps, valuesPixelLayers.begin(), valuesPixelLayers.end());
0165     filler.fill();
0166   }
0167   auto ret2 = std::make_unique<edm::ValueMap<unsigned int>>();
0168   {
0169     edm::ValueMap<unsigned int>::Filler filler(*ret2);
0170     filler.insert(htps, valuesStripMonoAndStereoLayers.begin(), valuesStripMonoAndStereoLayers.end());
0171     filler.fill();
0172   }
0173
0174   return std::make_tuple(std::move(ret0), std::move(ret1), std::move(ret2));
0175 }