SiPixelPhase1TrackingParticleV/src/SiPixelPhase1TrackingParticleV.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:     SiPixelPhase1TrackingParticleV
0004 // Class:       SiPixelPhase1TrackingParticleV
0005 //
0006
0007 // Original Author: Marcel Schneider
0008 // Additional Authors: Alexander Morton - modifying code for validation use
0009
0010 #include "FWCore/Framework/interface/MakerMacros.h"
0011 #include "Validation/SiPixelPhase1TrackingParticleV/interface/SiPixelPhase1TrackingParticleV.h"
0012
0013 #include "DataFormats/GeometryVector/interface/LocalPoint.h"
0014 #include "FWCore/Framework/interface/ESHandle.h"
0015
0016 #include "FWCore/Utilities/interface/isFinite.h"
0017 #include "TrackingTools/TrackAssociator/interface/DetIdAssociator.h"
0018
0019 #include "SimDataFormats/TrackingHit/interface/PSimHit.h"
0020 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
0021
0022 class TrackAssociatorByHits;
0023
0024 namespace {
0025   bool trackIdHitPairLess(const std::pair<unsigned int, const PSimHit *> &a,
0026                           const std::pair<unsigned int, const PSimHit *> &b) {
0027     return a.first < b.first;
0028   }
0029
0030   bool trackIdHitPairLessSort(const std::pair<unsigned int, const PSimHit *> &a,
0031                               const std::pair<unsigned int, const PSimHit *> &b) {
0032     if (a.first == b.first) {
0033       const auto atof = edm::isFinite(a.second->timeOfFlight())
0034                             ? a.second->timeOfFlight()
0035                             : std::numeric_limits<decltype(a.second->timeOfFlight())>::max();
0036       const auto btof = edm::isFinite(b.second->timeOfFlight())
0037                             ? b.second->timeOfFlight()
0038                             : std::numeric_limits<decltype(b.second->timeOfFlight())>::max();
0039       return atof < btof;
0040     }
0041     return a.first < b.first;
0042   }
0043 }  // namespace
0044
0045 SiPixelPhase1TrackingParticleV::SiPixelPhase1TrackingParticleV(const edm::ParameterSet &iConfig)
0046     : SiPixelPhase1Base(iConfig),
0047       vec_TrackingParticle_Token_(consumes<TrackingParticleCollection>(iConfig.getParameter<edm::InputTag>("src"))) {
0048   for (const auto &tag : iConfig.getParameter<std::vector<edm::InputTag>>("simHitToken")) {
0049     simHitTokens_.push_back(consumes<std::vector<PSimHit>>(tag));
0050   }
0051 }
0052
0053 void SiPixelPhase1TrackingParticleV::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0054   edm::Handle<TrackingParticleCollection> TruthTrackContainer;
0055   iEvent.getByToken(vec_TrackingParticle_Token_, TruthTrackContainer);
0056   const TrackingParticleCollection *tPC = TruthTrackContainer.product();
0057   std::vector<std::pair<unsigned int, const PSimHit *>> trackIdToHitPtr;
0058
0059   // A multimap linking SimTrack::trackId() to a pointer to PSimHit
0060   // Similar to TrackingTruthAccumulator
0061   for (const auto &simHitToken : simHitTokens_) {
0062     edm::Handle<std::vector<PSimHit>> hsimhits;
0063     iEvent.getByToken(simHitToken, hsimhits);
0064     trackIdToHitPtr.reserve(trackIdToHitPtr.size() + hsimhits->size());
0065     for (const auto &simHit : *hsimhits) {
0066       trackIdToHitPtr.emplace_back(simHit.trackId(), &simHit);
0067     }
0068   }
0069   std::stable_sort(trackIdToHitPtr.begin(), trackIdToHitPtr.end(), trackIdHitPairLessSort);
0070
0071   // Loop over TrackingParticle's
0072   for (TrackingParticleCollection::const_iterator t = tPC->begin(); t != tPC->end(); ++t) {
0073     // histo manager requires a det ID, use first tracker hit
0074
0075     bool isBpixtrack = false, isFpixtrack = false;
0076     DetId id;
0077
0078     for (const SimTrack &simTrack : t->g4Tracks()) {
0079       // Logic is from TrackingTruthAccumulator
0080       auto range = std::equal_range(trackIdToHitPtr.begin(),
0081                                     trackIdToHitPtr.end(),
0082                                     std::pair<unsigned int, const PSimHit *>(simTrack.trackId(), nullptr),
0083                                     trackIdHitPairLess);
0084       if (range.first == range.second)
0085         continue;
0086
0087       auto iHitPtr = range.first;
0088       for (; iHitPtr != range.second; ++iHitPtr) {
0089         const PSimHit &simHit = *(iHitPtr->second);
0090         if (simHit.eventId() != t->eventId())
0091           continue;
0092         id = DetId(simHit.detUnitId());
0093
0094         // check we are in pixel
0095         uint32_t subdetid = (id.subdetId());
0096         if (subdetid == PixelSubdetector::PixelBarrel)
0097           isBpixtrack = true;
0098         if (subdetid == PixelSubdetector::PixelEndcap)
0099           isFpixtrack = true;
0100         if (subdetid != PixelSubdetector::PixelBarrel && subdetid != PixelSubdetector::PixelEndcap)
0101           continue;
0102       }
0103     }
0104
0105     if (isBpixtrack || isFpixtrack) {
0106       histo[MASS].fill(t->mass(), id, &iEvent);
0107       histo[CHARGE].fill(t->charge(), id, &iEvent);
0108       histo[ID].fill(t->pdgId(), id, &iEvent);
0109       histo[NHITS].fill(t->numberOfTrackerHits(), id, &iEvent);
0110       histo[MATCHED].fill(t->numberOfTrackerLayers(), id, &iEvent);
0111       histo[PT].fill(sqrt(t->momentum().perp2()), id, &iEvent);
0112       histo[PHI].fill(t->momentum().Phi(), id, &iEvent);
0113       histo[ETA].fill(t->momentum().eta(), id, &iEvent);
0114       histo[VTX].fill(t->vx(), id, &iEvent);
0115       histo[VTY].fill(t->vy(), id, &iEvent);
0116       histo[VYZ].fill(t->vz(), id, &iEvent);
0117       histo[TIP].fill(sqrt(t->vertex().perp2()), id, &iEvent);
0118       histo[LIP].fill(t->vz(), id, &iEvent);
0119     }
0120   }
0121 }
0122
0123 DEFINE_FWK_MODULE(SiPixelPhase1TrackingParticleV);