SiTrackerPhase2V/plugins/Phase2ITValidateCluster.cc

0001 // -*- C++ -*-
0002 ///bookLayer
0003 // Package:    Phase2ITValidateCluster
0004 // Class:      Phase2ITValidateCluster
0005 //
0006 /**\class Phase2ITValidateCluster Phase2ITValidateCluster.cc
0007
0008  Description: Validation plots tracker clusters.
0009
0010 */
0011 //
0012 // Author: Gabriel Ramirez, Suvankar Roy Chowdhury
0013 // Date: May 23, 2020
0014 //
0015 #include <memory>
0016 #include "FWCore/Framework/interface/Frameworkfwd.h"
0017 #include "FWCore/Framework/interface/Event.h"
0018 #include "FWCore/Framework/interface/ESHandle.h"
0019 #include "FWCore/Framework/interface/MakerMacros.h"
0020 #include "FWCore/Framework/interface/ESWatcher.h"
0021 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0022 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0023 #include "FWCore/ServiceRegistry/interface/Service.h"
0024 #include "FWCore/Utilities/interface/InputTag.h"
0025
0026 #include "DataFormats/Common/interface/Handle.h"
0027 #include "DataFormats/DetId/interface/DetId.h"
0028 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0029 #include "DataFormats/Common/interface/DetSetVectorNew.h"
0030 #include "DataFormats/Common/interface/DetSetVector.h"
0031 #include "DataFormats/SiPixelCluster/interface/SiPixelCluster.h"
0032 #include "DataFormats/SiPixelDetId/interface/PixelChannelIdentifier.h"
0033
0034 #include "SimDataFormats/TrackerDigiSimLink/interface/PixelDigiSimLink.h"
0035 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
0036 #include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
0037 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
0038
0039 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0040 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0041 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0042 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0043
0044 // DQM Histograming
0045 #include "DQMServices/Core/interface/MonitorElement.h"
0046 #include "DQMServices/Core/interface/DQMEDAnalyzer.h"
0047 #include "DQMServices/Core/interface/DQMStore.h"
0048
0049 class Phase2ITValidateCluster : public DQMEDAnalyzer {
0050 public:
0051   typedef std::map<unsigned int, std::vector<PSimHit>> SimHitsMap;
0052   typedef std::map<unsigned int, SimTrack> SimTracksMap;
0053
0054   explicit Phase2ITValidateCluster(const edm::ParameterSet&);
0055   ~Phase2ITValidateCluster() override;
0056   void bookHistograms(DQMStore::IBooker& ibooker, edm::Run const& iRun, edm::EventSetup const& iSetup) override;
0057   void analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) override;
0058   void dqmBeginRun(const edm::Run& iRun, const edm::EventSetup& iSetup) override;
0059   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0060
0061 private:
0062   struct ClusterMEs {
0063     MonitorElement* deltaX_P = nullptr;
0064     MonitorElement* deltaY_P = nullptr;
0065     MonitorElement* deltaX_P_primary = nullptr;
0066     MonitorElement* deltaY_P_primary = nullptr;
0067   };
0068
0069   void fillITHistos(const edm::Event& iEvent,
0070                     const std::vector<edm::Handle<edm::PSimHitContainer>>& simHits,
0071                     const std::map<unsigned int, SimTrack>& simTracks);
0072   void bookLayerHistos(DQMStore::IBooker& ibooker, uint32_t det_it, const std::string& subdir);
0073   std::vector<unsigned int> getSimTrackId(const edm::Handle<edm::DetSetVector<PixelDigiSimLink>>& pixelSimLinks,
0074                                           const DetId& detId,
0075                                           unsigned int channel);
0076
0077   std::map<std::string, ClusterMEs> layerMEs_;
0078
0079   edm::ParameterSet config_;
0080   double simtrackminpt_;
0081   std::vector<edm::EDGetTokenT<edm::PSimHitContainer>> simHitTokens_;
0082   edm::EDGetTokenT<edm::DetSetVector<PixelDigiSimLink>> simITLinksToken_;
0083   edm::EDGetTokenT<edm::SimTrackContainer> simTracksToken_;
0084   edm::EDGetTokenT<edmNew::DetSetVector<SiPixelCluster>> clustersToken_;
0085   std::vector<edm::InputTag> pSimHitSrc_;
0086   const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> geomToken_;
0087   const edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> topoToken_;
0088   const TrackerGeometry* tkGeom_ = nullptr;
0089   const TrackerTopology* tTopo_ = nullptr;
0090 };
0091 #include "Validation/SiTrackerPhase2V/interface/TrackerPhase2ValidationUtil.h"
0092 #include "DQM/SiTrackerPhase2/interface/TrackerPhase2DQMUtil.h"
0093
0094 //
0095 // constructors
0096 //
0097
0098 Phase2ITValidateCluster::Phase2ITValidateCluster(const edm::ParameterSet& iConfig)
0099     : config_(iConfig),
0100       simtrackminpt_(config_.getParameter<double>("SimTrackMinPt")),
0101       simITLinksToken_(consumes<edm::DetSetVector<PixelDigiSimLink>>(
0102           config_.getParameter<edm::InputTag>("InnerTrackerDigiSimLinkSource"))),
0103       simTracksToken_(consumes<edm::SimTrackContainer>(config_.getParameter<edm::InputTag>("simtracks"))),
0104       clustersToken_(
0105           consumes<edmNew::DetSetVector<SiPixelCluster>>(config_.getParameter<edm::InputTag>("ClusterSource"))),
0106       pSimHitSrc_(config_.getParameter<std::vector<edm::InputTag>>("PSimHitSource")),
0107       geomToken_(esConsumes<TrackerGeometry, TrackerDigiGeometryRecord, edm::Transition::BeginRun>()),
0108       topoToken_(esConsumes<TrackerTopology, TrackerTopologyRcd, edm::Transition::BeginRun>()) {
0109   edm::LogInfo("Phase2ITValidateCluster") << ">>> Construct Phase2ITValidateCluster ";
0110   for (const auto& itag : pSimHitSrc_)
0111     simHitTokens_.push_back(consumes<edm::PSimHitContainer>(itag));
0112 }
0113
0114 Phase2ITValidateCluster::~Phase2ITValidateCluster() {
0115   // do anything here that needs to be done at desctruction time
0116   // (e.g. close files, deallocate resources etc.)
0117   edm::LogInfo("Phase2ITValidateCluster") << ">>> Destroy Phase2ITValidateCluster ";
0118 }
0119 //
0120 // -- DQM Begin Run
0121 //
0122 void Phase2ITValidateCluster::dqmBeginRun(const edm::Run& iRun, const edm::EventSetup& iSetup) {
0123   tkGeom_ = &iSetup.getData(geomToken_);
0124   tTopo_ = &iSetup.getData(topoToken_);
0125 }
0126
0127 // -- Analyze
0128 //
0129 void Phase2ITValidateCluster::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0130   // Getting simHits
0131   std::vector<edm::Handle<edm::PSimHitContainer>> simHits;
0132   for (const auto& itoken : simHitTokens_) {
0133     const auto& simHitHandle = iEvent.getHandle(itoken);
0134     if (!simHitHandle.isValid())
0135       continue;
0136     simHits.emplace_back(simHitHandle);
0137   }
0138   // Get the SimTracks
0139   const auto& simTracksRaw = iEvent.getHandle(simTracksToken_);
0140   // Rearrange the simTracks for ease of use <simTrackID, simTrack>
0141   SimTracksMap simTracks;
0142   for (edm::SimTrackContainer::const_iterator simTrackIt(simTracksRaw->begin()); simTrackIt != simTracksRaw->end();
0143        ++simTrackIt) {
0144     if (simTrackIt->momentum().pt() > simtrackminpt_) {
0145       simTracks.emplace(simTrackIt->trackId(), *simTrackIt);
0146     }
0147   }
0148   fillITHistos(iEvent, simHits, simTracks);
0149 }
0150
0151 void Phase2ITValidateCluster::fillITHistos(const edm::Event& iEvent,
0152                                            const std::vector<edm::Handle<edm::PSimHitContainer>>& simHits,
0153                                            const std::map<unsigned int, SimTrack>& simTracks) {
0154   // Getting the clusters
0155   const auto& clusterHandle = iEvent.getHandle(clustersToken_);
0156   if (!clusterHandle.isValid()) {
0157     edm::LogWarning("Phase2ITValidateCluster") << "No SiPixelCluster Collection found in the event. Skipping!";
0158     return;
0159   }
0160
0161   // Getting PixelDigiSimLinks
0162   const auto& pixelSimLinksHandle = iEvent.getHandle(simITLinksToken_);
0163   if (!pixelSimLinksHandle.isValid()) {
0164     edm::LogWarning("Phase2ITValidateCluster") << "No PixelDigiSimLinks Collection found in the event. Skipping!";
0165     return;
0166   }
0167
0168   for (const auto& DSVItr : *clusterHandle) {
0169     // Getting the id of detector unit
0170     uint32_t rawid = DSVItr.detId();
0171     DetId detId(rawid);
0172     const GeomDetUnit* geomDetUnit(tkGeom_->idToDetUnit(detId));
0173     if (!geomDetUnit)
0174       continue;
0175
0176     std::string folderkey = phase2tkutil::getITHistoId(detId, tTopo_);
0177     for (const auto& clusterItr : DSVItr) {
0178       MeasurementPoint mpCluster(clusterItr.x(), clusterItr.y());
0179       Local3DPoint localPosCluster = geomDetUnit->topology().localPosition(mpCluster);
0180       // Get simTracks from the cluster
0181       std::vector<unsigned int> clusterSimTrackIds;
0182       for (int irow = clusterItr.minPixelRow(); irow <= clusterItr.maxPixelRow(); ++irow) {
0183         for (int icol = clusterItr.minPixelCol(); icol <= clusterItr.maxPixelCol(); ++icol) {
0184           uint32_t channel = PixelChannelIdentifier::pixelToChannel(irow, icol);
0185           std::vector<unsigned int> simTrackIds(getSimTrackId(pixelSimLinksHandle, detId, channel));
0186           for (auto it : simTrackIds) {
0187             bool add = true;
0188             for (unsigned int j = 0; j < clusterSimTrackIds.size(); ++j) {
0189               // only save simtrackids that are not present yet
0190               if (it == clusterSimTrackIds.at(j))
0191                 add = false;
0192             }
0193             if (add)
0194               clusterSimTrackIds.push_back(it);
0195           }
0196         }
0197       }
0198       std::sort(clusterSimTrackIds.begin(), clusterSimTrackIds.end());
0199       const PSimHit* closestSimHit = nullptr;
0200       float mind = 10000.;
0201       // Get the SimHit
0202       for (const auto& psimhitCont : simHits) {
0203         for (const auto& simhitIt : *psimhitCont) {
0204           if (rawid == simhitIt.detUnitId()) {
0205             auto it = std::lower_bound(clusterSimTrackIds.begin(), clusterSimTrackIds.end(), simhitIt.trackId());
0206             if (it != clusterSimTrackIds.end() && *it == simhitIt.trackId()) {
0207               float dx = simhitIt.localPosition().x() - localPosCluster.x();
0208               float dy = simhitIt.localPosition().y() - localPosCluster.y();
0209               float dist = dx * dx + dy * dy;
0210               if (!closestSimHit || dist < mind) {
0211                 mind = dist;
0212                 closestSimHit = &simhitIt;
0213               }
0214             }
0215           }
0216         }  //end loop over PSimhitcontainers
0217       }  //end loop over simHits
0218
0219       if (!closestSimHit)
0220         continue;
0221       // only look at simhits from highpT tracks
0222       auto simTrackIt(simTracks.find(closestSimHit->trackId()));
0223       if (simTrackIt == simTracks.end())
0224         continue;
0225       Local3DPoint localPosSimHit(closestSimHit->localPosition());
0226       const double deltaX = phase2tkutil::cmtomicron * (localPosCluster.x() - localPosSimHit.x());
0227       const double deltaY = phase2tkutil::cmtomicron * (localPosCluster.y() - localPosSimHit.y());
0228
0229       auto layerMEIt = layerMEs_.find(folderkey);
0230       if (layerMEIt == layerMEs_.end())
0231         continue;
0232
0233       ClusterMEs& local_mes = layerMEIt->second;
0234       local_mes.deltaX_P->Fill(deltaX);
0235       local_mes.deltaY_P->Fill(deltaY);
0236       // Primary particles only
0237       if (phase2tkutil::isPrimary(simTrackIt->second, closestSimHit)) {
0238         local_mes.deltaX_P_primary->Fill(deltaX);
0239         local_mes.deltaY_P_primary->Fill(deltaY);
0240       }
0241     }
0242   }
0243 }
0244
0245 //
0246 // -- Book Histograms
0247 //
0248 void Phase2ITValidateCluster::bookHistograms(DQMStore::IBooker& ibooker,
0249                                              edm::Run const& iRun,
0250                                              edm::EventSetup const& iSetup) {
0251   std::string top_folder = config_.getParameter<std::string>("TopFolderName");
0252   edm::LogInfo("Phase2ITValidateCluster") << " Booking Histograms in: " << top_folder;
0253
0254   edm::ESWatcher<TrackerDigiGeometryRecord> theTkDigiGeomWatcher;
0255   if (theTkDigiGeomWatcher.check(iSetup)) {
0256     for (auto const& det_u : tkGeom_->detUnits()) {
0257       //Always check TrackerNumberingBuilder before changing this part
0258       if (!(det_u->subDetector() == GeomDetEnumerators::SubDetector::P2PXB ||
0259             det_u->subDetector() == GeomDetEnumerators::SubDetector::P2PXEC))
0260         continue;  //continue if not Pixel
0261       uint32_t detId_raw = det_u->geographicalId().rawId();
0262       bookLayerHistos(ibooker, detId_raw, top_folder);
0263     }
0264   }
0265 }
0266
0267 //////////////////Layer Histo/////////////////////////////////
0268 void Phase2ITValidateCluster::bookLayerHistos(DQMStore::IBooker& ibooker, uint32_t det_id, const std::string& subdir) {
0269   std::string folderName = phase2tkutil::getITHistoId(det_id, tTopo_);
0270   if (folderName.empty()) {
0271     edm::LogWarning("Phase2ITValidateCluster") << ">>>> Invalid histo_id ";
0272     return;
0273   }
0274
0275   if (layerMEs_.find(folderName) == layerMEs_.end()) {
0276     ibooker.cd();
0277     ibooker.setCurrentFolder(subdir + '/' + folderName);
0278     edm::LogInfo("Phase2TrackerValidateDigi") << " Booking Histograms in: " << subdir + '/' + folderName;
0279     ClusterMEs local_mes;
0280
0281     local_mes.deltaX_P =
0282         phase2tkutil::book1DFromPSet(config_.getParameter<edm::ParameterSet>("Delta_X_Pixel"), ibooker);
0283
0284     local_mes.deltaY_P =
0285         phase2tkutil::book1DFromPSet(config_.getParameter<edm::ParameterSet>("Delta_Y_Pixel"), ibooker);
0286
0287     // Puting primary digis in a subfolder
0288     ibooker.setCurrentFolder(subdir + '/' + folderName + "/PrimarySimHits");
0289
0290     local_mes.deltaX_P_primary =
0291         phase2tkutil::book1DFromPSet(config_.getParameter<edm::ParameterSet>("Delta_X_Pixel_Primary"), ibooker);
0292
0293     local_mes.deltaY_P_primary =
0294         phase2tkutil::book1DFromPSet(config_.getParameter<edm::ParameterSet>("Delta_Y_Pixel_Primary"), ibooker);
0295     layerMEs_.emplace(folderName, local_mes);
0296   }
0297 }
0298
0299 std::vector<unsigned int> Phase2ITValidateCluster::getSimTrackId(
0300     const edm::Handle<edm::DetSetVector<PixelDigiSimLink>>& pixelSimLinks, const DetId& detId, unsigned int channel) {
0301   std::vector<unsigned int> retvec;
0302   edm::DetSetVector<PixelDigiSimLink>::const_iterator DSViter(pixelSimLinks->find(detId));
0303   if (DSViter == pixelSimLinks->end())
0304     return retvec;
0305   for (edm::DetSet<PixelDigiSimLink>::const_iterator it = DSViter->data.begin(); it != DSViter->data.end(); ++it) {
0306     if (channel == it->channel()) {
0307       retvec.push_back(it->SimTrackId());
0308     }
0309   }
0310   return retvec;
0311 }
0312
0313 void Phase2ITValidateCluster::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0314   edm::ParameterSetDescription desc;
0315   //for macro-pixel sensors
0316   {
0317     edm::ParameterSetDescription psd0;
0318     psd0.add<std::string>("name", "Delta_X_Pixel");
0319     psd0.add<std::string>("title", "#Delta X;Cluster resolution X coordinate [#mum]");
0320     psd0.add<bool>("switch", true);
0321     psd0.add<double>("xmax", 250.);
0322     psd0.add<double>("xmin", -250.);
0323     psd0.add<int>("NxBins", 100);
0324     desc.add<edm::ParameterSetDescription>("Delta_X_Pixel", psd0);
0325   }
0326   {
0327     edm::ParameterSetDescription psd0;
0328     psd0.add<std::string>("name", "Delta_Y_Pixel");
0329     psd0.add<std::string>("title", "#Delta Y ;Cluster resolution Y coordinate [#mum]");
0330     psd0.add<double>("xmin", -250.);
0331     psd0.add<double>("xmax", 250.);
0332     psd0.add<bool>("switch", true);
0333     psd0.add<int>("NxBins", 100);
0334     desc.add<edm::ParameterSetDescription>("Delta_Y_Pixel", psd0);
0335   }
0336   {
0337     edm::ParameterSetDescription psd0;
0338     psd0.add<std::string>("name", "Delta_X_Pixel_Primary");
0339     psd0.add<std::string>("title", "#Delta X ;cluster resolution X coordinate [#mum]");
0340     psd0.add<double>("xmin", -250.);
0341     psd0.add<double>("xmax", 250.);
0342     psd0.add<bool>("switch", true);
0343     psd0.add<int>("NxBins", 100);
0344     desc.add<edm::ParameterSetDescription>("Delta_X_Pixel_Primary", psd0);
0345   }
0346   {
0347     edm::ParameterSetDescription psd0;
0348     psd0.add<std::string>("name", "Delta_Y_Pixel_Primary");
0349     psd0.add<std::string>("title", "#Delta Y ;cluster resolution Y coordinate [#mum]");
0350     psd0.add<double>("xmin", -250.);
0351     psd0.add<double>("xmax", 250.);
0352     psd0.add<bool>("switch", true);
0353     psd0.add<int>("NxBins", 100);
0354     desc.add<edm::ParameterSetDescription>("Delta_Y_Pixel_Primary", psd0);
0355   }
0356
0357   desc.add<std::string>("TopFolderName", "TrackerPhase2ITClusterV");
0358   desc.add<edm::InputTag>("ClusterSource", edm::InputTag("siPixelClusters"));
0359   desc.add<edm::InputTag>("InnerTrackerDigiSimLinkSource", edm::InputTag("simSiPixelDigis", "Pixel"));
0360   desc.add<edm::InputTag>("simtracks", edm::InputTag("g4SimHits"));
0361   desc.add<double>("SimTrackMinPt", 0.0);
0362   desc.add<std::vector<edm::InputTag>>("PSimHitSource",
0363                                        {
0364                                            edm::InputTag("g4SimHits:TrackerHitsPixelBarrelLowTof"),
0365                                            edm::InputTag("g4SimHits:TrackerHitsPixelBarrelHighTof"),
0366                                            edm::InputTag("g4SimHits:TrackerHitsPixelEndcapLowTof"),
0367                                            edm::InputTag("g4SimHits:TrackerHitsPixelEndcapHighTof"),
0368                                        });
0369   descriptions.add("Phase2ITValidateCluster", desc);
0370 }
0371 DEFINE_FWK_MODULE(Phase2ITValidateCluster);