Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-08-30 02:10:35 

0001 // -*- C++ -*-
0002 // Package:    SiPixelCompareTrackSoA
0003 // Class:      SiPixelCompareTrackSoA
0004 //
0005 /**\class SiPixelCompareTrackSoA SiPixelCompareTrackSoA.cc
0006 */
0007 //
0008 // Author: Suvankar Roy Chowdhury
0009 //
0010 #include "DataFormats/Common/interface/Handle.h"
0011 #include "DataFormats/Math/interface/deltaR.h"
0012 #include "DataFormats/Math/interface/deltaPhi.h"
0013 #include "FWCore/Framework/interface/Event.h"
0014 #include "FWCore/Framework/interface/Frameworkfwd.h"
0015 #include "FWCore/Framework/interface/MakerMacros.h"
0016 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0017 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0018 #include "FWCore/Utilities/interface/InputTag.h"
0019 // DQM Histograming
0020 #include "DQMServices/Core/interface/MonitorElement.h"
0021 #include "DQMServices/Core/interface/DQMEDAnalyzer.h"
0022 #include "DQMServices/Core/interface/DQMStore.h"
0023 #include "CUDADataFormats/Track/interface/TrackSoAHeterogeneousHost.h"
0024 #include "CUDADataFormats/Track/interface/TrackSoAHeterogeneousDevice.h"
0025 // for string manipulations
0026 #include <fmt/printf.h>
0027 
0028 namespace {
0029   // same logic used for the MTV:
0030   // cf https://github.com/cms-sw/cmssw/blob/master/Validation/RecoTrack/src/MTVHistoProducerAlgoForTracker.cc
0031   typedef dqm::reco::DQMStore DQMStore;
0032 
0033   void setBinLog(TAxis* axis) {
0034     int bins = axis->GetNbins();
0035     float from = axis->GetXmin();
0036     float to = axis->GetXmax();
0037     float width = (to - from) / bins;
0038     std::vector<float> new_bins(bins + 1, 0);
0039     for (int i = 0; i <= bins; i++) {
0040       new_bins[i] = TMath::Power(10, from + i * width);
0041     }
0042     axis->Set(bins, new_bins.data());
0043   }
0044 
0045   void setBinLogX(TH1* h) {
0046     TAxis* axis = h->GetXaxis();
0047     setBinLog(axis);
0048   }
0049   void setBinLogY(TH1* h) {
0050     TAxis* axis = h->GetYaxis();
0051     setBinLog(axis);
0052   }
0053 
0054   template <typename... Args>
0055   dqm::reco::MonitorElement* make2DIfLog(DQMStore::IBooker& ibook, bool logx, bool logy, Args&&... args) {
0056     auto h = std::make_unique<TH2I>(std::forward<Args>(args)...);
0057     if (logx)
0058       setBinLogX(h.get());
0059     if (logy)
0060       setBinLogY(h.get());
0061     const auto& name = h->GetName();
0062     return ibook.book2I(name, h.release());
0063   }
0064 }  // namespace
0065 
0066 template <typename T>
0067 class SiPixelCompareTrackSoA : public DQMEDAnalyzer {
0068 public:
0069   using PixelTrackSoA = TrackSoAHeterogeneousHost<T>;
0070 
0071   explicit SiPixelCompareTrackSoA(const edm::ParameterSet&);
0072   ~SiPixelCompareTrackSoA() override = default;
0073   void bookHistograms(DQMStore::IBooker& ibooker, edm::Run const& iRun, edm::EventSetup const& iSetup) override;
0074   void analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) override;
0075   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0076 
0077 private:
0078   const edm::EDGetTokenT<PixelTrackSoA> tokenSoATrackCPU_;
0079   const edm::EDGetTokenT<PixelTrackSoA> tokenSoATrackGPU_;
0080   const std::string topFolderName_;
0081   const bool useQualityCut_;
0082   const pixelTrack::Quality minQuality_;
0083   const float dr2cut_;
0084   MonitorElement* hnTracks_;
0085   MonitorElement* hnLooseAndAboveTracks_;
0086   MonitorElement* hnLooseAndAboveTracks_matched_;
0087   MonitorElement* hDeltaNTracks_;
0088   MonitorElement* hDeltaNLooseAndAboveTracks_;
0089   MonitorElement* hDeltaNLooseAndAboveTracks_matched_;
0090   MonitorElement* hnHits_;
0091   MonitorElement* hnHitsVsPhi_;
0092   MonitorElement* hnHitsVsEta_;
0093   MonitorElement* hnLayers_;
0094   MonitorElement* hnLayersVsPhi_;
0095   MonitorElement* hnLayersVsEta_;
0096   MonitorElement* hCharge_;
0097   MonitorElement* hchi2_;
0098   MonitorElement* hChi2VsPhi_;
0099   MonitorElement* hChi2VsEta_;
0100   MonitorElement* hpt_;
0101   MonitorElement* hCurvature_;
0102   MonitorElement* hptLogLog_;
0103   MonitorElement* heta_;
0104   MonitorElement* hphi_;
0105   MonitorElement* hz_;
0106   MonitorElement* htip_;
0107   MonitorElement* hquality_;
0108   //1D differences
0109   MonitorElement* hptdiffMatched_;
0110   MonitorElement* hCurvdiffMatched_;
0111   MonitorElement* hetadiffMatched_;
0112   MonitorElement* hphidiffMatched_;
0113   MonitorElement* hzdiffMatched_;
0114   MonitorElement* htipdiffMatched_;
0115 
0116   //for matching eff vs region: derive the ratio at harvesting
0117   MonitorElement* hpt_eta_tkAllRef_;
0118   MonitorElement* hpt_eta_tkAllRefMatched_;
0119   MonitorElement* hphi_z_tkAllRef_;
0120   MonitorElement* hphi_z_tkAllRefMatched_;
0121 };
0122 
0123 //
0124 // constructors
0125 //
0126 
0127 template <typename T>
0128 SiPixelCompareTrackSoA<T>::SiPixelCompareTrackSoA(const edm::ParameterSet& iConfig)
0129     : tokenSoATrackCPU_(consumes<PixelTrackSoA>(iConfig.getParameter<edm::InputTag>("pixelTrackSrcCPU"))),
0130       tokenSoATrackGPU_(consumes<PixelTrackSoA>(iConfig.getParameter<edm::InputTag>("pixelTrackSrcGPU"))),
0131       topFolderName_(iConfig.getParameter<std::string>("topFolderName")),
0132       useQualityCut_(iConfig.getParameter<bool>("useQualityCut")),
0133       minQuality_(pixelTrack::qualityByName(iConfig.getParameter<std::string>("minQuality"))),
0134       dr2cut_(iConfig.getParameter<double>("deltaR2cut")) {}
0135 
0136 //
0137 // -- Analyze
0138 //
0139 template <typename T>
0140 void SiPixelCompareTrackSoA<T>::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0141   using helper = TracksUtilities<T>;
0142   const auto& tsoaHandleCPU = iEvent.getHandle(tokenSoATrackCPU_);
0143   const auto& tsoaHandleGPU = iEvent.getHandle(tokenSoATrackGPU_);
0144   if (not tsoaHandleCPU or not tsoaHandleGPU) {
0145     edm::LogWarning out("SiPixelCompareTrackSoA");
0146     if (not tsoaHandleCPU) {
0147       out << "reference (cpu) tracks not found; ";
0148     }
0149     if (not tsoaHandleGPU) {
0150       out << "target (gpu) tracks not found; ";
0151     }
0152     out << "the comparison will not run.";
0153     return;
0154   }
0155 
0156   auto const& tsoaCPU = *tsoaHandleCPU;
0157   auto const& tsoaGPU = *tsoaHandleGPU;
0158   auto maxTracksCPU = tsoaCPU.view().metadata().size();  //this should be same for both?
0159   auto maxTracksGPU = tsoaGPU.view().metadata().size();  //this should be same for both?
0160   auto const* qualityCPU = tsoaCPU.view().quality();
0161   auto const* qualityGPU = tsoaGPU.view().quality();
0162   int32_t nTracksCPU = 0;
0163   int32_t nTracksGPU = 0;
0164   int32_t nLooseAndAboveTracksCPU = 0;
0165   int32_t nLooseAndAboveTracksCPU_matchedGPU = 0;
0166   int32_t nLooseAndAboveTracksGPU = 0;
0167 
0168   //Loop over GPU tracks and store the indices of the loose tracks. Whats happens if useQualityCut_ is false?
0169   std::vector<int32_t> looseTrkidxGPU;
0170   for (int32_t jt = 0; jt < maxTracksGPU; ++jt) {
0171     if (helper::nHits(tsoaGPU.view(), jt) == 0)
0172       break;  // this is a guard
0173     if (!(tsoaGPU.view()[jt].pt() > 0.))
0174       continue;
0175     nTracksGPU++;
0176     if (useQualityCut_ && qualityGPU[jt] < minQuality_)
0177       continue;
0178     nLooseAndAboveTracksGPU++;
0179     looseTrkidxGPU.emplace_back(jt);
0180   }
0181 
0182   //Now loop over CPU tracks//nested loop for loose gPU tracks
0183   for (int32_t it = 0; it < maxTracksCPU; ++it) {
0184     int nHitsCPU = helper::nHits(tsoaCPU.view(), it);
0185 
0186     if (nHitsCPU == 0)
0187       break;  // this is a guard
0188 
0189     float ptCPU = tsoaCPU.view()[it].pt();
0190     float etaCPU = tsoaCPU.view()[it].eta();
0191     float phiCPU = helper::phi(tsoaCPU.view(), it);
0192     float zipCPU = helper::zip(tsoaCPU.view(), it);
0193     float tipCPU = helper::tip(tsoaCPU.view(), it);
0194     auto qCPU = helper::charge(tsoaCPU.view(), it);
0195 
0196     if (!(ptCPU > 0.))
0197       continue;
0198     nTracksCPU++;
0199     if (useQualityCut_ && qualityCPU[it] < minQuality_)
0200       continue;
0201     nLooseAndAboveTracksCPU++;
0202     //Now loop over loose GPU trk and find the closest in DeltaR//do we need pt cut?
0203     const int32_t notFound = -1;
0204     int32_t closestTkidx = notFound;
0205     float mindr2 = dr2cut_;
0206 
0207     for (auto gid : looseTrkidxGPU) {
0208       float etaGPU = tsoaGPU.view()[gid].eta();
0209       float phiGPU = helper::phi(tsoaGPU.view(), gid);
0210       float dr2 = reco::deltaR2(etaCPU, phiCPU, etaGPU, phiGPU);
0211       if (dr2 > dr2cut_)
0212         continue;  // this is arbitrary
0213       if (mindr2 > dr2) {
0214         mindr2 = dr2;
0215         closestTkidx = gid;
0216       }
0217     }
0218 
0219     hpt_eta_tkAllRef_->Fill(etaCPU, ptCPU);  //all CPU tk
0220     hphi_z_tkAllRef_->Fill(phiCPU, zipCPU);
0221     if (closestTkidx == notFound)
0222       continue;
0223     nLooseAndAboveTracksCPU_matchedGPU++;
0224 
0225     hchi2_->Fill(tsoaCPU.view()[it].chi2(), tsoaGPU.view()[closestTkidx].chi2());
0226     hCharge_->Fill(qCPU, helper::charge(tsoaGPU.view(), closestTkidx));
0227     hnHits_->Fill(helper::nHits(tsoaCPU.view(), it), helper::nHits(tsoaGPU.view(), closestTkidx));
0228     hnLayers_->Fill(tsoaCPU.view()[it].nLayers(), tsoaGPU.view()[closestTkidx].nLayers());
0229     hpt_->Fill(ptCPU, tsoaGPU.view()[closestTkidx].pt());
0230     hCurvature_->Fill(qCPU / ptCPU, helper::charge(tsoaGPU.view(), closestTkidx) / tsoaGPU.view()[closestTkidx].pt());
0231     hptLogLog_->Fill(ptCPU, tsoaGPU.view()[closestTkidx].pt());
0232     heta_->Fill(etaCPU, tsoaGPU.view()[closestTkidx].eta());
0233     hphi_->Fill(phiCPU, helper::phi(tsoaGPU.view(), closestTkidx));
0234     hz_->Fill(zipCPU, helper::zip(tsoaGPU.view(), closestTkidx));
0235     htip_->Fill(tipCPU, helper::tip(tsoaGPU.view(), closestTkidx));
0236     hptdiffMatched_->Fill(ptCPU - tsoaGPU.view()[closestTkidx].pt());
0237     hCurvdiffMatched_->Fill((helper::charge(tsoaCPU.view(), it) / tsoaCPU.view()[it].pt()) -
0238                             (helper::charge(tsoaGPU.view(), closestTkidx) / tsoaGPU.view()[closestTkidx].pt()));
0239     hetadiffMatched_->Fill(etaCPU - tsoaGPU.view()[closestTkidx].eta());
0240     hphidiffMatched_->Fill(reco::deltaPhi(phiCPU, helper::phi(tsoaGPU.view(), closestTkidx)));
0241     hzdiffMatched_->Fill(zipCPU - helper::zip(tsoaGPU.view(), closestTkidx));
0242     htipdiffMatched_->Fill(tipCPU - helper::tip(tsoaGPU.view(), closestTkidx));
0243     hpt_eta_tkAllRefMatched_->Fill(etaCPU, tsoaCPU.view()[it].pt());  //matched to gpu
0244     hphi_z_tkAllRefMatched_->Fill(etaCPU, zipCPU);
0245   }
0246 
0247   // Define a lambda function for filling the histograms
0248   auto fillHistogram = [](auto& histogram, auto xValue, auto yValue) { histogram->Fill(xValue, yValue); };
0249 
0250   // Define a lambda for filling delta histograms
0251   auto fillDeltaHistogram = [](auto& histogram, int cpuValue, int gpuValue) {
0252     histogram->Fill(std::min(cpuValue, 1000), std::clamp(gpuValue - cpuValue, -100, 100));
0253   };
0254 
0255   // Fill the histograms
0256   fillHistogram(hnTracks_, nTracksCPU, nTracksGPU);
0257   fillHistogram(hnLooseAndAboveTracks_, nLooseAndAboveTracksCPU, nLooseAndAboveTracksGPU);
0258   fillHistogram(hnLooseAndAboveTracks_matched_, nLooseAndAboveTracksCPU, nLooseAndAboveTracksCPU_matchedGPU);
0259 
0260   fillDeltaHistogram(hDeltaNTracks_, nTracksCPU, nTracksGPU);
0261   fillDeltaHistogram(hDeltaNLooseAndAboveTracks_, nLooseAndAboveTracksCPU, nLooseAndAboveTracksGPU);
0262   fillDeltaHistogram(hDeltaNLooseAndAboveTracks_matched_, nLooseAndAboveTracksCPU, nLooseAndAboveTracksCPU_matchedGPU);
0263 }
0264 
0265 //
0266 // -- Book Histograms
0267 //
0268 template <typename T>
0269 void SiPixelCompareTrackSoA<T>::bookHistograms(DQMStore::IBooker& iBook,
0270                                                edm::Run const& iRun,
0271                                                edm::EventSetup const& iSetup) {
0272   iBook.cd();
0273   iBook.setCurrentFolder(topFolderName_);
0274 
0275   // Define a helper function for booking histograms
0276   std::string toRep = "Number of tracks";
0277   auto bookTracksTH2I = [&](const std::string& name,
0278                             const std::string& title,
0279                             int xBins,
0280                             double xMin,
0281                             double xMax,
0282                             int yBins,
0283                             double yMin,
0284                             double yMax) {
0285     return iBook.book2I(name, fmt::sprintf(title, toRep), xBins, xMin, xMax, yBins, yMin, yMax);
0286   };
0287 
0288   // Define common parameters for different histogram types
0289   constexpr int xBins = 501;
0290   constexpr double xMin = -0.5;
0291   constexpr double xMax = 1001.5;
0292 
0293   constexpr int dXBins = 1001;
0294   constexpr double dXMin = -0.5;
0295   constexpr double dXMax = 1000.5;
0296 
0297   constexpr int dYBins = 201;
0298   constexpr double dYMin = -100.5;
0299   constexpr double dYMax = 100.5;
0300 
0301   // FIXME: all the 2D correlation plots are quite heavy in terms of memory consumption, so a as soon as DQM supports THnSparse
0302   // these should be moved to a less resource consuming format
0303 
0304   // Book histograms using the helper function
0305   // clang-format off
0306   hnTracks_ = bookTracksTH2I("nTracks", "%s per event; Reference; Target", xBins, xMin, xMax, xBins, xMin, xMax);
0307   hnLooseAndAboveTracks_ = bookTracksTH2I("nLooseAndAboveTracks", "%s (quality #geq loose) per event; Reference; Target", xBins, xMin, xMax, xBins, xMin, xMax);
0308   hnLooseAndAboveTracks_matched_ = bookTracksTH2I("nLooseAndAboveTracks_matched", "%s (quality #geq loose) per event; Reference; Target", xBins, xMin, xMax, xBins, xMin, xMax);
0309 
0310   hDeltaNTracks_ = bookTracksTH2I("deltaNTracks", "%s per event; Reference; Target - Reference", dXBins, dXMin, dXMax, dYBins, dYMin, dYMax);
0311   hDeltaNLooseAndAboveTracks_ = bookTracksTH2I("deltaNLooseAndAboveTracks", "%s (quality #geq loose) per event; Reference; Target - Reference", dXBins, dXMin, dXMax, dYBins, dYMin, dYMax);
0312   hDeltaNLooseAndAboveTracks_matched_ = bookTracksTH2I("deltaNLooseAndAboveTracks_matched", "%s (quality #geq loose) per event; Reference; Target - Reference", dXBins, dXMin, dXMax, dYBins, dYMin, dYMax);
0313 
0314   toRep = "Number of all RecHits per track (quality #geq loose)";
0315   hnHits_ = iBook.book2I("nRecHits", fmt::sprintf("%s;CPU;GPU",toRep), 15, -0.5, 14.5, 15, -0.5, 14.5);
0316 
0317   toRep = "Number of all layers per track (quality #geq loose)";
0318   hnLayers_ = iBook.book2I("nLayers", fmt::sprintf("%s;CPU;GPU",toRep), 15, -0.5, 14.5, 15, -0.5, 14.5);
0319 
0320   toRep = "Track (quality #geq loose) #chi^{2}/ndof";
0321   hchi2_ = iBook.book2I("nChi2ndof", fmt::sprintf("%s;CPU;GPU",toRep), 40, 0., 20., 40, 0., 20.);
0322 
0323   toRep = "Track (quality #geq loose) charge";
0324   hCharge_ = iBook.book2I("charge",fmt::sprintf("%s;CPU;GPU",toRep),3, -1.5, 1.5, 3, -1.5, 1.5);
0325 
0326   hpt_ = iBook.book2I("pt", "Track (quality #geq loose) p_{T} [GeV];CPU;GPU", 200, 0., 200., 200, 0., 200.);
0327   hCurvature_ = iBook.book2I("curvature", "Track (quality #geq loose) q/p_{T} [GeV^{-1}];CPU;GPU",  60,- 3., 3., 60, -3., 3. );
0328   hptLogLog_ = make2DIfLog(iBook, true, true, "ptLogLog", "Track (quality #geq loose) p_{T} [GeV];CPU;GPU", 200, log10(0.5), log10(200.), 200, log10(0.5), log10(200.));
0329   heta_ = iBook.book2I("eta", "Track (quality #geq loose) #eta;CPU;GPU", 30, -3., 3., 30, -3., 3.);
0330   hphi_ = iBook.book2I("phi", "Track (quality #geq loose) #phi;CPU;GPU", 30, -M_PI, M_PI, 30, -M_PI, M_PI);
0331   hz_ = iBook.book2I("z", "Track (quality #geq loose) z [cm];CPU;GPU", 30, -30., 30., 30, -30., 30.);
0332   htip_ = iBook.book2I("tip", "Track (quality #geq loose) TIP [cm];CPU;GPU", 100, -0.5, 0.5, 100, -0.5, 0.5);
0333   //1D difference plots
0334   hptdiffMatched_ = iBook.book1D("ptdiffmatched", " p_{T} diff [GeV] between matched tracks; #Delta p_{T} [GeV]", 61, -30.5, 30.5);
0335   hCurvdiffMatched_ = iBook.book1D("curvdiffmatched", "q/p_{T} diff [GeV^{-1}] between matched tracks; #Delta q/p_{T} [GeV^{-1}]", 61, -3.05, 3.05);
0336   hetadiffMatched_ = iBook.book1D("etadiffmatched", " #eta diff between matched tracks; #Delta #eta", 161, -0.045 ,0.045);
0337   hphidiffMatched_ = iBook.book1D("phidiffmatched", " #phi diff between matched tracks; #Delta #phi",  161, -0.045 ,0.045);
0338   hzdiffMatched_ = iBook.book1D("zdiffmatched", " z diff between matched tracks; #Delta z [cm]", 301, -1.55, 1.55);
0339   htipdiffMatched_ = iBook.book1D("tipdiffmatched", " TIP diff between matched tracks; #Delta TIP [cm]", 301, -1.55, 1.55);
0340   //2D plots for eff
0341   hpt_eta_tkAllRef_ = iBook.book2I("ptetatrkAllReference", "Track (quality #geq loose) on CPU; #eta; p_{T} [GeV];", 30, -M_PI, M_PI, 200, 0., 200.);
0342   hpt_eta_tkAllRefMatched_ = iBook.book2I("ptetatrkAllReferencematched", "Track (quality #geq loose) on CPU matched to GPU track; #eta; p_{T} [GeV];", 30, -M_PI, M_PI, 200, 0., 200.);
0343 
0344   hphi_z_tkAllRef_ = iBook.book2I("phiztrkAllReference", "Track (quality #geq loose) on CPU; #phi; z [cm];",  30, -M_PI, M_PI, 30, -30., 30.);
0345   hphi_z_tkAllRefMatched_ = iBook.book2I("phiztrkAllReferencematched", "Track (quality #geq loose) on CPU; #phi; z [cm];", 30, -M_PI, M_PI, 30, -30., 30.);
0346 
0347 }
0348 
0349 template<typename T>
0350 void SiPixelCompareTrackSoA<T>::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0351   // monitorpixelTrackSoA
0352   edm::ParameterSetDescription desc;
0353   desc.add<edm::InputTag>("pixelTrackSrcCPU", edm::InputTag("pixelTracksSoA@cpu"));
0354   desc.add<edm::InputTag>("pixelTrackSrcGPU", edm::InputTag("pixelTracksSoA@cuda"));
0355   desc.add<std::string>("topFolderName", "SiPixelHeterogeneous/PixelTrackCompareGPUvsCPU");
0356   desc.add<bool>("useQualityCut", true);
0357   desc.add<std::string>("minQuality", "loose");
0358   desc.add<double>("deltaR2cut", 0.02 * 0.02)->setComment("deltaR2 cut between track on CPU and GPU");
0359   descriptions.addWithDefaultLabel(desc);
0360 }
0361 
0362 using SiPixelPhase1CompareTrackSoA = SiPixelCompareTrackSoA<pixelTopology::Phase1>;
0363 using SiPixelPhase2CompareTrackSoA = SiPixelCompareTrackSoA<pixelTopology::Phase2>;
0364 using SiPixelHIonPhase1CompareTrackSoA = SiPixelCompareTrackSoA<pixelTopology::HIonPhase1>;
0365 
0366 DEFINE_FWK_MODULE(SiPixelPhase1CompareTrackSoA);
0367 DEFINE_FWK_MODULE(SiPixelPhase2CompareTrackSoA);
0368 DEFINE_FWK_MODULE(SiPixelHIonPhase1CompareTrackSoA);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!