TrackerHitAssociation/plugins/SimDoubletsProducer.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    SimTracker/TrackerHitAssociation
0004 // Class:      SimDoubletsProducer
0005 //
0006
0007 // user include files
0008 #include "FWCore/Framework/interface/MakerMacros.h"
0009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0010
0011 #include "FWCore/Framework/interface/Frameworkfwd.h"
0012 #include "FWCore/Framework/interface/Event.h"
0013 #include "FWCore/Framework/interface/stream/EDProducer.h"
0014 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0015 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0016 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0017 #include "FWCore/Framework/interface/EventSetup.h"
0018 #include "FWCore/Utilities/interface/InputTag.h"
0019 #include "FWCore/Utilities/interface/IndexSet.h"
0020 #include "FWCore/Utilities/interface/ESGetToken.h"
0021
0022 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0023 #include "Geometry/CommonTopologies/interface/SimplePixelTopology.h"
0024
0025 #include "SimTracker/Common/interface/TrackingParticleSelector.h"
0026 #include "SimTracker/TrackerHitAssociation/interface/ClusterTPAssociation.h"
0027
0028 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0029 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHitCollection.h"
0030 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0031 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0032 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0033 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
0034 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHit.h"
0035 #include "DataFormats/TrackerRecHit2D/interface/OmniClusterRef.h"
0036 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
0037 #include "SimDataFormats/TrackingAnalysis/interface/SimDoublets.h"
0038
0039 #include <cstddef>
0040 #include <utility>
0041 #include <vector>
0042 #include <memory>
0043 #include <typeinfo>
0044
0045 /** Class: SimDoubletsProducer
0046  *
0047  * @brief Produces SimDoublets (MC-info based PixelRecHit doublets) for selected TrackingParticles.
0048  *
0049  * SimDoublets represent the true doublets of RecHits that a simulated particle (TrackingParticle)
0050  * created in the pixel detector. They can be used to analyze cuts which are applied in the reconstruction
0051  * when producing doublets as the first part of patatrack pixel tracking.
0052  *
0053  * The SimDoublets are produced in the following way:
0054  * 1. We select reasonable TrackingParticles according to the criteria given in the config file as
0055  *    "TrackingParticleSelectionConfig".
0056  * 2. For each selected particle, we create and append a new SimDoublets object to the SimDoubletsCollection.
0057  * 3. We loop over all RecHits in the pixel tracker and check if the given RecHit is associated to one of
0058  *    the selected particles (association via TP to cluster association). If it is, we add a RecHit reference
0059  *    to the respective SimDoublet.
0060  * 4. In the end, we sort the RecHits in each SimDoublets object according to their global position.
0061  *
0062  * @author Jan Schulz (jan.gerrit.schulz@cern.ch)
0063  * @date January 2025
0064  */
0065 template <typename TrackerTraits>
0066 class SimDoubletsProducer : public edm::stream::EDProducer<> {
0067 public:
0068   explicit SimDoubletsProducer(const edm::ParameterSet&);
0069   static void fillDescriptions(edm::ConfigurationDescriptions&);
0070
0071   void produce(edm::Event&, const edm::EventSetup&) override;
0072   void beginRun(const edm::Run&, const edm::EventSetup&) override;
0073
0074 private:
0075   TrackingParticleSelector trackingParticleSelector;
0076   const TrackerTopology* trackerTopology_ = nullptr;
0077
0078   const edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> topology_getToken_;
0079   const edm::EDGetTokenT<ClusterTPAssociation> clusterTPAssociation_getToken_;
0080   const edm::EDGetTokenT<TrackingParticleCollection> trackingParticles_getToken_;
0081   const edm::EDGetTokenT<SiPixelRecHitCollection> pixelRecHits_getToken_;
0082   const edm::EDGetTokenT<reco::BeamSpot> beamSpot_getToken_;
0083   const edm::EDPutTokenT<SimDoubletsCollection> simDoublets_putToken_;
0084 };
0085
0086 namespace simdoublets {
0087   // function that determines the layerId from the detId for Phase 1 and 2
0088   template <typename TrackerTraits>
0089   unsigned int getLayerId(DetId const& detId, const TrackerTopology* trackerTopology) {
0090     // number of barrel layers
0091     constexpr unsigned int numBarrelLayers{4};
0092     // number of disks per endcap
0093     constexpr unsigned int numEndcapDisks = (TrackerTraits::numberOfLayers - numBarrelLayers) / 2;
0094
0095     // set default to 999 (invalid)
0096     unsigned int layerId{999};
0097
0098     if (detId.subdetId() == PixelSubdetector::PixelBarrel) {
0099       // subtract 1 in the barrel to get, e.g. for Phase 2, from (1,4) to (0,3)
0100       layerId = trackerTopology->pxbLayer(detId) - 1;
0101     } else if (detId.subdetId() == PixelSubdetector::PixelEndcap) {
0102       if (trackerTopology->pxfSide(detId) == 1) {
0103         // add offset in the backward endcap to get, e.g. for Phase 2, from (1,12) to (16,27)
0104         layerId = trackerTopology->pxfDisk(detId) + numBarrelLayers + numEndcapDisks - 1;
0105       } else {
0106         // add offest in the forward endcap to get, e.g. for Phase 2, from (1,12) to (4,15)
0107         layerId = trackerTopology->pxfDisk(detId) + numBarrelLayers - 1;
0108       }
0109     }
0110     // return the determined Id
0111     return layerId;
0112   }
0113 }  // namespace simdoublets
0114
0115 // constructor
0116 template <typename TrackerTraits>
0117 SimDoubletsProducer<TrackerTraits>::SimDoubletsProducer(const edm::ParameterSet& pSet)
0118     : topology_getToken_(esConsumes<TrackerTopology, TrackerTopologyRcd>()),
0119       clusterTPAssociation_getToken_(
0120           consumes<ClusterTPAssociation>(pSet.getParameter<edm::InputTag>("clusterTPAssociationSrc"))),
0121       trackingParticles_getToken_(consumes(pSet.getParameter<edm::InputTag>("trackingParticleSrc"))),
0122       pixelRecHits_getToken_(consumes(pSet.getParameter<edm::InputTag>("pixelRecHitSrc"))),
0123       beamSpot_getToken_(consumes(pSet.getParameter<edm::InputTag>("beamSpotSrc"))),
0124       simDoublets_putToken_(produces<SimDoubletsCollection>()) {
0125   // initialize the selector for TrackingParticles used to create SimDoublets
0126   const edm::ParameterSet& pSetTPSel = pSet.getParameter<edm::ParameterSet>("TrackingParticleSelectionConfig");
0127   trackingParticleSelector = TrackingParticleSelector(pSetTPSel.getParameter<double>("ptMin"),
0128                                                       pSetTPSel.getParameter<double>("ptMax"),
0129                                                       pSetTPSel.getParameter<double>("minRapidity"),
0130                                                       pSetTPSel.getParameter<double>("maxRapidity"),
0131                                                       pSetTPSel.getParameter<double>("tip"),
0132                                                       pSetTPSel.getParameter<double>("lip"),
0133                                                       pSetTPSel.getParameter<int>("minHit"),
0134                                                       pSetTPSel.getParameter<bool>("signalOnly"),
0135                                                       pSetTPSel.getParameter<bool>("intimeOnly"),
0136                                                       pSetTPSel.getParameter<bool>("chargedOnly"),
0137                                                       pSetTPSel.getParameter<bool>("stableOnly"),
0138                                                       pSetTPSel.getParameter<std::vector<int>>("pdgId"),
0139                                                       pSetTPSel.getParameter<bool>("invertRapidityCut"),
0140                                                       pSetTPSel.getParameter<double>("minPhi"),
0141                                                       pSetTPSel.getParameter<double>("maxPhi"));
0142 }
0143
0144 // dummy fillDescription
0145 template <typename TrackerTraits>
0146 void SimDoubletsProducer<TrackerTraits>::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {}
0147
0148 // Phase 1 fillDescription
0149 template <>
0150 void SimDoubletsProducer<pixelTopology::Phase1>::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0151   edm::ParameterSetDescription desc;
0152
0153   // sources for cluster-TrackingParticle association, TrackingParticles, RecHits and the beamspot
0154   desc.add<edm::InputTag>("clusterTPAssociationSrc", edm::InputTag("hltTPClusterProducer"));
0155   desc.add<edm::InputTag>("trackingParticleSrc", edm::InputTag("mix", "MergedTrackTruth"));
0156   desc.add<edm::InputTag>("pixelRecHitSrc", edm::InputTag("hltSiPixelRecHits"));
0157   desc.add<edm::InputTag>("beamSpotSrc", edm::InputTag("hltOnlineBeamSpot"));
0158
0159   // parameter set for the selection of TrackingParticles that will be used for SimHitDoublets
0160   edm::ParameterSetDescription descTPSelector;
0161   descTPSelector.add<double>("ptMin", 0.9);
0162   descTPSelector.add<double>("ptMax", 1e100);
0163   descTPSelector.add<double>("minRapidity", -3.);
0164   descTPSelector.add<double>("maxRapidity", 3.);
0165   descTPSelector.add<double>("tip", 60.);
0166   descTPSelector.add<double>("lip", 30.);
0167   descTPSelector.add<int>("minHit", 0);
0168   descTPSelector.add<bool>("signalOnly", true);
0169   descTPSelector.add<bool>("intimeOnly", false);
0170   descTPSelector.add<bool>("chargedOnly", true);
0171   descTPSelector.add<bool>("stableOnly", false);
0172   descTPSelector.add<std::vector<int>>("pdgId", {});
0173   descTPSelector.add<bool>("invertRapidityCut", false);
0174   descTPSelector.add<double>("minPhi", -3.2);
0175   descTPSelector.add<double>("maxPhi", 3.2);
0176   desc.add<edm::ParameterSetDescription>("TrackingParticleSelectionConfig", descTPSelector);
0177
0178   descriptions.addWithDefaultLabel(desc);
0179 }
0180
0181 // Phase 2 fillDescription
0182 template <>
0183 void SimDoubletsProducer<pixelTopology::Phase2>::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0184   edm::ParameterSetDescription desc;
0185
0186   // sources for cluster-TrackingParticle association, TrackingParticles, RecHits and the beamspot
0187   desc.add<edm::InputTag>("clusterTPAssociationSrc", edm::InputTag("hltTPClusterProducer"));
0188   desc.add<edm::InputTag>("trackingParticleSrc", edm::InputTag("mix", "MergedTrackTruth"));
0189   desc.add<edm::InputTag>("pixelRecHitSrc", edm::InputTag("hltSiPixelRecHits"));
0190   desc.add<edm::InputTag>("beamSpotSrc", edm::InputTag("hltOnlineBeamSpot"));
0191
0192   // parameter set for the selection of TrackingParticles that will be used for SimHitDoublets
0193   edm::ParameterSetDescription descTPSelector;
0194   descTPSelector.add<double>("ptMin", 0.9);
0195   descTPSelector.add<double>("ptMax", 1e100);
0196   descTPSelector.add<double>("minRapidity", -4.5);
0197   descTPSelector.add<double>("maxRapidity", 4.5);
0198   descTPSelector.add<double>("tip", 2.);  // NOTE: differs from HLT MultiTrackValidator
0199   descTPSelector.add<double>("lip", 30.);
0200   descTPSelector.add<int>("minHit", 0);
0201   descTPSelector.add<bool>("signalOnly", true);
0202   descTPSelector.add<bool>("intimeOnly", false);
0203   descTPSelector.add<bool>("chargedOnly", true);
0204   descTPSelector.add<bool>("stableOnly", false);
0205   descTPSelector.add<std::vector<int>>("pdgId", {});
0206   descTPSelector.add<bool>("invertRapidityCut", false);
0207   descTPSelector.add<double>("minPhi", -3.2);
0208   descTPSelector.add<double>("maxPhi", 3.2);
0209   desc.add<edm::ParameterSetDescription>("TrackingParticleSelectionConfig", descTPSelector);
0210
0211   descriptions.addWithDefaultLabel(desc);
0212 }
0213
0214 template <typename TrackerTraits>
0215 void SimDoubletsProducer<TrackerTraits>::beginRun(const edm::Run& run, const edm::EventSetup& eventSetup) {}
0216
0217 template <typename TrackerTraits>
0218 void SimDoubletsProducer<TrackerTraits>::produce(edm::Event& event, const edm::EventSetup& eventSetup) {
0219   // get TrackerTopology
0220   trackerTopology_ = &eventSetup.getData(topology_getToken_);
0221
0222   // get cluster to TrackingParticle association
0223   ClusterTPAssociation const& clusterTPAssociation = event.get(clusterTPAssociation_getToken_);
0224
0225   // get the pixel RecHit collection from the event
0226   edm::Handle<SiPixelRecHitCollection> hits;
0227   event.getByToken(pixelRecHits_getToken_, hits);
0228   if (!hits.isValid()) {
0229     return;
0230   }
0231
0232   // get TrackingParticles from the event
0233   edm::Handle<TrackingParticleCollection> trackingParticles;
0234   event.getByToken(trackingParticles_getToken_, trackingParticles);
0235   if (!trackingParticles.isValid()) {
0236     return;
0237   }
0238
0239   // get beamspot from the event
0240   edm::Handle<reco::BeamSpot> beamSpot;
0241   event.getByToken(beamSpot_getToken_, beamSpot);
0242   if (!beamSpot.isValid()) {
0243     return;
0244   }
0245
0246   // create collection of SimDoublets
0247   // each element will correspond to one selected TrackingParticle
0248   SimDoubletsCollection simDoubletsCollection;
0249
0250   // loop over TrackingParticles
0251   for (size_t i = 0; i < trackingParticles->size(); ++i) {
0252     TrackingParticle const& trackingParticle = trackingParticles->at(i);
0253
0254     // select reasonable TrackingParticles for the study (e.g., only signal)
0255     if (trackingParticleSelector(trackingParticle)) {
0256       simDoubletsCollection.push_back(SimDoublets(TrackingParticleRef(trackingParticles, i), *beamSpot));
0257     }
0258   }
0259
0260   // create a set of the keys of the selected TrackingParticles
0261   edm::IndexSet selectedTrackingParticleKeys;
0262   selectedTrackingParticleKeys.reserve(simDoubletsCollection.size());
0263   for (const auto& simDoublets : simDoubletsCollection) {
0264     TrackingParticleRef trackingParticleRef = simDoublets.trackingParticle();
0265     selectedTrackingParticleKeys.insert(trackingParticleRef.key());
0266   }
0267
0268   // initialize a couple of counters
0269   int count_associatedRecHits{0}, count_RecHitsInSimDoublets{0};
0270
0271   // loop over pixel RecHit collections of the different pixel modules
0272   for (const auto& detSet : *hits) {
0273     // determine layer Id
0274     unsigned int layerId = simdoublets::getLayerId<TrackerTraits>(detSet.detId(), trackerTopology_);
0275
0276     // loop over RecHits
0277     for (auto const& hit : detSet) {
0278       // find associated TrackingParticles
0279       auto range = clusterTPAssociation.equal_range(OmniClusterRef(hit.cluster()));
0280
0281       // if the RecHit has associated TrackingParticles
0282       if (range.first != range.second) {
0283         for (auto assocTrackingParticleIter = range.first; assocTrackingParticleIter != range.second;
0284              assocTrackingParticleIter++) {
0285           const TrackingParticleRef assocTrackingParticle = (assocTrackingParticleIter->second);
0286
0287           // if the associated TrackingParticle is among the selected ones
0288           if (selectedTrackingParticleKeys.has(assocTrackingParticle.key())) {
0289             SiPixelRecHitRef hitRef = edmNew::makeRefTo(hits, &hit);
0290             count_associatedRecHits++;
0291             // loop over collection of SimDoublets and find the one of the associated TrackingParticle
0292             for (auto& simDoublets : simDoubletsCollection) {
0293               TrackingParticleRef trackingParticleRef = simDoublets.trackingParticle();
0294               if (assocTrackingParticle.key() == trackingParticleRef.key()) {
0295                 simDoublets.addRecHit(hitRef, layerId);
0296                 count_RecHitsInSimDoublets++;
0297               }
0298             }
0299           }
0300         }
0301       }
0302     }  // end loop over RecHits
0303   }  // end loop over pixel RecHit collections of the different pixel modules
0304
0305   // loop over collection of SimDoublets and sort the RecHits according to their position
0306   for (auto& simDoublets : simDoubletsCollection) {
0307     simDoublets.sortRecHits();
0308   }
0309
0310   LogDebug("SimDoubletsProducer") << "Size of SiPixelRecHitCollection : " << hits->size() << std::endl;
0311   LogDebug("SimDoubletsProducer") << count_associatedRecHits << " of " << hits->size()
0312                                   << " RecHits are associated to selected TrackingParticles ("
0313                                   << count_RecHitsInSimDoublets - count_associatedRecHits
0314                                   << " of them were associated multiple times)." << std::endl;
0315   LogDebug("SimDoubletsProducer") << "Number of selected TrackingParticles : " << simDoubletsCollection.size()
0316                                   << std::endl;
0317   LogDebug("SimDoubletsProducer") << "Size of TrackingParticle Collection  : " << trackingParticles->size()
0318                                   << std::endl;
0319
0320   // put the produced SimDoublets collection in the event
0321   event.emplace(simDoublets_putToken_, std::move(simDoubletsCollection));
0322 }
0323
0324 // define two plugins for Phase 1 and 2
0325 using SimDoubletsProducerPhase1 = SimDoubletsProducer<pixelTopology::Phase1>;
0326 using SimDoubletsProducerPhase2 = SimDoubletsProducer<pixelTopology::Phase2>;
0327
0328 DEFINE_FWK_MODULE(SimDoubletsProducerPhase1);
0329 DEFINE_FWK_MODULE(SimDoubletsProducerPhase2);