Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-05-26 22:39:26

0001 // -*- C++ -*-
0002 // Package:    SiPixelPhase1CompareTrackSoA
0003 // Class:      SiPixelPhase1CompareTrackSoA
0004 //
0005 /**\class SiPixelPhase1CompareTrackSoA SiPixelPhase1CompareTrackSoA.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/PixelTrackHeterogeneous.h"
0024 // for string manipulations
0025 #include <fmt/printf.h>
0026 
0027 namespace {
0028   // same logic used for the MTV:
0029   // cf https://github.com/cms-sw/cmssw/blob/master/Validation/RecoTrack/src/MTVHistoProducerAlgoForTracker.cc
0030   typedef dqm::reco::DQMStore DQMStore;
0031 
0032   void setBinLog(TAxis* axis) {
0033     int bins = axis->GetNbins();
0034     float from = axis->GetXmin();
0035     float to = axis->GetXmax();
0036     float width = (to - from) / bins;
0037     std::vector<float> new_bins(bins + 1, 0);
0038     for (int i = 0; i <= bins; i++) {
0039       new_bins[i] = TMath::Power(10, from + i * width);
0040     }
0041     axis->Set(bins, new_bins.data());
0042   }
0043 
0044   void setBinLogX(TH1* h) {
0045     TAxis* axis = h->GetXaxis();
0046     setBinLog(axis);
0047   }
0048   void setBinLogY(TH1* h) {
0049     TAxis* axis = h->GetYaxis();
0050     setBinLog(axis);
0051   }
0052 
0053   template <typename... Args>
0054   dqm::reco::MonitorElement* make2DIfLog(DQMStore::IBooker& ibook, bool logx, bool logy, Args&&... args) {
0055     auto h = std::make_unique<TH2I>(std::forward<Args>(args)...);
0056     if (logx)
0057       setBinLogX(h.get());
0058     if (logy)
0059       setBinLogY(h.get());
0060     const auto& name = h->GetName();
0061     return ibook.book2I(name, h.release());
0062   }
0063 }  // namespace
0064 
0065 class SiPixelPhase1CompareTrackSoA : public DQMEDAnalyzer {
0066 public:
0067   explicit SiPixelPhase1CompareTrackSoA(const edm::ParameterSet&);
0068   ~SiPixelPhase1CompareTrackSoA() override = default;
0069   void bookHistograms(DQMStore::IBooker& ibooker, edm::Run const& iRun, edm::EventSetup const& iSetup) override;
0070   void analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) override;
0071   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0072 
0073 private:
0074   const edm::EDGetTokenT<PixelTrackHeterogeneous> tokenSoATrackCPU_;
0075   const edm::EDGetTokenT<PixelTrackHeterogeneous> tokenSoATrackGPU_;
0076   const std::string topFolderName_;
0077   const bool useQualityCut_;
0078   const pixelTrack::Quality minQuality_;
0079   const float dr2cut_;
0080   MonitorElement* hnTracks_;
0081   MonitorElement* hnLooseAndAboveTracks_;
0082   MonitorElement* hnLooseAndAboveTracks_matched_;
0083   MonitorElement* hnHits_;
0084   MonitorElement* hnHitsVsPhi_;
0085   MonitorElement* hnHitsVsEta_;
0086   MonitorElement* hnLayers_;
0087   MonitorElement* hnLayersVsPhi_;
0088   MonitorElement* hnLayersVsEta_;
0089   MonitorElement* hchi2_;
0090   MonitorElement* hChi2VsPhi_;
0091   MonitorElement* hChi2VsEta_;
0092   MonitorElement* hpt_;
0093   MonitorElement* hptLogLog_;
0094   MonitorElement* heta_;
0095   MonitorElement* hphi_;
0096   MonitorElement* hz_;
0097   MonitorElement* htip_;
0098   MonitorElement* hquality_;
0099   //1D differences
0100   MonitorElement* hptdiffMatched_;
0101   MonitorElement* hetadiffMatched_;
0102   MonitorElement* hphidiffMatched_;
0103   MonitorElement* hzdiffMatched_;
0104   //for matching eff vs region:dervie ration at harvesting
0105   MonitorElement* hpt_eta_tkAllCPU_;
0106   MonitorElement* hpt_eta_tkAllCPUMatched_;
0107   MonitorElement* hphi_z_tkAllCPU_;
0108   MonitorElement* hphi_z_tkAllCPUMatched_;
0109 };
0110 
0111 //
0112 // constructors
0113 //
0114 
0115 SiPixelPhase1CompareTrackSoA::SiPixelPhase1CompareTrackSoA(const edm::ParameterSet& iConfig)
0116     : tokenSoATrackCPU_(consumes<PixelTrackHeterogeneous>(iConfig.getParameter<edm::InputTag>("pixelTrackSrcCPU"))),
0117       tokenSoATrackGPU_(consumes<PixelTrackHeterogeneous>(iConfig.getParameter<edm::InputTag>("pixelTrackSrcGPU"))),
0118       topFolderName_(iConfig.getParameter<std::string>("topFolderName")),
0119       useQualityCut_(iConfig.getParameter<bool>("useQualityCut")),
0120       minQuality_(pixelTrack::qualityByName(iConfig.getParameter<std::string>("minQuality"))),
0121       dr2cut_(iConfig.getParameter<double>("deltaR2cut")) {}
0122 
0123 //
0124 // -- Analyze
0125 //
0126 void SiPixelPhase1CompareTrackSoA::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0127   const auto& tsoaHandleCPU = iEvent.getHandle(tokenSoATrackCPU_);
0128   const auto& tsoaHandleGPU = iEvent.getHandle(tokenSoATrackGPU_);
0129   if (not tsoaHandleCPU or not tsoaHandleGPU) {
0130     edm::LogWarning out("SiPixelPhase1CompareTrackSoA");
0131     if (not tsoaHandleCPU) {
0132       out << "reference (cpu) tracks not found; ";
0133     }
0134     if (not tsoaHandleGPU) {
0135       out << "target (gpu) tracks not found; ";
0136     }
0137     out << "the comparison will not run.";
0138     return;
0139   }
0140 
0141   auto const& tsoaCPU = *tsoaHandleCPU->get();
0142   auto const& tsoaGPU = *tsoaHandleGPU->get();
0143   auto maxTracksCPU = tsoaCPU.stride();  //this should be same for both?
0144   auto maxTracksGPU = tsoaGPU.stride();  //this should be same for both?
0145   auto const* qualityCPU = tsoaCPU.qualityData();
0146   auto const* qualityGPU = tsoaGPU.qualityData();
0147   int32_t nTracksCPU = 0;
0148   int32_t nTracksGPU = 0;
0149   int32_t nLooseAndAboveTracksCPU = 0;
0150   int32_t nLooseAndAboveTracksCPU_matchedGPU = 0;
0151   int32_t nLooseAndAboveTracksGPU = 0;
0152 
0153   //Loop over GPU tracks and store the indices of the loose tracks. Whats happens if useQualityCut_ is false?
0154   std::vector<int32_t> looseTrkidxGPU;
0155   for (int32_t jt = 0; jt < maxTracksGPU; ++jt) {
0156     if (tsoaGPU.nHits(jt) == 0)
0157       break;  // this is a guard
0158     if (!(tsoaGPU.pt(jt) > 0.))
0159       continue;
0160     nTracksGPU++;
0161     if (useQualityCut_ && qualityGPU[jt] < minQuality_)
0162       continue;
0163     nLooseAndAboveTracksGPU++;
0164     looseTrkidxGPU.emplace_back(jt);
0165   }
0166 
0167   //Now loop over CPU tracks//nested loop for loose gPU tracks
0168   for (int32_t it = 0; it < maxTracksCPU; ++it) {
0169     if (tsoaCPU.nHits(it) == 0)
0170       break;  // this is a guard
0171     if (!(tsoaCPU.pt(it) > 0.))
0172       continue;
0173     nTracksCPU++;
0174     if (useQualityCut_ && qualityCPU[it] < minQuality_)
0175       continue;
0176     nLooseAndAboveTracksCPU++;
0177     //Now loop over loose GPU trk and find the closest in DeltaR//do we need pt cut?
0178     const int32_t notFound = -1;
0179     int32_t closestTkidx = notFound;
0180     float mindr2 = dr2cut_;
0181     float etacpu = tsoaCPU.eta(it);
0182     float phicpu = tsoaCPU.phi(it);
0183     for (auto gid : looseTrkidxGPU) {
0184       float etagpu = tsoaGPU.eta(gid);
0185       float phigpu = tsoaGPU.phi(gid);
0186       float dr2 = reco::deltaR2(etacpu, phicpu, etagpu, phigpu);
0187       if (dr2 > dr2cut_)
0188         continue;  // this is arbitrary
0189       if (mindr2 > dr2) {
0190         mindr2 = dr2;
0191         closestTkidx = gid;
0192       }
0193     }
0194 
0195     hpt_eta_tkAllCPU_->Fill(etacpu, tsoaCPU.pt(it));  //all CPU tk
0196     hphi_z_tkAllCPU_->Fill(phicpu, tsoaCPU.zip(it));
0197     if (closestTkidx == notFound)
0198       continue;
0199     nLooseAndAboveTracksCPU_matchedGPU++;
0200 
0201     hchi2_->Fill(tsoaCPU.chi2(it), tsoaGPU.chi2(closestTkidx));
0202     hnHits_->Fill(tsoaCPU.nHits(it), tsoaGPU.nHits(closestTkidx));
0203     hnLayers_->Fill(tsoaCPU.nLayers(it), tsoaGPU.nLayers(closestTkidx));
0204     hpt_->Fill(tsoaCPU.pt(it), tsoaGPU.pt(closestTkidx));
0205     hptLogLog_->Fill(tsoaCPU.pt(it), tsoaGPU.pt(closestTkidx));
0206     heta_->Fill(etacpu, tsoaGPU.eta(closestTkidx));
0207     hphi_->Fill(phicpu, tsoaGPU.phi(closestTkidx));
0208     hz_->Fill(tsoaCPU.zip(it), tsoaGPU.zip(closestTkidx));
0209     htip_->Fill(tsoaCPU.tip(it), tsoaGPU.tip(closestTkidx));
0210     hptdiffMatched_->Fill(tsoaCPU.pt(it) - tsoaGPU.pt(closestTkidx));
0211     hetadiffMatched_->Fill(etacpu - tsoaGPU.eta(closestTkidx));
0212     hphidiffMatched_->Fill(reco::deltaPhi(phicpu, tsoaGPU.phi(closestTkidx)));
0213     hzdiffMatched_->Fill(tsoaCPU.zip(it) - tsoaGPU.zip(closestTkidx));
0214     hpt_eta_tkAllCPUMatched_->Fill(etacpu, tsoaCPU.pt(it));  //matched to gpu
0215     hphi_z_tkAllCPUMatched_->Fill(phicpu, tsoaCPU.zip(it));
0216   }
0217   hnTracks_->Fill(nTracksCPU, nTracksGPU);
0218   hnLooseAndAboveTracks_->Fill(nLooseAndAboveTracksCPU, nLooseAndAboveTracksGPU);
0219   hnLooseAndAboveTracks_matched_->Fill(nLooseAndAboveTracksCPU, nLooseAndAboveTracksCPU_matchedGPU);
0220 }
0221 
0222 //
0223 // -- Book Histograms
0224 //
0225 void SiPixelPhase1CompareTrackSoA::bookHistograms(DQMStore::IBooker& iBook,
0226                                                   edm::Run const& iRun,
0227                                                   edm::EventSetup const& iSetup) {
0228   iBook.cd();
0229   iBook.setCurrentFolder(topFolderName_);
0230 
0231   // clang-format off
0232   std::string toRep = "Number of tracks";
0233   // FIXME: all the 2D correlation plots are quite heavy in terms of memory consumption, so a as soon as DQM supports either TH2I or THnSparse
0234   // these should be moved to a less resource consuming format
0235   hnTracks_ = iBook.book2I("nTracks", fmt::sprintf("%s per event; CPU; GPU",toRep), 501, -0.5, 500.5, 501, -0.5, 500.5);
0236   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);
0237   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);
0238 
0239   toRep = "Number of all RecHits per track (quality #geq loose)";
0240   hnHits_ = iBook.book2I("nRecHits", fmt::sprintf("%s;CPU;GPU",toRep), 15, -0.5, 14.5, 15, -0.5, 14.5);
0241 
0242   toRep = "Number of all layers per track (quality #geq loose)";
0243   hnLayers_ = iBook.book2I("nLayers", fmt::sprintf("%s;CPU;GPU",toRep), 15, -0.5, 14.5, 15, -0.5, 14.5);
0244 
0245   toRep = "Track (quality #geq loose) #chi^{2}/ndof";
0246   hchi2_ = iBook.book2I("nChi2ndof", fmt::sprintf("%s;CPU;GPU",toRep), 40, 0., 20., 40, 0., 20.);
0247 
0248   hpt_ = iBook.book2I("pt", "Track (quality #geq loose) p_{T} [GeV];CPU;GPU", 200, 0., 200., 200, 0., 200.);
0249   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.));
0250   heta_ = iBook.book2I("eta", "Track (quality #geq loose) #eta;CPU;GPU", 30, -3., 3., 30, -3., 3.);
0251   hphi_ = iBook.book2I("phi", "Track (quality #geq loose) #phi;CPU;GPU", 30, -M_PI, M_PI, 30, -M_PI, M_PI);
0252   hz_ = iBook.book2I("z", "Track (quality #geq loose) z [cm];CPU;GPU", 30, -30., 30., 30, -30., 30.);
0253   htip_ = iBook.book2I("tip", "Track (quality #geq loose) TIP [cm];CPU;GPU", 100, -0.5, 0.5, 100, -0.5, 0.5);
0254   //1D difference plots
0255   hptdiffMatched_ = iBook.book1D("ptdiffmatched", " p_{T} diff [GeV] between matched tracks; #Delta p_{T} [GeV]", 60, -30., 30.);
0256   hetadiffMatched_ = iBook.book1D("etadiffmatched", " #eta diff between matched tracks; #Delta #eta", 160, -0.04 ,0.04);
0257   hphidiffMatched_ = iBook.book1D("phidiffmatched", " #phi diff between matched tracks; #Delta #phi",  160, -0.04 ,0.04);
0258   hzdiffMatched_ = iBook.book1D("zdiffmatched", " z diff between matched tracks; #Delta z [cm]", 300, -1.5, 1.5);
0259   //2D plots for eff
0260   hpt_eta_tkAllCPU_ = iBook.book2I("ptetatrkAllCPU", "Track (quality #geq loose) on CPU; #eta; p_{T} [GeV];", 30, -M_PI, M_PI, 200, 0., 200.);
0261   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.);
0262 
0263   hphi_z_tkAllCPU_ = iBook.book2I("phiztrkAllCPU", "Track (quality #geq loose) on CPU; #phi; z [cm];",  30, -M_PI, M_PI, 30, -30., 30.);
0264   hphi_z_tkAllCPUMatched_ = iBook.book2I("phiztrkAllCPUmatched", "Track (quality #geq loose) on CPU; #phi; z [cm];", 30, -M_PI, M_PI, 30, -30., 30.);
0265 
0266 }
0267 
0268 void SiPixelPhase1CompareTrackSoA::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0269   // monitorpixelTrackSoA
0270   edm::ParameterSetDescription desc;
0271   desc.add<edm::InputTag>("pixelTrackSrcCPU", edm::InputTag("pixelTracksSoA@cpu"));
0272   desc.add<edm::InputTag>("pixelTrackSrcGPU", edm::InputTag("pixelTracksSoA@cuda"));
0273   desc.add<std::string>("topFolderName", "SiPixelHeterogeneous/PixelTrackCompareGPUvsCPU");
0274   desc.add<bool>("useQualityCut", true);
0275   desc.add<std::string>("minQuality", "loose");
0276   desc.add<double>("deltaR2cut", 0.04);
0277   descriptions.addWithDefaultLabel(desc);
0278 }
0279 DEFINE_FWK_MODULE(SiPixelPhase1CompareTrackSoA);