CMS 3D CMS Logo

SiPixelCompareTrackSoAAlpaka.cc
Go to the documentation of this file.
1 // for string manipulations
2 #include <fmt/printf.h>
12 // DQM Histograming
16 // DataFormats
19 
20 namespace {
21  // same logic used for the MTV:
22  // cf https://github.com/cms-sw/cmssw/blob/master/Validation/RecoTrack/src/MTVHistoProducerAlgoForTracker.cc
24 
25  void setBinLog(TAxis* axis) {
26  int bins = axis->GetNbins();
27  float from = axis->GetXmin();
28  float to = axis->GetXmax();
29  float width = (to - from) / bins;
30  std::vector<float> new_bins(bins + 1, 0);
31  for (int i = 0; i <= bins; i++) {
32  new_bins[i] = TMath::Power(10, from + i * width);
33  }
34  axis->Set(bins, new_bins.data());
35  }
36 
37  void setBinLogX(TH1* h) {
38  TAxis* axis = h->GetXaxis();
39  setBinLog(axis);
40  }
41  void setBinLogY(TH1* h) {
42  TAxis* axis = h->GetYaxis();
43  setBinLog(axis);
44  }
45 
46  template <typename... Args>
47  dqm::reco::MonitorElement* make2DIfLog(DQMStore::IBooker& ibook, bool logx, bool logy, Args&&... args) {
48  auto h = std::make_unique<TH2I>(std::forward<Args>(args)...);
49  if (logx)
50  setBinLogX(h.get());
51  if (logy)
52  setBinLogY(h.get());
53  const auto& name = h->GetName();
54  return ibook.book2I(name, h.release());
55  }
56 } // namespace
57 
58 template <typename T>
60 public:
62 
64  ~SiPixelCompareTrackSoAAlpaka() override = default;
65  void bookHistograms(DQMStore::IBooker& ibooker, edm::Run const& iRun, edm::EventSetup const& iSetup) override;
66  void analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) override;
67  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
68 
69 private:
73  const bool useQualityCut_;
75  const float dr2cut_;
96  //1D differences
103 
104  //for matching eff vs region: derive the ratio at harvesting
109 };
110 
111 //
112 // constructors
113 //
114 
115 template <typename T>
117  : tokenSoATrackHost_(consumes<PixelTrackSoA>(iConfig.getParameter<edm::InputTag>("pixelTrackSrcHost"))),
118  tokenSoATrackDevice_(consumes<PixelTrackSoA>(iConfig.getParameter<edm::InputTag>("pixelTrackSrcDevice"))),
119  topFolderName_(iConfig.getParameter<std::string>("topFolderName")),
120  useQualityCut_(iConfig.getParameter<bool>("useQualityCut")),
121  minQuality_(pixelTrack::qualityByName(iConfig.getParameter<std::string>("minQuality"))),
122  dr2cut_(iConfig.getParameter<double>("deltaR2cut")) {}
123 
124 //
125 // -- Analyze
126 //
127 template <typename T>
129  using helper = TracksUtilities<T>;
130  const auto& tsoaHandleHost = iEvent.getHandle(tokenSoATrackHost_);
131  const auto& tsoaHandleDevice = iEvent.getHandle(tokenSoATrackDevice_);
132  if (not tsoaHandleHost or not tsoaHandleDevice) {
133  edm::LogWarning out("SiPixelCompareTrackSoAAlpaka");
134  if (not tsoaHandleHost) {
135  out << "reference (cpu) tracks not found; ";
136  }
137  if (not tsoaHandleDevice) {
138  out << "target (gpu) tracks not found; ";
139  }
140  out << "the comparison will not run.";
141  return;
142  }
143 
144  auto const& tsoaHost = *tsoaHandleHost;
145  auto const& tsoaDevice = *tsoaHandleDevice;
146  auto maxTracksHost = tsoaHost.view().metadata().size(); //this should be same for both?
147  auto maxTracksDevice = tsoaDevice.view().metadata().size(); //this should be same for both?
148  auto const* qualityHost = tsoaHost.view().quality();
149  auto const* qualityDevice = tsoaDevice.view().quality();
150  int32_t nTracksHost = 0;
151  int32_t nTracksDevice = 0;
152  int32_t nLooseAndAboveTracksHost = 0;
153  int32_t nLooseAndAboveTracksHost_matchedDevice = 0;
154  int32_t nLooseAndAboveTracksDevice = 0;
155 
156  //Loop over Device tracks and store the indices of the loose tracks. Whats happens if useQualityCut_ is false?
157  std::vector<int32_t> looseTrkidxDevice;
158  for (int32_t jt = 0; jt < maxTracksDevice; ++jt) {
159  if (helper::nHits(tsoaDevice.view(), jt) == 0)
160  break; // this is a guard
161  if (!(tsoaDevice.view()[jt].pt() > 0.))
162  continue;
163  nTracksDevice++;
164  if (useQualityCut_ && qualityDevice[jt] < minQuality_)
165  continue;
166  nLooseAndAboveTracksDevice++;
167  looseTrkidxDevice.emplace_back(jt);
168  }
169 
170  //Now loop over Host tracks//nested loop for loose gPU tracks
171  for (int32_t it = 0; it < maxTracksHost; ++it) {
172  int nHitsHost = helper::nHits(tsoaHost.view(), it);
173 
174  if (nHitsHost == 0)
175  break; // this is a guard
176 
177  float ptHost = tsoaHost.view()[it].pt();
178  float etaHost = tsoaHost.view()[it].eta();
179  float phiHost = reco::phi(tsoaHost.view(), it);
180  float zipHost = reco::zip(tsoaHost.view(), it);
181  float tipHost = reco::tip(tsoaHost.view(), it);
182 
183  if (!(ptHost > 0.))
184  continue;
185  nTracksHost++;
186  if (useQualityCut_ && qualityHost[it] < minQuality_)
187  continue;
188  nLooseAndAboveTracksHost++;
189  //Now loop over loose Device trk and find the closest in DeltaR//do we need pt cut?
190  const int32_t notFound = -1;
191  int32_t closestTkidx = notFound;
192  float mindr2 = dr2cut_;
193 
194  for (auto gid : looseTrkidxDevice) {
195  float etaDevice = tsoaDevice.view()[gid].eta();
196  float phiDevice = reco::phi(tsoaDevice.view(), gid);
197  float dr2 = reco::deltaR2(etaHost, phiHost, etaDevice, phiDevice);
198  if (dr2 > dr2cut_)
199  continue; // this is arbitrary
200  if (mindr2 > dr2) {
201  mindr2 = dr2;
202  closestTkidx = gid;
203  }
204  }
205 
206  hpt_eta_tkAllHost_->Fill(etaHost, ptHost); //all Host tk
207  hphi_z_tkAllHost_->Fill(phiHost, zipHost);
208  if (closestTkidx == notFound)
209  continue;
210  nLooseAndAboveTracksHost_matchedDevice++;
211 
212  hchi2_->Fill(tsoaHost.view()[it].chi2(), tsoaDevice.view()[closestTkidx].chi2());
213  hCharge_->Fill(reco::charge(tsoaHost.view(), it), reco::charge(tsoaDevice.view(), closestTkidx));
214  hnHits_->Fill(helper::nHits(tsoaHost.view(), it), helper::nHits(tsoaDevice.view(), closestTkidx));
215  hnLayers_->Fill(tsoaHost.view()[it].nLayers(), tsoaDevice.view()[closestTkidx].nLayers());
216  hpt_->Fill(tsoaHost.view()[it].pt(), tsoaDevice.view()[closestTkidx].pt());
217  hptLogLog_->Fill(tsoaHost.view()[it].pt(), tsoaDevice.view()[closestTkidx].pt());
218  heta_->Fill(etaHost, tsoaDevice.view()[closestTkidx].eta());
219  hphi_->Fill(phiHost, reco::phi(tsoaDevice.view(), closestTkidx));
220  hz_->Fill(zipHost, reco::zip(tsoaDevice.view(), closestTkidx));
221  htip_->Fill(tipHost, reco::tip(tsoaDevice.view(), closestTkidx));
222  hptdiffMatched_->Fill(ptHost - tsoaDevice.view()[closestTkidx].pt());
223  hCurvdiffMatched_->Fill((reco::charge(tsoaHost.view(), it) / tsoaHost.view()[it].pt()) -
224  (reco::charge(tsoaDevice.view(), closestTkidx) / tsoaDevice.view()[closestTkidx].pt()));
225  hetadiffMatched_->Fill(etaHost - tsoaDevice.view()[closestTkidx].eta());
226  hphidiffMatched_->Fill(reco::deltaPhi(phiHost, reco::phi(tsoaDevice.view(), closestTkidx)));
227  hzdiffMatched_->Fill(zipHost - reco::zip(tsoaDevice.view(), closestTkidx));
228  htipdiffMatched_->Fill(tipHost - reco::tip(tsoaDevice.view(), closestTkidx));
229  hpt_eta_tkAllHostMatched_->Fill(etaHost, tsoaHost.view()[it].pt()); //matched to gpu
230  hphi_z_tkAllHostMatched_->Fill(etaHost, zipHost);
231  }
232  hnTracks_->Fill(nTracksHost, nTracksDevice);
233  hnLooseAndAboveTracks_->Fill(nLooseAndAboveTracksHost, nLooseAndAboveTracksDevice);
234  hnLooseAndAboveTracks_matched_->Fill(nLooseAndAboveTracksHost, nLooseAndAboveTracksHost_matchedDevice);
235 }
236 
237 //
238 // -- Book Histograms
239 //
240 template <typename T>
242  edm::Run const& iRun,
243  edm::EventSetup const& iSetup) {
244  iBook.cd();
245  iBook.setCurrentFolder(topFolderName_);
246 
247  // clang-format off
248  std::string toRep = "Number of tracks";
249  // FIXME: all the 2D correlation plots are quite heavy in terms of memory consumption, so a as soon as DQM supports THnSparse
250  // these should be moved to a less resource consuming format
251  hnTracks_ = iBook.book2I("nTracks", fmt::format("{} per event; Host; Device",toRep), 501, -0.5, 500.5, 501, -0.5, 500.5);
252  hnLooseAndAboveTracks_ = iBook.book2I("nLooseAndAboveTracks", fmt::format("{} (quality #geq loose) per event; Host; Device",toRep), 501, -0.5, 500.5, 501, -0.5, 500.5);
253  hnLooseAndAboveTracks_matched_ = iBook.book2I("nLooseAndAboveTracks_matched", fmt::format("{} (quality #geq loose) per event; Host; Device",toRep), 501, -0.5, 500.5, 501, -0.5, 500.5);
254 
255  toRep = "Number of all RecHits per track (quality #geq loose)";
256  hnHits_ = iBook.book2I("nRecHits", fmt::format("{};Host;Device",toRep), 15, -0.5, 14.5, 15, -0.5, 14.5);
257 
258  toRep = "Number of all layers per track (quality #geq loose)";
259  hnLayers_ = iBook.book2I("nLayers", fmt::format("{};Host;Device",toRep), 15, -0.5, 14.5, 15, -0.5, 14.5);
260 
261  toRep = "Track (quality #geq loose) #chi^{2}/ndof";
262  hchi2_ = iBook.book2I("nChi2ndof", fmt::format("{};Host;Device",toRep), 40, 0., 20., 40, 0., 20.);
263 
264  toRep = "Track (quality #geq loose) charge";
265  hCharge_ = iBook.book2I("charge",fmt::format("{};Host;Device",toRep),3, -1.5, 1.5, 3, -1.5, 1.5);
266 
267  hpt_ = iBook.book2I("pt", "Track (quality #geq loose) p_{T} [GeV];Host;Device", 200, 0., 200., 200, 0., 200.);
268  hptLogLog_ = make2DIfLog(iBook, true, true, "ptLogLog", "Track (quality #geq loose) p_{T} [GeV];Host;Device", 200, log10(0.5), log10(200.), 200, log10(0.5), log10(200.));
269  heta_ = iBook.book2I("eta", "Track (quality #geq loose) #eta;Host;Device", 30, -3., 3., 30, -3., 3.);
270  hphi_ = iBook.book2I("phi", "Track (quality #geq loose) #phi;Host;Device", 30, -M_PI, M_PI, 30, -M_PI, M_PI);
271  hz_ = iBook.book2I("z", "Track (quality #geq loose) z [cm];Host;Device", 30, -30., 30., 30, -30., 30.);
272  htip_ = iBook.book2I("tip", "Track (quality #geq loose) TIP [cm];Host;Device", 100, -0.5, 0.5, 100, -0.5, 0.5);
273  //1D difference plots
274  hptdiffMatched_ = iBook.book1D("ptdiffmatched", " p_{T} diff [GeV] between matched tracks; #Delta p_{T} [GeV]", 60, -30., 30.);
275  hCurvdiffMatched_ = iBook.book1D("curvdiffmatched", "q/p_{T} diff [GeV] between matched tracks; #Delta q/p_{T} [GeV]", 60, -30., 30.);
276  hetadiffMatched_ = iBook.book1D("etadiffmatched", " #eta diff between matched tracks; #Delta #eta", 160, -0.04 ,0.04);
277  hphidiffMatched_ = iBook.book1D("phidiffmatched", " #phi diff between matched tracks; #Delta #phi", 160, -0.04 ,0.04);
278  hzdiffMatched_ = iBook.book1D("zdiffmatched", " z diff between matched tracks; #Delta z [cm]", 300, -1.5, 1.5);
279  htipdiffMatched_ = iBook.book1D("tipdiffmatched", " TIP diff between matched tracks; #Delta TIP [cm]", 300, -1.5, 1.5);
280  //2D plots for eff
281  hpt_eta_tkAllHost_ = iBook.book2I("ptetatrkAllHost", "Track (quality #geq loose) on Host; #eta; p_{T} [GeV];", 30, -M_PI, M_PI, 200, 0., 200.);
282  hpt_eta_tkAllHostMatched_ = iBook.book2I("ptetatrkAllHostmatched", "Track (quality #geq loose) on Host matched to Device track; #eta; p_{T} [GeV];", 30, -M_PI, M_PI, 200, 0., 200.);
283 
284  hphi_z_tkAllHost_ = iBook.book2I("phiztrkAllHost", "Track (quality #geq loose) on Host; #phi; z [cm];", 30, -M_PI, M_PI, 30, -30., 30.);
285  hphi_z_tkAllHostMatched_ = iBook.book2I("phiztrkAllHostmatched", "Track (quality #geq loose) on Host; #phi; z [cm];", 30, -M_PI, M_PI, 30, -30., 30.);
286 
287 }
288 
289 template<typename T>
291  // monitorpixelTrackSoA
293  desc.add<edm::InputTag>("pixelTrackSrcHost", edm::InputTag("pixelTracksAlpakaSerial"));
294  desc.add<edm::InputTag>("pixelTrackSrcDevice", edm::InputTag("pixelTracksAlpaka"));
295  desc.add<std::string>("topFolderName", "SiPixelHeterogeneous/PixelTrackCompareDeviceVSHost");
296  desc.add<bool>("useQualityCut", true);
297  desc.add<std::string>("minQuality", "loose");
298  desc.add<double>("deltaR2cut", 0.04);
299  descriptions.addWithDefaultLabel(desc);
300 }
301 
305 
constexpr double deltaPhi(double phi1, double phi2)
Definition: deltaPhi.h:26
void addWithDefaultLabel(ParameterSetDescription const &psetDescription)
Definition: helper.py:1
virtual void setCurrentFolder(std::string const &fullpath)
Definition: DQMStore.cc:36
Quality qualityByName(std::string const &name)
ALPAKA_FN_HOST_ACC ALPAKA_FN_INLINE constexpr float zip(ConstView const &tracks, int32_t i)
Definition: TracksSoA.h:90
T eta() const
Definition: PV3DBase.h:73
~SiPixelCompareTrackSoAAlpaka() override=default
ALPAKA_FN_HOST_ACC ALPAKA_FN_INLINE constexpr float charge(ConstView const &tracks, int32_t i)
Definition: TracksSoA.h:73
dqm::reco::DQMStore DQMStore
ALPAKA_FN_HOST_ACC ALPAKA_FN_INLINE constexpr float tip(ConstView const &tracks, int32_t i)
Definition: TracksSoA.h:85
int iEvent
Definition: GenABIO.cc:224
The Signals That Services Can Subscribe To This is based on ActivityRegistry and is current per Services can connect to the signals distributed by the ActivityRegistry in order to monitor the activity of the application Each possible callback has some defined which we here list in angle e< void, edm::EventID const &, edm::Timestamp const & > We also list in braces which AR_WATCH_USING_METHOD_ is used for those or
Definition: Activities.doc:12
void bookHistograms(DQMStore::IBooker &ibooker, edm::Run const &iRun, edm::EventSetup const &iSetup) override
ALPAKA_FN_HOST_ACC ALPAKA_FN_INLINE constexpr float phi(ConstView const &tracks, int32_t i)
Definition: TracksSoA.h:80
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
void analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) override
SiPixelCompareTrackSoAAlpaka(const edm::ParameterSet &)
#define M_PI
constexpr auto deltaR2(const T1 &t1, const T2 &t2) -> decltype(t1.eta())
Definition: deltaR.h:16
const edm::EDGetTokenT< PixelTrackSoA > tokenSoATrackHost_
static const GlobalPoint notFound(0, 0, 0)
HLT enums.
MonitorElement * book1D(TString const &name, TString const &title, int const nchX, double const lowX, double const highX, FUNC onbooking=NOOP())
Definition: DQMStore.h:98
TupleMultiplicity< TrackerTraits > const *__restrict__ uint32_t nHits
The Signals That Services Can Subscribe To This is based on ActivityRegistry h
Helper function to determine trigger accepts.
Definition: Activities.doc:4
const edm::EDGetTokenT< PixelTrackSoA > tokenSoATrackDevice_
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
MonitorElement * book2I(TString const &name, TString const &title, int nchX, double lowX, double highX, int nchY, double lowY, double highY, FUNC onbooking=NOOP())
Definition: DQMStore.h:296
Definition: Run.h:45