File indexing completed on 2024-08-30 02:10:34
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
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
0026 #include "DQMServices/Core/interface/MonitorElement.h"
0027 #include "DQMServices/Core/interface/DQMEDAnalyzer.h"
0028 #include "DQMServices/Core/interface/DQMStore.h"
0029
0030 #include "DataFormats/TrackSoA/interface/TracksHost.h"
0031 #include "DataFormats/TrackSoA/interface/alpaka/TrackUtilities.h"
0032
0033 namespace {
0034
0035
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 }
0070
0071
0072 template <typename T>
0073 class SiPixelCompareTracks : public DQMEDAnalyzer {
0074 public:
0075 using PixelTrackSoA = TracksHost<T>;
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
0082
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
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
0120 MonitorElement* hptdiffMatched_;
0121 MonitorElement* hCurvdiffMatched_;
0122 MonitorElement* hetadiffMatched_;
0123 MonitorElement* hphidiffMatched_;
0124 MonitorElement* hzdiffMatched_;
0125 MonitorElement* htipdiffMatched_;
0126
0127
0128 MonitorElement* hpt_eta_tkAllRef_;
0129 MonitorElement* hpt_eta_tkAllRefMatched_;
0130 MonitorElement* hphi_z_tkAllRef_;
0131 MonitorElement* hphi_z_tkAllRefMatched_;
0132 };
0133
0134
0135
0136
0137
0138 template <typename T>
0139 SiPixelCompareTracks<T>::SiPixelCompareTracks(const edm::ParameterSet& iConfig)
0140 : tokenSoATrackReference_(consumes<PixelTrackSoA>(iConfig.getParameter<edm::InputTag>("pixelTrackReferenceSoA"))),
0141 tokenSoATrackTarget_(consumes<PixelTrackSoA>(iConfig.getParameter<edm::InputTag>("pixelTrackTargetSoA"))),
0142 topFolderName_(iConfig.getParameter<std::string>("topFolderName")),
0143 useQualityCut_(iConfig.getParameter<bool>("useQualityCut")),
0144 minQuality_(pixelTrack::qualityByName(iConfig.getParameter<std::string>("minQuality"))),
0145 dr2cut_(iConfig.getParameter<double>("deltaR2cut")) {}
0146
0147 template <typename T>
0148 template <typename U, typename V>
0149 void SiPixelCompareTracks<T>::analyzeSeparate(U tokenRef, V tokenTar, const edm::Event& iEvent) {
0150 using helper = TracksUtilities<T>;
0151
0152 const auto& tsoaHandleRef = iEvent.getHandle(tokenRef);
0153 const auto& tsoaHandleTar = iEvent.getHandle(tokenTar);
0154
0155 if (not tsoaHandleRef or not tsoaHandleTar) {
0156 edm::LogWarning out("SiPixelCompareTracks");
0157 if (not tsoaHandleRef) {
0158 out << "reference tracks not found; ";
0159 }
0160 if (not tsoaHandleTar) {
0161 out << "target tracks not found; ";
0162 }
0163 out << "the comparison will not run.";
0164 return;
0165 }
0166
0167 auto const& tsoaRef = *tsoaHandleRef;
0168 auto const& tsoaTar = *tsoaHandleTar;
0169
0170 auto maxTracksRef = tsoaRef.view().metadata().size();
0171 auto maxTracksTar = tsoaTar.view().metadata().size();
0172
0173 auto const* qualityRef = tsoaRef.view().quality();
0174 auto const* qualityTar = tsoaTar.view().quality();
0175
0176 int32_t nTracksRef = 0;
0177 int32_t nTracksTar = 0;
0178 int32_t nLooseAndAboveTracksRef = 0;
0179 int32_t nLooseAndAboveTracksRef_matchedTar = 0;
0180 int32_t nLooseAndAboveTracksTar = 0;
0181
0182
0183 std::vector<int32_t> looseTrkidxTar;
0184 for (int32_t jt = 0; jt < maxTracksTar; ++jt) {
0185 if (helper::nHits(tsoaTar.view(), jt) == 0)
0186 break;
0187 if (!(tsoaTar.view()[jt].pt() > 0.))
0188 continue;
0189 nTracksTar++;
0190 if (useQualityCut_ && qualityTar[jt] < minQuality_)
0191 continue;
0192 nLooseAndAboveTracksTar++;
0193 looseTrkidxTar.emplace_back(jt);
0194 }
0195
0196
0197 for (int32_t it = 0; it < maxTracksRef; ++it) {
0198 int nHitsRef = helper::nHits(tsoaRef.view(), it);
0199
0200 if (nHitsRef == 0)
0201 break;
0202
0203 float ptRef = tsoaRef.view()[it].pt();
0204 float etaRef = tsoaRef.view()[it].eta();
0205 float phiRef = reco::phi(tsoaRef.view(), it);
0206 float zipRef = reco::zip(tsoaRef.view(), it);
0207 float tipRef = reco::tip(tsoaRef.view(), it);
0208 auto qRef = reco::charge(tsoaRef.view(), it);
0209
0210 if (!(ptRef > 0.))
0211 continue;
0212 nTracksRef++;
0213 if (useQualityCut_ && qualityRef[it] < minQuality_)
0214 continue;
0215 nLooseAndAboveTracksRef++;
0216
0217 const int32_t notFound = -1;
0218 int32_t closestTkidx = notFound;
0219 float mindr2 = dr2cut_;
0220
0221 for (auto gid : looseTrkidxTar) {
0222 float etaTar = tsoaTar.view()[gid].eta();
0223 float phiTar = reco::phi(tsoaTar.view(), gid);
0224 float dr2 = reco::deltaR2(etaRef, phiRef, etaTar, phiTar);
0225 if (dr2 > dr2cut_)
0226 continue;
0227 if (mindr2 > dr2) {
0228 mindr2 = dr2;
0229 closestTkidx = gid;
0230 }
0231 }
0232
0233 hpt_eta_tkAllRef_->Fill(etaRef, ptRef);
0234 hphi_z_tkAllRef_->Fill(phiRef, zipRef);
0235 if (closestTkidx == notFound)
0236 continue;
0237 nLooseAndAboveTracksRef_matchedTar++;
0238
0239 hchi2_->Fill(tsoaRef.view()[it].chi2(), tsoaTar.view()[closestTkidx].chi2());
0240 hCharge_->Fill(qRef, reco::charge(tsoaTar.view(), closestTkidx));
0241 hnHits_->Fill(helper::nHits(tsoaRef.view(), it), helper::nHits(tsoaTar.view(), closestTkidx));
0242 hnLayers_->Fill(tsoaRef.view()[it].nLayers(), tsoaTar.view()[closestTkidx].nLayers());
0243 hpt_->Fill(ptRef, tsoaTar.view()[closestTkidx].pt());
0244 hCurvature_->Fill(qRef / ptRef, reco::charge(tsoaTar.view(), closestTkidx) / tsoaTar.view()[closestTkidx].pt());
0245 hptLogLog_->Fill(ptRef, tsoaTar.view()[closestTkidx].pt());
0246 heta_->Fill(etaRef, tsoaTar.view()[closestTkidx].eta());
0247 hphi_->Fill(phiRef, reco::phi(tsoaTar.view(), closestTkidx));
0248 hz_->Fill(zipRef, reco::zip(tsoaTar.view(), closestTkidx));
0249 htip_->Fill(tipRef, reco::tip(tsoaTar.view(), closestTkidx));
0250 hptdiffMatched_->Fill(ptRef - tsoaTar.view()[closestTkidx].pt());
0251 hCurvdiffMatched_->Fill(qRef / ptRef -
0252 (reco::charge(tsoaTar.view(), closestTkidx) / tsoaTar.view()[closestTkidx].pt()));
0253 hetadiffMatched_->Fill(etaRef - tsoaTar.view()[closestTkidx].eta());
0254 hphidiffMatched_->Fill(reco::deltaPhi(phiRef, reco::phi(tsoaTar.view(), closestTkidx)));
0255 hzdiffMatched_->Fill(zipRef - reco::zip(tsoaTar.view(), closestTkidx));
0256 htipdiffMatched_->Fill(tipRef - reco::tip(tsoaTar.view(), closestTkidx));
0257 hpt_eta_tkAllRefMatched_->Fill(etaRef, tsoaRef.view()[it].pt());
0258 hphi_z_tkAllRefMatched_->Fill(etaRef, zipRef);
0259 }
0260
0261
0262 auto fillHistogram = [](auto& histogram, auto xValue, auto yValue) { histogram->Fill(xValue, yValue); };
0263
0264
0265 auto fillDeltaHistogram = [](auto& histogram, int cpuValue, int gpuValue) {
0266 histogram->Fill(std::min(cpuValue, 1000), std::clamp(gpuValue - cpuValue, -100, 100));
0267 };
0268
0269
0270 fillHistogram(hnTracks_, nTracksRef, nTracksTar);
0271 fillHistogram(hnLooseAndAboveTracks_, nLooseAndAboveTracksRef, nLooseAndAboveTracksTar);
0272 fillHistogram(hnLooseAndAboveTracks_matched_, nLooseAndAboveTracksRef, nLooseAndAboveTracksRef_matchedTar);
0273
0274 fillDeltaHistogram(hDeltaNTracks_, nTracksRef, nTracksTar);
0275 fillDeltaHistogram(hDeltaNLooseAndAboveTracks_, nLooseAndAboveTracksRef, nLooseAndAboveTracksTar);
0276 fillDeltaHistogram(hDeltaNLooseAndAboveTracks_matched_, nLooseAndAboveTracksRef, nLooseAndAboveTracksRef_matchedTar);
0277 }
0278
0279
0280
0281
0282 template <typename T>
0283 void SiPixelCompareTracks<T>::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0284
0285
0286 analyzeSeparate(tokenSoATrackReference_, tokenSoATrackTarget_, iEvent);
0287 }
0288
0289
0290
0291
0292 template <typename T>
0293 void SiPixelCompareTracks<T>::bookHistograms(DQMStore::IBooker& iBook,
0294 edm::Run const& iRun,
0295 edm::EventSetup const& iSetup) {
0296 iBook.cd();
0297 iBook.setCurrentFolder(topFolderName_);
0298
0299
0300 std::string toRep = "Number of tracks";
0301 auto bookTracksTH2I = [&](const std::string& name,
0302 const std::string& title,
0303 int xBins,
0304 double xMin,
0305 double xMax,
0306 int yBins,
0307 double yMin,
0308 double yMax) {
0309 return iBook.book2I(name, fmt::sprintf(title, toRep), xBins, xMin, xMax, yBins, yMin, yMax);
0310 };
0311
0312
0313 constexpr int xBins = 501;
0314 constexpr double xMin = -0.5;
0315 constexpr double xMax = 1001.5;
0316
0317 constexpr int dXBins = 1001;
0318 constexpr double dXMin = -0.5;
0319 constexpr double dXMax = 1000.5;
0320
0321 constexpr int dYBins = 201;
0322 constexpr double dYMin = -100.5;
0323 constexpr double dYMax = 100.5;
0324
0325
0326
0327
0328
0329
0330 hnTracks_ = bookTracksTH2I("nTracks", "%s per event; Reference; Target", xBins, xMin, xMax, xBins, xMin, xMax);
0331 hnLooseAndAboveTracks_ = bookTracksTH2I("nLooseAndAboveTracks", "%s (quality #geq loose) per event; Reference; Target", xBins, xMin, xMax, xBins, xMin, xMax);
0332 hnLooseAndAboveTracks_matched_ = bookTracksTH2I("nLooseAndAboveTracks_matched", "%s (quality #geq loose) per event; Reference; Target", xBins, xMin, xMax, xBins, xMin, xMax);
0333
0334 hDeltaNTracks_ = bookTracksTH2I("deltaNTracks", "%s per event; Reference; Target - Reference", dXBins, dXMin, dXMax, dYBins, dYMin, dYMax);
0335 hDeltaNLooseAndAboveTracks_ = bookTracksTH2I("deltaNLooseAndAboveTracks", "%s (quality #geq loose) per event; Reference; Target - Reference", dXBins, dXMin, dXMax, dYBins, dYMin, dYMax);
0336 hDeltaNLooseAndAboveTracks_matched_ = bookTracksTH2I("deltaNLooseAndAboveTracks_matched", "%s (quality #geq loose) per event; Reference; Target - Reference", dXBins, dXMin, dXMax, dYBins, dYMin, dYMax);
0337
0338 toRep = "Number of all RecHits per track (quality #geq loose)";
0339 hnHits_ = iBook.book2I("nRecHits", fmt::sprintf("%s;Reference;Target",toRep), 15, -0.5, 14.5, 15, -0.5, 14.5);
0340
0341 toRep = "Number of all layers per track (quality #geq loose)";
0342 hnLayers_ = iBook.book2I("nLayers", fmt::sprintf("%s;Reference;Target",toRep), 15, -0.5, 14.5, 15, -0.5, 14.5);
0343
0344 toRep = "Track (quality #geq loose) #chi^{2}/ndof";
0345 hchi2_ = iBook.book2I("nChi2ndof", fmt::sprintf("%s;Reference;Target",toRep), 40, 0., 20., 40, 0., 20.);
0346
0347 toRep = "Track (quality #geq loose) charge";
0348 hCharge_ = iBook.book2I("charge",fmt::sprintf("%s;Reference;Target",toRep),3, -1.5, 1.5, 3, -1.5, 1.5);
0349
0350 hpt_ = iBook.book2I("pt", "Track (quality #geq loose) p_{T} [GeV];Reference;Target", 200, 0., 200., 200, 0., 200.);
0351 hCurvature_ = iBook.book2I("curvature", "Track (quality #geq loose) q/p_{T} [GeV^{-1}];Reference;Target", 60,- 3., 3., 60, -3., 3. );
0352 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.));
0353 heta_ = iBook.book2I("eta", "Track (quality #geq loose) #eta;Reference;Target", 30, -3., 3., 30, -3., 3.);
0354 hphi_ = iBook.book2I("phi", "Track (quality #geq loose) #phi;Reference;Target", 30, -M_PI, M_PI, 30, -M_PI, M_PI);
0355 hz_ = iBook.book2I("z", "Track (quality #geq loose) z [cm];Reference;Target", 30, -30., 30., 30, -30., 30.);
0356 htip_ = iBook.book2I("tip", "Track (quality #geq loose) TIP [cm];Reference;Target", 100, -0.5, 0.5, 100, -0.5, 0.5);
0357
0358
0359 hptdiffMatched_ = iBook.book1D("ptdiffmatched", " p_{T} diff [GeV] between matched tracks; #Delta p_{T} [GeV]", 61, -30.5, 30.5);
0360 hCurvdiffMatched_ = iBook.book1D("curvdiffmatched", "q/p_{T} diff [GeV^{-1}] between matched tracks; #Delta q/p_{T} [GeV^{-1}]", 61, -3.05, 3.05);
0361 hetadiffMatched_ = iBook.book1D("etadiffmatched", " #eta diff between matched tracks; #Delta #eta", 161, -0.045 ,0.045);
0362 hphidiffMatched_ = iBook.book1D("phidiffmatched", " #phi diff between matched tracks; #Delta #phi", 161, -0.045 ,0.045);
0363 hzdiffMatched_ = iBook.book1D("zdiffmatched", " z diff between matched tracks; #Delta z [cm]", 301, -1.55, 1.55);
0364 htipdiffMatched_ = iBook.book1D("tipdiffmatched", " TIP diff between matched tracks; #Delta TIP [cm]", 301, -1.55, 1.55);
0365
0366 hpt_eta_tkAllRef_ = iBook.book2I("ptetatrkAllReference", "Track (quality #geq loose) on Reference; #eta; p_{T} [GeV];", 30, -M_PI, M_PI, 200, 0., 200.);
0367 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.);
0368
0369 hphi_z_tkAllRef_ = iBook.book2I("phiztrkAllReference", "Track (quality #geq loose) on Reference; #phi; z [cm];", 30, -M_PI, M_PI, 30, -30., 30.);
0370 hphi_z_tkAllRefMatched_ = iBook.book2I("phiztrkAllReferencematched", "Track (quality #geq loose) on Reference; #phi; z [cm];", 30, -M_PI, M_PI, 30, -30., 30.);
0371
0372 }
0373
0374 template<typename T>
0375 void SiPixelCompareTracks<T>::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0376
0377 edm::ParameterSetDescription desc;
0378 desc.add<edm::InputTag>("pixelTrackReferenceSoA", edm::InputTag("pixelTracksAlpakaSerial"));
0379 desc.add<edm::InputTag>("pixelTrackTargetSoA", edm::InputTag("pixelTracksAlpaka"));
0380 desc.add<std::string>("topFolderName", "SiPixelHeterogeneous/PixelTrackCompareDeviceVSHost");
0381 desc.add<bool>("useQualityCut", true);
0382 desc.add<std::string>("minQuality", "loose");
0383 desc.add<double>("deltaR2cut", 0.02 * 0.02)->setComment("deltaR2 cut between track on device and host");
0384 descriptions.addWithDefaultLabel(desc);
0385 }
0386
0387
0388
0389 using SiPixelPhase1CompareTracks = SiPixelCompareTracks<pixelTopology::Phase1>;
0390 using SiPixelPhase2CompareTracks = SiPixelCompareTracks<pixelTopology::Phase2>;
0391 using SiPixelHIonPhase1CompareTracks = SiPixelCompareTracks<pixelTopology::HIonPhase1>;
0392
0393 DEFINE_FWK_MODULE(SiPixelPhase1CompareTracks);
0394 DEFINE_FWK_MODULE(SiPixelPhase2CompareTracks);
0395 DEFINE_FWK_MODULE(SiPixelHIonPhase1CompareTracks);
0396