Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-03-17 10:56:00 

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* hnHits_;
0088   MonitorElement* hnHitsVsPhi_;
0089   MonitorElement* hnHitsVsEta_;
0090   MonitorElement* hnLayers_;
0091   MonitorElement* hnLayersVsPhi_;
0092   MonitorElement* hnLayersVsEta_;
0093   MonitorElement* hCharge_;
0094   MonitorElement* hchi2_;
0095   MonitorElement* hChi2VsPhi_;
0096   MonitorElement* hChi2VsEta_;
0097   MonitorElement* hpt_;
0098   MonitorElement* hptLogLog_;
0099   MonitorElement* heta_;
0100   MonitorElement* hphi_;
0101   MonitorElement* hz_;
0102   MonitorElement* htip_;
0103   MonitorElement* hquality_;
0104   //1D differences
0105   MonitorElement* hptdiffMatched_;
0106   MonitorElement* hCurvdiffMatched_;
0107   MonitorElement* hetadiffMatched_;
0108   MonitorElement* hphidiffMatched_;
0109   MonitorElement* hzdiffMatched_;
0110   MonitorElement* htipdiffMatched_;
0111 
0112   //for matching eff vs region: derive the ratio at harvesting
0113   MonitorElement* hpt_eta_tkAllCPU_;
0114   MonitorElement* hpt_eta_tkAllCPUMatched_;
0115   MonitorElement* hphi_z_tkAllCPU_;
0116   MonitorElement* hphi_z_tkAllCPUMatched_;
0117 };
0118 
0119 //
0120 // constructors
0121 //
0122 
0123 template <typename T>
0124 SiPixelCompareTrackSoA<T>::SiPixelCompareTrackSoA(const edm::ParameterSet& iConfig)
0125     : tokenSoATrackCPU_(consumes<PixelTrackSoA>(iConfig.getParameter<edm::InputTag>("pixelTrackSrcCPU"))),
0126       tokenSoATrackGPU_(consumes<PixelTrackSoA>(iConfig.getParameter<edm::InputTag>("pixelTrackSrcGPU"))),
0127       topFolderName_(iConfig.getParameter<std::string>("topFolderName")),
0128       useQualityCut_(iConfig.getParameter<bool>("useQualityCut")),
0129       minQuality_(pixelTrack::qualityByName(iConfig.getParameter<std::string>("minQuality"))),
0130       dr2cut_(iConfig.getParameter<double>("deltaR2cut")) {}
0131 
0132 //
0133 // -- Analyze
0134 //
0135 template <typename T>
0136 void SiPixelCompareTrackSoA<T>::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0137   using helper = TracksUtilities<T>;
0138   const auto& tsoaHandleCPU = iEvent.getHandle(tokenSoATrackCPU_);
0139   const auto& tsoaHandleGPU = iEvent.getHandle(tokenSoATrackGPU_);
0140   if (not tsoaHandleCPU or not tsoaHandleGPU) {
0141     edm::LogWarning out("SiPixelCompareTrackSoA");
0142     if (not tsoaHandleCPU) {
0143       out << "reference (cpu) tracks not found; ";
0144     }
0145     if (not tsoaHandleGPU) {
0146       out << "target (gpu) tracks not found; ";
0147     }
0148     out << "the comparison will not run.";
0149     return;
0150   }
0151 
0152   auto const& tsoaCPU = *tsoaHandleCPU;
0153   auto const& tsoaGPU = *tsoaHandleGPU;
0154   auto maxTracksCPU = tsoaCPU.view().metadata().size();  //this should be same for both?
0155   auto maxTracksGPU = tsoaGPU.view().metadata().size();  //this should be same for both?
0156   auto const* qualityCPU = tsoaCPU.view().quality();
0157   auto const* qualityGPU = tsoaGPU.view().quality();
0158   int32_t nTracksCPU = 0;
0159   int32_t nTracksGPU = 0;
0160   int32_t nLooseAndAboveTracksCPU = 0;
0161   int32_t nLooseAndAboveTracksCPU_matchedGPU = 0;
0162   int32_t nLooseAndAboveTracksGPU = 0;
0163 
0164   //Loop over GPU tracks and store the indices of the loose tracks. Whats happens if useQualityCut_ is false?
0165   std::vector<int32_t> looseTrkidxGPU;
0166   for (int32_t jt = 0; jt < maxTracksGPU; ++jt) {
0167     if (helper::nHits(tsoaGPU.view(), jt) == 0)
0168       break;  // this is a guard
0169     if (!(tsoaGPU.view()[jt].pt() > 0.))
0170       continue;
0171     nTracksGPU++;
0172     if (useQualityCut_ && qualityGPU[jt] < minQuality_)
0173       continue;
0174     nLooseAndAboveTracksGPU++;
0175     looseTrkidxGPU.emplace_back(jt);
0176   }
0177 
0178   //Now loop over CPU tracks//nested loop for loose gPU tracks
0179   for (int32_t it = 0; it < maxTracksCPU; ++it) {
0180     int nHitsCPU = helper::nHits(tsoaCPU.view(), it);
0181 
0182     if (nHitsCPU == 0)
0183       break;  // this is a guard
0184 
0185     float ptCPU = tsoaCPU.view()[it].pt();
0186     float etaCPU = tsoaCPU.view()[it].eta();
0187     float phiCPU = helper::phi(tsoaCPU.view(), it);
0188     float zipCPU = helper::zip(tsoaCPU.view(), it);
0189     float tipCPU = helper::tip(tsoaCPU.view(), it);
0190 
0191     if (!(ptCPU > 0.))
0192       continue;
0193     nTracksCPU++;
0194     if (useQualityCut_ && qualityCPU[it] < minQuality_)
0195       continue;
0196     nLooseAndAboveTracksCPU++;
0197     //Now loop over loose GPU trk and find the closest in DeltaR//do we need pt cut?
0198     const int32_t notFound = -1;
0199     int32_t closestTkidx = notFound;
0200     float mindr2 = dr2cut_;
0201 
0202     for (auto gid : looseTrkidxGPU) {
0203       float etaGPU = tsoaGPU.view()[gid].eta();
0204       float phiGPU = helper::phi(tsoaGPU.view(), gid);
0205       float dr2 = reco::deltaR2(etaCPU, phiCPU, etaGPU, phiGPU);
0206       if (dr2 > dr2cut_)
0207         continue;  // this is arbitrary
0208       if (mindr2 > dr2) {
0209         mindr2 = dr2;
0210         closestTkidx = gid;
0211       }
0212     }
0213 
0214     hpt_eta_tkAllCPU_->Fill(etaCPU, ptCPU);  //all CPU tk
0215     hphi_z_tkAllCPU_->Fill(phiCPU, zipCPU);
0216     if (closestTkidx == notFound)
0217       continue;
0218     nLooseAndAboveTracksCPU_matchedGPU++;
0219 
0220     hchi2_->Fill(tsoaCPU.view()[it].chi2(), tsoaGPU.view()[closestTkidx].chi2());
0221     hCharge_->Fill(helper::charge(tsoaCPU.view(), it), helper::charge(tsoaGPU.view(), closestTkidx));
0222     hnHits_->Fill(helper::nHits(tsoaCPU.view(), it), helper::nHits(tsoaGPU.view(), closestTkidx));
0223     hnLayers_->Fill(tsoaCPU.view()[it].nLayers(), tsoaGPU.view()[closestTkidx].nLayers());
0224     hpt_->Fill(tsoaCPU.view()[it].pt(), tsoaGPU.view()[closestTkidx].pt());
0225     hptLogLog_->Fill(tsoaCPU.view()[it].pt(), tsoaGPU.view()[closestTkidx].pt());
0226     heta_->Fill(etaCPU, tsoaGPU.view()[closestTkidx].eta());
0227     hphi_->Fill(etaCPU, helper::phi(tsoaGPU.view(), closestTkidx));
0228     hz_->Fill(zipCPU, helper::zip(tsoaGPU.view(), closestTkidx));
0229     htip_->Fill(tipCPU, helper::tip(tsoaGPU.view(), closestTkidx));
0230     hptdiffMatched_->Fill(tsoaCPU.view()[it].pt() - tsoaGPU.view()[closestTkidx].pt());
0231     hCurvdiffMatched_->Fill((helper::charge(tsoaCPU.view(), it) / tsoaCPU.view()[it].pt()) -
0232                             (helper::charge(tsoaGPU.view(), closestTkidx) / tsoaGPU.view()[closestTkidx].pt()));
0233     hetadiffMatched_->Fill(etaCPU - tsoaGPU.view()[closestTkidx].eta());
0234     hphidiffMatched_->Fill(reco::deltaPhi(etaCPU, helper::phi(tsoaGPU.view(), closestTkidx)));
0235     hzdiffMatched_->Fill(zipCPU - helper::zip(tsoaGPU.view(), closestTkidx));
0236     htipdiffMatched_->Fill(tipCPU - helper::tip(tsoaGPU.view(), closestTkidx));
0237     hpt_eta_tkAllCPUMatched_->Fill(etaCPU, tsoaCPU.view()[it].pt());  //matched to gpu
0238     hphi_z_tkAllCPUMatched_->Fill(etaCPU, zipCPU);
0239   }
0240   hnTracks_->Fill(nTracksCPU, nTracksGPU);
0241   hnLooseAndAboveTracks_->Fill(nLooseAndAboveTracksCPU, nLooseAndAboveTracksGPU);
0242   hnLooseAndAboveTracks_matched_->Fill(nLooseAndAboveTracksCPU, nLooseAndAboveTracksCPU_matchedGPU);
0243 }
0244 
0245 //
0246 // -- Book Histograms
0247 //
0248 template <typename T>
0249 void SiPixelCompareTrackSoA<T>::bookHistograms(DQMStore::IBooker& iBook,
0250                                                edm::Run const& iRun,
0251                                                edm::EventSetup const& iSetup) {
0252   iBook.cd();
0253   iBook.setCurrentFolder(topFolderName_);
0254 
0255   // clang-format off
0256   std::string toRep = "Number of tracks";
0257   // FIXME: all the 2D correlation plots are quite heavy in terms of memory consumption, so a as soon as DQM supports THnSparse
0258   // these should be moved to a less resource consuming format
0259   hnTracks_ = iBook.book2I("nTracks", fmt::sprintf("%s per event; CPU; GPU",toRep), 501, -0.5, 500.5, 501, -0.5, 500.5);
0260   hnLooseAndAboveTracks_ = iBook.book2I("nLooseAndAboveTracks", fmt::sprintf("%s (quality #geq loose) per event; CPU; GPU",toRep), 501, -0.5, 500.5, 501, -0.5, 500.5);
0261   hnLooseAndAboveTracks_matched_ = iBook.book2I("nLooseAndAboveTracks_matched", fmt::sprintf("%s (quality #geq loose) per event; CPU; GPU",toRep), 501, -0.5, 500.5, 501, -0.5, 500.5);
0262 
0263   toRep = "Number of all RecHits per track (quality #geq loose)";
0264   hnHits_ = iBook.book2I("nRecHits", fmt::sprintf("%s;CPU;GPU",toRep), 15, -0.5, 14.5, 15, -0.5, 14.5);
0265 
0266   toRep = "Number of all layers per track (quality #geq loose)";
0267   hnLayers_ = iBook.book2I("nLayers", fmt::sprintf("%s;CPU;GPU",toRep), 15, -0.5, 14.5, 15, -0.5, 14.5);
0268 
0269   toRep = "Track (quality #geq loose) #chi^{2}/ndof";
0270   hchi2_ = iBook.book2I("nChi2ndof", fmt::sprintf("%s;CPU;GPU",toRep), 40, 0., 20., 40, 0., 20.);
0271 
0272   toRep = "Track (quality #geq loose) charge";
0273   hCharge_ = iBook.book2I("charge",fmt::sprintf("%s;CPU;GPU",toRep),3, -1.5, 1.5, 3, -1.5, 1.5);
0274 
0275   hpt_ = iBook.book2I("pt", "Track (quality #geq loose) p_{T} [GeV];CPU;GPU", 200, 0., 200., 200, 0., 200.);
0276   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.));
0277   heta_ = iBook.book2I("eta", "Track (quality #geq loose) #eta;CPU;GPU", 30, -3., 3., 30, -3., 3.);
0278   hphi_ = iBook.book2I("phi", "Track (quality #geq loose) #phi;CPU;GPU", 30, -M_PI, M_PI, 30, -M_PI, M_PI);
0279   hz_ = iBook.book2I("z", "Track (quality #geq loose) z [cm];CPU;GPU", 30, -30., 30., 30, -30., 30.);
0280   htip_ = iBook.book2I("tip", "Track (quality #geq loose) TIP [cm];CPU;GPU", 100, -0.5, 0.5, 100, -0.5, 0.5);
0281   //1D difference plots
0282   hptdiffMatched_ = iBook.book1D("ptdiffmatched", " p_{T} diff [GeV] between matched tracks; #Delta p_{T} [GeV]", 60, -30., 30.);
0283   hCurvdiffMatched_ = iBook.book1D("curvdiffmatched", "q/p_{T} diff [GeV] between matched tracks; #Delta q/p_{T} [GeV]", 60, -30., 30.);
0284   hetadiffMatched_ = iBook.book1D("etadiffmatched", " #eta diff between matched tracks; #Delta #eta", 160, -0.04 ,0.04);
0285   hphidiffMatched_ = iBook.book1D("phidiffmatched", " #phi diff between matched tracks; #Delta #phi",  160, -0.04 ,0.04);
0286   hzdiffMatched_ = iBook.book1D("zdiffmatched", " z diff between matched tracks; #Delta z [cm]", 300, -1.5, 1.5);
0287   htipdiffMatched_ = iBook.book1D("tipdiffmatched", " TIP diff between matched tracks; #Delta TIP [cm]", 300, -1.5, 1.5);
0288   //2D plots for eff
0289   hpt_eta_tkAllCPU_ = iBook.book2I("ptetatrkAllCPU", "Track (quality #geq loose) on CPU; #eta; p_{T} [GeV];", 30, -M_PI, M_PI, 200, 0., 200.);
0290   hpt_eta_tkAllCPUMatched_ = iBook.book2I("ptetatrkAllCPUmatched", "Track (quality #geq loose) on CPU matched to GPU track; #eta; p_{T} [GeV];", 30, -M_PI, M_PI, 200, 0., 200.);
0291 
0292   hphi_z_tkAllCPU_ = iBook.book2I("phiztrkAllCPU", "Track (quality #geq loose) on CPU; #phi; z [cm];",  30, -M_PI, M_PI, 30, -30., 30.);
0293   hphi_z_tkAllCPUMatched_ = iBook.book2I("phiztrkAllCPUmatched", "Track (quality #geq loose) on CPU; #phi; z [cm];", 30, -M_PI, M_PI, 30, -30., 30.);
0294 
0295 }
0296 
0297 template<typename T>
0298 void SiPixelCompareTrackSoA<T>::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0299   // monitorpixelTrackSoA
0300   edm::ParameterSetDescription desc;
0301   desc.add<edm::InputTag>("pixelTrackSrcCPU", edm::InputTag("pixelTracksSoA@cpu"));
0302   desc.add<edm::InputTag>("pixelTrackSrcGPU", edm::InputTag("pixelTracksSoA@cuda"));
0303   desc.add<std::string>("topFolderName", "SiPixelHeterogeneous/PixelTrackCompareGPUvsCPU");
0304   desc.add<bool>("useQualityCut", true);
0305   desc.add<std::string>("minQuality", "loose");
0306   desc.add<double>("deltaR2cut", 0.04);
0307   descriptions.addWithDefaultLabel(desc);
0308 }
0309 
0310 using SiPixelPhase1CompareTrackSoA = SiPixelCompareTrackSoA<pixelTopology::Phase1>;
0311 using SiPixelPhase2CompareTrackSoA = SiPixelCompareTrackSoA<pixelTopology::Phase2>;
0312 
0313 DEFINE_FWK_MODULE(SiPixelPhase1CompareTrackSoA);
0314 DEFINE_FWK_MODULE(SiPixelPhase2CompareTrackSoA);
[ 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!