Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-07-03 00:42:14 

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