TrackerHitAssociation/plugins/ClusterTPAssociationProducer.cc

0001 #include <memory>
0002 #include <vector>
0003 #include <utility>
0004
0005 #include "FWCore/Framework/interface/Frameworkfwd.h"
0006 #include "FWCore/Framework/interface/global/EDProducer.h"
0007 #include "FWCore/Framework/interface/Event.h"
0008 #include "FWCore/Utilities/interface/InputTag.h"
0009 #include "FWCore/Utilities/interface/EDGetToken.h"
0010 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0011 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0012 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0013
0014 #include "DataFormats/Common/interface/Handle.h"
0015 #include "DataFormats/Common/interface/DetSetVector.h"
0016 #include "DataFormats/Common/interface/DetSetVectorNew.h"
0017 #include "DataFormats/DetId/interface/DetId.h"
0018 #include "DataFormats/SiPixelDetId/interface/PixelChannelIdentifier.h"
0019 #include "DataFormats/TrackerRecHit2D/interface/OmniClusterRef.h"
0020 #include "DataFormats/SiPixelCluster/interface/SiPixelCluster.h"
0021 #include "DataFormats/SiStripCluster/interface/SiStripCluster.h"
0022 #include "DataFormats/Phase2TrackerCluster/interface/Phase2TrackerCluster1D.h"
0023
0024 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
0025 #include "SimDataFormats/TrackerDigiSimLink/interface/StripDigiSimLink.h"
0026 #include "SimDataFormats/TrackerDigiSimLink/interface/PixelDigiSimLink.h"
0027 #include "DataFormats/Phase2TrackerDigi/interface/Phase2TrackerDigi.h"
0028 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
0029 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticleFwd.h"
0030 #include "SimDataFormats/TrackingAnalysis/interface/UniqueSimTrackId.h"
0031 #include "SimTracker/TrackerHitAssociation/interface/ClusterTPAssociation.h"
0032
0033 namespace {
0034
0035   template <typename T>
0036   void getSimTrackId(std::vector<UniqueSimTrackId>& simTrkId,
0037                      const edm::Handle<edm::DetSetVector<T> >& simLinks,
0038                      const DetId& detId,
0039                      uint32_t channel) {
0040     auto isearch = simLinks->find(detId);
0041     if (isearch != simLinks->end()) {
0042       // Loop over DigiSimLink in this det unit
0043       edm::DetSet<T> link_detset = (*isearch);
0044       for (typename edm::DetSet<T>::const_iterator it = link_detset.data.begin(); it != link_detset.data.end(); ++it) {
0045         if (channel == it->channel()) {
0046           simTrkId.emplace_back(it->SimTrackId(), it->eventId());
0047         }
0048       }
0049     }
0050   }
0051
0052 }  // namespace
0053
0054 class ClusterTPAssociationProducer final : public edm::global::EDProducer<> {
0055 public:
0056   using OmniClusterCollection = std::vector<OmniClusterRef>;
0057
0058   explicit ClusterTPAssociationProducer(const edm::ParameterSet&);
0059   ~ClusterTPAssociationProducer() override = default;
0060
0061   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0062
0063 private:
0064   void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0065
0066   edm::EDGetTokenT<edm::DetSetVector<PixelDigiSimLink> > sipixelSimLinksToken_;
0067   edm::EDGetTokenT<edm::DetSetVector<StripDigiSimLink> > sistripSimLinksToken_;
0068   edm::EDGetTokenT<edm::DetSetVector<PixelDigiSimLink> > siphase2OTSimLinksToken_;
0069   edm::EDGetTokenT<edmNew::DetSetVector<SiPixelCluster> > pixelClustersToken_;
0070   edm::EDGetTokenT<edmNew::DetSetVector<SiStripCluster> > stripClustersToken_;
0071   edm::EDGetTokenT<edmNew::DetSetVector<Phase2TrackerCluster1D> > phase2OTClustersToken_;
0072   edm::EDGetTokenT<TrackingParticleCollection> trackingParticleToken_;
0073   bool throwOnMissingCollections_;
0074 };
0075
0076 ClusterTPAssociationProducer::ClusterTPAssociationProducer(const edm::ParameterSet& cfg)
0077     : sipixelSimLinksToken_(
0078           consumes<edm::DetSetVector<PixelDigiSimLink> >(cfg.getParameter<edm::InputTag>("pixelSimLinkSrc"))),
0079       sistripSimLinksToken_(
0080           consumes<edm::DetSetVector<StripDigiSimLink> >(cfg.getParameter<edm::InputTag>("stripSimLinkSrc"))),
0081       siphase2OTSimLinksToken_(
0082           consumes<edm::DetSetVector<PixelDigiSimLink> >(cfg.getParameter<edm::InputTag>("phase2OTSimLinkSrc"))),
0083       pixelClustersToken_(
0084           consumes<edmNew::DetSetVector<SiPixelCluster> >(cfg.getParameter<edm::InputTag>("pixelClusterSrc"))),
0085       stripClustersToken_(
0086           consumes<edmNew::DetSetVector<SiStripCluster> >(cfg.getParameter<edm::InputTag>("stripClusterSrc"))),
0087       phase2OTClustersToken_(consumes<edmNew::DetSetVector<Phase2TrackerCluster1D> >(
0088           cfg.getParameter<edm::InputTag>("phase2OTClusterSrc"))),
0089       trackingParticleToken_(
0090           consumes<TrackingParticleCollection>(cfg.getParameter<edm::InputTag>("trackingParticleSrc"))),
0091       throwOnMissingCollections_(cfg.getParameter<bool>("throwOnMissingCollections")) {
0092   produces<ClusterTPAssociation>();
0093 }
0094
0095 void ClusterTPAssociationProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0096   edm::ParameterSetDescription desc;
0097   desc.add<edm::InputTag>("simTrackSrc", edm::InputTag("g4SimHits"));
0098   desc.add<edm::InputTag>("pixelSimLinkSrc", edm::InputTag("simSiPixelDigis"));
0099   desc.add<edm::InputTag>("stripSimLinkSrc", edm::InputTag("simSiStripDigis"));
0100   desc.add<edm::InputTag>("phase2OTSimLinkSrc", edm::InputTag("simSiPixelDigis", "Tracker"));
0101   desc.add<edm::InputTag>("pixelClusterSrc", edm::InputTag("siPixelClusters"));
0102   desc.add<edm::InputTag>("stripClusterSrc", edm::InputTag("siStripClusters"));
0103   desc.add<edm::InputTag>("phase2OTClusterSrc", edm::InputTag("siPhase2Clusters"));
0104   desc.add<edm::InputTag>("trackingParticleSrc", edm::InputTag("mix", "MergedTrackTruth"));
0105   desc.add<bool>("throwOnMissingCollections", true);
0106   descriptions.add("tpClusterProducerDefault", desc);
0107 }
0108
0109 void ClusterTPAssociationProducer::produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& es) const {
0110   // Pixel Cluster
0111   edm::Handle<edmNew::DetSetVector<SiPixelCluster> > pixelClusters;
0112   bool foundPixelClusters = iEvent.getByToken(pixelClustersToken_, pixelClusters);
0113
0114   // Strip Cluster
0115   edm::Handle<edmNew::DetSetVector<SiStripCluster> > stripClusters;
0116   bool foundStripClusters = iEvent.getByToken(stripClustersToken_, stripClusters);
0117
0118   // Phase2 Cluster
0119   edm::Handle<edmNew::DetSetVector<Phase2TrackerCluster1D> > phase2OTClusters;
0120   bool foundPhase2OTClusters = iEvent.getByToken(phase2OTClustersToken_, phase2OTClusters);
0121
0122   // Pixel DigiSimLink
0123   edm::Handle<edm::DetSetVector<PixelDigiSimLink> > sipixelSimLinks;
0124   auto pixelSimLinksFound = iEvent.getByToken(sipixelSimLinksToken_, sipixelSimLinks);
0125   if (not throwOnMissingCollections_ and foundPixelClusters and not pixelSimLinksFound) {
0126     auto clusterTPList = std::make_unique<ClusterTPAssociation>();
0127     iEvent.put(std::move(clusterTPList));
0128     return;
0129   }
0130
0131   // SiStrip DigiSimLink
0132   edm::Handle<edm::DetSetVector<StripDigiSimLink> > sistripSimLinks;
0133   auto stripSimLinksFound = iEvent.getByToken(sistripSimLinksToken_, sistripSimLinks);
0134   if (not throwOnMissingCollections_ and foundStripClusters and not stripSimLinksFound) {
0135     auto clusterTPList = std::make_unique<ClusterTPAssociation>();
0136     iEvent.put(std::move(clusterTPList));
0137     return;
0138   }
0139
0140   // Phase2 OT DigiSimLink
0141   edm::Handle<edm::DetSetVector<PixelDigiSimLink> > siphase2OTSimLinks;
0142   auto phase2OTSimLinksFound = iEvent.getByToken(siphase2OTSimLinksToken_, siphase2OTSimLinks);
0143   if (not throwOnMissingCollections_ and foundPhase2OTClusters and not phase2OTSimLinksFound) {
0144     auto clusterTPList = std::make_unique<ClusterTPAssociation>();
0145     iEvent.put(std::move(clusterTPList));
0146     return;
0147   }
0148
0149   // TrackingParticle
0150   edm::Handle<TrackingParticleCollection> TPCollectionH;
0151   auto tpFound = iEvent.getByToken(trackingParticleToken_, TPCollectionH);
0152   if (not throwOnMissingCollections_ and not tpFound) {
0153     auto clusterTPList = std::make_unique<ClusterTPAssociation>();
0154     iEvent.put(std::move(clusterTPList));
0155     return;
0156   }
0157
0158   auto clusterTPList = std::make_unique<ClusterTPAssociation>(TPCollectionH);
0159
0160   // prepare temporary map between SimTrackId and TrackingParticle index
0161   std::unordered_map<UniqueSimTrackId, TrackingParticleRef, UniqueSimTrackIdHash> mapping;
0162   auto const& tpColl = *TPCollectionH.product();
0163   for (TrackingParticleCollection::size_type itp = 0; itp < tpColl.size(); ++itp) {
0164     TrackingParticleRef trackingParticleRef(TPCollectionH, itp);
0165     auto const& trackingParticle = tpColl[itp];
0166     // SimTracks inside TrackingParticle
0167     EncodedEventId eid(trackingParticle.eventId());
0168     //size_t index = 0;
0169     for (auto const& trk : trackingParticle.g4Tracks()) {
0170       UniqueSimTrackId trkid(trk.trackId(), eid);
0171       //std::cout << "creating map for id: " << trkid.first << " with tp: " << trackingParticle.key() << std::endl;
0172       mapping.insert(std::make_pair(trkid, trackingParticleRef));
0173     }
0174   }
0175
0176   std::unordered_set<UniqueSimTrackId, UniqueSimTrackIdHash> simTkIds;
0177   std::vector<UniqueSimTrackId> trkid;
0178   if (foundPixelClusters) {
0179     // Pixel Clusters
0180     clusterTPList->addKeyID(pixelClusters.id());
0181     for (edmNew::DetSetVector<SiPixelCluster>::const_iterator iter = pixelClusters->begin();
0182          iter != pixelClusters->end();
0183          ++iter) {
0184       uint32_t detid = iter->id();
0185       DetId detId(detid);
0186       edmNew::DetSet<SiPixelCluster> link_pixel = (*iter);
0187       for (edmNew::DetSet<SiPixelCluster>::const_iterator di = link_pixel.begin(); di != link_pixel.end(); ++di) {
0188         const SiPixelCluster& cluster = (*di);
0189         edm::Ref<edmNew::DetSetVector<SiPixelCluster>, SiPixelCluster> c_ref = edmNew::makeRefTo(pixelClusters, di);
0190
0191         simTkIds.clear();
0192         for (int irow = cluster.minPixelRow(); irow <= cluster.maxPixelRow(); ++irow) {
0193           for (int icol = cluster.minPixelCol(); icol <= cluster.maxPixelCol(); ++icol) {
0194             uint32_t channel = PixelChannelIdentifier::pixelToChannel(irow, icol);
0195             trkid.clear();
0196             getSimTrackId<PixelDigiSimLink>(trkid, sipixelSimLinks, detId, channel);
0197             simTkIds.insert(trkid.begin(), trkid.end());
0198           }
0199         }
0200         for (auto iset = simTkIds.begin(); iset != simTkIds.end(); iset++) {
0201           auto ipos = mapping.find(*iset);
0202           if (ipos != mapping.end()) {
0203             //std::cout << "cluster in detid: " << detid << " from tp: " << ipos->second.key() << " " << iset->first << std::endl;
0204             clusterTPList->emplace_back(OmniClusterRef(c_ref), ipos->second);
0205           }
0206         }
0207       }
0208     }
0209   }
0210
0211   if (foundStripClusters) {
0212     // Strip Clusters
0213     clusterTPList->addKeyID(stripClusters.id());
0214     for (edmNew::DetSetVector<SiStripCluster>::const_iterator iter = stripClusters->begin(false),
0215                                                               eter = stripClusters->end(false);
0216          iter != eter;
0217          ++iter) {
0218       if (!(*iter).isValid())
0219         continue;
0220       uint32_t detid = iter->id();
0221       DetId detId(detid);
0222       edmNew::DetSet<SiStripCluster> link_strip = (*iter);
0223       for (edmNew::DetSet<SiStripCluster>::const_iterator di = link_strip.begin(); di != link_strip.end(); di++) {
0224         const SiStripCluster& cluster = (*di);
0225         edm::Ref<edmNew::DetSetVector<SiStripCluster>, SiStripCluster> c_ref = edmNew::makeRefTo(stripClusters, di);
0226
0227         simTkIds.clear();
0228         int first = cluster.firstStrip();
0229         int last = first + cluster.amplitudes().size();
0230
0231         for (int istr = first; istr < last; ++istr) {
0232           trkid.clear();
0233           getSimTrackId<StripDigiSimLink>(trkid, sistripSimLinks, detId, istr);
0234           simTkIds.insert(trkid.begin(), trkid.end());
0235         }
0236         for (auto iset = simTkIds.begin(); iset != simTkIds.end(); iset++) {
0237           auto ipos = mapping.find(*iset);
0238           if (ipos != mapping.end()) {
0239             //std::cout << "cluster in detid: " << detid << " from tp: " << ipos->second.key() << " " << iset->first << std::endl;
0240             clusterTPList->emplace_back(OmniClusterRef(c_ref), ipos->second);
0241           }
0242         }
0243       }
0244     }
0245   }
0246
0247   if (foundPhase2OTClusters) {
0248     // Phase2 Clusters
0249     clusterTPList->addKeyID(phase2OTClusters.id());
0250     if (phase2OTClusters.isValid()) {
0251       for (edmNew::DetSetVector<Phase2TrackerCluster1D>::const_iterator iter = phase2OTClusters->begin(false),
0252                                                                         eter = phase2OTClusters->end(false);
0253            iter != eter;
0254            ++iter) {
0255         if (!(*iter).isValid())
0256           continue;
0257         uint32_t detid = iter->id();
0258         DetId detId(detid);
0259         edmNew::DetSet<Phase2TrackerCluster1D> link_phase2 = (*iter);
0260         for (edmNew::DetSet<Phase2TrackerCluster1D>::const_iterator di = link_phase2.begin(); di != link_phase2.end();
0261              di++) {
0262           const Phase2TrackerCluster1D& cluster = (*di);
0263           edm::Ref<edmNew::DetSetVector<Phase2TrackerCluster1D>, Phase2TrackerCluster1D> c_ref =
0264               edmNew::makeRefTo(phase2OTClusters, di);
0265
0266           simTkIds.clear();
0267
0268           for (unsigned int istr(0); istr < cluster.size(); ++istr) {
0269             uint32_t channel = Phase2TrackerDigi::pixelToChannel(cluster.firstRow() + istr, cluster.column());
0270             trkid.clear();
0271             getSimTrackId<PixelDigiSimLink>(trkid, siphase2OTSimLinks, detId, channel);
0272             simTkIds.insert(trkid.begin(), trkid.end());
0273           }
0274
0275           for (auto iset = simTkIds.begin(); iset != simTkIds.end(); iset++) {
0276             auto ipos = mapping.find(*iset);
0277             if (ipos != mapping.end()) {
0278               clusterTPList->emplace_back(OmniClusterRef(c_ref), ipos->second);
0279             }
0280           }
0281         }
0282       }
0283     }
0284   }
0285   clusterTPList->sortAndUnique();
0286   iEvent.put(std::move(clusterTPList));
0287 }
0288
0289 #include "FWCore/PluginManager/interface/ModuleDef.h"
0290 #include "FWCore/Framework/interface/MakerMacros.h"
0291
0292 DEFINE_FWK_MODULE(ClusterTPAssociationProducer);