LDMX Software
TrackingRecoDQM.cxx
1#include "Tracking/dqm/TrackingRecoDQM.h"
2
3namespace tracking::dqm {
4
5void TrackingRecoDQM::configure(framework::config::Parameters& parameters) {
6 track_collection_ = parameters.get<std::string>("track_collection");
7 truth_collection_ = parameters.get<std::string>("truth_collection");
8 measurement_collection_ =
9 parameters.get<std::string>("measurement_collection");
10 measurement_passname_ = parameters.get<std::string>("measurement_passname");
11
12 ecal_sp_coll_name_ = parameters.get<std::string>("ecal_sp_coll_name");
13 ecal_sp_passname_ = parameters.get<std::string>("ecal_sp_passname");
14 target_sp_coll_name_ = parameters.get<std::string>("target_sp_coll_name");
15 target_sp_passname_ = parameters.get<std::string>("target_sp_passname");
16 truth_passname_ = parameters.get<std::string>("truth_passname");
17 truth_events_passname_ = parameters.get<std::string>("truth_events_passname");
18 track_passname_ = parameters.get<std::string>("track_passname");
19 track_collection_events_passname_ =
20 parameters.get<std::string>("track_collection_events_passname");
21
22 title_ = parameters.get<std::string>("title", "tagger_trk_");
23 track_prob_cut_ = parameters.get<double>("trackProb_cut", 0.5);
24 subdetector_ = parameters.get<std::string>("subdetector", "Tagger");
25 track_states_ = parameters.get<std::vector<std::string>>("track_states", {});
26
27 pidmap_[-321] = PIDBins::kminus;
28 pidmap_[321] = PIDBins::kplus;
29 pidmap_[-211] = PIDBins::piminus;
30 pidmap_[211] = PIDBins::piplus;
31 pidmap_[11] = PIDBins::electron;
32 pidmap_[-11] = PIDBins::positron;
33 pidmap_[2212] = PIDBins::proton;
34 pidmap_[-2212] = PIDBins::antiproton;
35}
36
37void TrackingRecoDQM::analyze(const framework::Event& event) {
38 ldmx_log(trace) << "DQM Reading in:" << track_collection_;
39
40 if (!event.exists(track_collection_, track_collection_events_passname_)) {
41 ldmx_log(error) << "TrackCollection " << track_collection_
42 << " with pass = " << track_collection_events_passname_
43 << " not in event";
44 return;
45 }
46
47 auto tracks{
48 event.getCollection<ldmx::Track>(track_collection_, track_passname_)};
49
50 if (!event.exists(measurement_collection_, measurement_passname_)) {
51 ldmx_log(error) << "Measurement collection " << measurement_collection_
52 << " with pass = " << measurement_passname_
53 << " not in event";
54 return;
55 }
56
57 auto measurements{event.getCollection<ldmx::Measurement>(
58 measurement_collection_, measurement_passname_)};
59
60 // The truth track collection
61 if (event.exists(truth_collection_, truth_events_passname_)) {
62 truth_track_collection_ = std::make_shared<std::vector<ldmx::Track>>(
63 event.getCollection<ldmx::Track>(truth_collection_, truth_passname_));
64 do_truth_comparison_ = true;
65 }
66
67 // The scoring plane hits_
68 if (event.exists(ecal_sp_coll_name_, ecal_sp_passname_)) {
69 ecal_scoring_hits_ = std::make_shared<std::vector<ldmx::SimTrackerHit>>(
70 event.getCollection<ldmx::SimTrackerHit>(ecal_sp_coll_name_,
71 ecal_sp_passname_));
72 }
73
74 if (event.exists(target_sp_coll_name_, target_sp_passname_)) {
75 target_scoring_hits_ = std::make_shared<std::vector<ldmx::SimTrackerHit>>(
76 event.getCollection<ldmx::SimTrackerHit>(target_sp_coll_name_,
77 target_sp_passname_));
78 }
79
80 ldmx_log(debug) << "Do truth comparison::" << do_truth_comparison_;
81
82 if (do_truth_comparison_) {
83 sortTracks(tracks, unique_tracks_, duplicate_tracks_, fake_tracks_);
84 } else {
85 unique_tracks_ = tracks;
86 }
87
88 ldmx_log(debug) << "Filling histograms for " << tracks.size() << " tracks";
89
90 // General Plots
91 histograms_.fill(title_ + "N_tracks", tracks.size());
92
93 if (!unique_tracks_.empty()) {
94 ldmx_log(debug) << "Track Monitoring on " << unique_tracks_.size()
95 << " Unique Tracks";
96 trackMonitoring(unique_tracks_, measurements, title_, true,
97 do_truth_comparison_);
98 }
99
100 // Fakes and duplicates
101 if (!duplicate_tracks_.empty()) {
102 ldmx_log(debug) << "Track Monitoring on " << duplicate_tracks_.size()
103 << " duplicates";
104 trackMonitoring(duplicate_tracks_, measurements, title_ + "dup_", false,
105 false);
106 }
107 if (!fake_tracks_.empty()) {
108 ldmx_log(debug) << "Track Monitoring on " << fake_tracks_.size()
109 << " fakes";
110 trackMonitoring(fake_tracks_, measurements, title_ + "fake_", false, false);
111 }
112
113 // Track Extrapolation to Ecal Monitoring
114 // trackStateMonitoring requires truth_track_collection_ to be available
115 ldmx_log(trace) << "Track Extrapolation to Ecal Monitoring";
116 if (do_truth_comparison_) {
117 if (std::find(track_states_.begin(), track_states_.end(), "target") !=
118 track_states_.end()) {
119 trackStateMonitoring(tracks, ldmx::AtTarget, "target");
120 }
121
122 if (std::find(track_states_.begin(), track_states_.end(), "ecal") !=
123 track_states_.end()) {
124 trackStateMonitoring(tracks, ldmx::AtECAL, "ecal");
125 }
126
127 if (std::find(track_states_.begin(), track_states_.end(), "beamOrigin") !=
128 track_states_.end()) {
129 trackStateMonitoring(tracks, ldmx::AtBeamOrigin, "beamOrigin");
130 }
131 }
132
133 // Technical Efficiency plots
134 if (do_truth_comparison_) {
135 ldmx_log(trace) << "Technical Efficiency plots";
136 efficiencyPlots(tracks, measurements, title_);
137 }
138
139 // Tagger Recoil Matching
140
141 // Clear the vectors
142 ldmx_log(trace) << "Clear the vectors";
143 unique_tracks_.clear();
144 duplicate_tracks_.clear();
145 fake_tracks_.clear();
146}
147
148void TrackingRecoDQM::onProcessEnd() {
149 // Produce the efficiency plots. (TODO::Switch to TEfficiency instead)
150}
151
152void TrackingRecoDQM::efficiencyPlots(
153 const std::vector<ldmx::Track>& tracks,
154 const std::vector<ldmx::Measurement>& measurements,
155 const std::string& title) {
156 // Do all truth track plots - denominator
157
158 histograms_.fill(title + "truth_N_tracks", truth_track_collection_->size());
159 for (auto& truth_trk : *(truth_track_collection_)) {
160 auto truth_phi = truth_trk.getPhi();
161 auto truth_d0 = truth_trk.getD0();
162 auto truth_z0 = truth_trk.getZ0();
163 auto truth_theta = truth_trk.getTheta();
164 auto truth_qop = truth_trk.getQoP();
165 auto truth_p = 1000. / abs(truth_trk.getQoP()); // MeV
166 auto truth_n_hits = truth_trk.getNhits();
167
168 auto truth_mom = truth_trk.getMomentumAtTarget();
169 double truth_pt_beam{0.}, truth_beam_angle{0.};
170 if (truth_mom.size() == 3) {
171 truth_pt_beam =
172 std::sqrt(truth_mom[0] * truth_mom[0] + truth_mom[1] * truth_mom[1]);
173 truth_beam_angle = std::atan2(truth_pt_beam, truth_mom[2]);
174 }
175
176 histograms_.fill(title + "truth_nHits", truth_n_hits);
177 histograms_.fill(title + "truth_d0", truth_d0);
178 histograms_.fill(title + "truth_z0", truth_z0);
179 histograms_.fill(title + "truth_phi", truth_phi);
180 histograms_.fill(title + "truth_theta", truth_theta);
181 histograms_.fill(title + "truth_qop", truth_qop);
182 histograms_.fill(title + "truth_p", truth_p);
183 histograms_.fill(title + "truth_beam_angle", truth_beam_angle);
184
185 if (pidmap_.count(truth_trk.getPdgID()) != 0) {
186 histograms_.fill(title + "truth_PID", pidmap_[truth_trk.getPdgID()]);
187
188 // TODO do this properly.
189
190 if (pidmap_[truth_trk.getPdgID()] == PIDBins::kminus) {
191 histograms_.fill(title + "truth_kminus_p", truth_p);
192 }
193
194 if (pidmap_[truth_trk.getPdgID()] == PIDBins::kplus) {
195 histograms_.fill(title + "truth_kplus_p", truth_p);
196 }
197
198 if (pidmap_[truth_trk.getPdgID()] == PIDBins::piminus) {
199 histograms_.fill(title + "truth_piminus_p", truth_p);
200 }
201
202 if (pidmap_[truth_trk.getPdgID()] == PIDBins::piplus) {
203 histograms_.fill(title + "truth_piplus_p", truth_p);
204 }
205
206 if (pidmap_[truth_trk.getPdgID()] == PIDBins::electron) {
207 histograms_.fill(title + "truth_electron_p", truth_p);
208 }
209
210 if (pidmap_[truth_trk.getPdgID()] == PIDBins::positron) {
211 histograms_.fill(title + "truth_positron_p", truth_p);
212 }
213
214 if (pidmap_[truth_trk.getPdgID()] == PIDBins::proton) {
215 histograms_.fill(title + "truth_proton_p", truth_p);
216 }
217 }
218
219 } // loop on truth tracks
220
221 for (auto& track : tracks) {
222 // Match the tracks to truth
223 ldmx::Track* truth_trk = nullptr;
224
225 auto it = std::find_if(truth_track_collection_->begin(),
226 truth_track_collection_->end(),
227 [&](const ldmx::Track& tt) {
228 return tt.getTrackID() == track.getTrackID();
229 });
230
231 double track_truth_prob = track.getTruthProb();
232
233 if (it != truth_track_collection_->end() &&
234 track_truth_prob >= track_prob_cut_)
235 truth_trk = &(*it);
236
237 // Match not found
238 if (!truth_trk) return;
239
240 auto truth_phi = truth_trk->getPhi();
241 auto truth_d0 = truth_trk->getD0();
242 auto truth_z0 = truth_trk->getZ0();
243 auto truth_theta = truth_trk->getTheta();
244 auto truth_qop = truth_trk->getQoP();
245 auto truth_p = 1000. / abs(truth_trk->getQoP()); // MeV
246
247 auto truth_mom = truth_trk->getMomentumAtTarget();
248 double truth_pt_beam{0.}, truth_beam_angle{0.};
249 if (truth_mom.size() == 3) {
250 truth_pt_beam =
251 std::sqrt(truth_mom[0] * truth_mom[0] + truth_mom[1] * truth_mom[1]);
252 truth_beam_angle = std::atan2(truth_pt_beam, truth_mom[2]);
253 }
254
255 // Fill reco plots for efficiencies - numerator. The quantities are truth
256 histograms_.fill(title + "match_prob", track_truth_prob);
257 histograms_.fill(title + "match_d0", truth_d0);
258 histograms_.fill(title + "match_z0", truth_z0);
259 histograms_.fill(title + "match_phi", truth_phi);
260 histograms_.fill(title + "match_theta", truth_theta);
261 histograms_.fill(title + "match_p", truth_p);
262 histograms_.fill(title + "match_qop", truth_qop);
263 histograms_.fill(title + "match_beam_angle", truth_beam_angle);
264 histograms_.fill(title + "match_nHits", measurements.size());
265 auto dedx_measurements = track.getDedxMeasurements();
266 auto measurement_idxs = track.getMeasurementsIdxs();
267 for (size_t i = 0; i < measurement_idxs.size(); ++i) {
268 histograms_.fill(title + "match_layers_hit",
269 measurements.at(measurement_idxs[i]).getLayer());
270 // Add histogram of the measurement dE/dx
271 if (i < dedx_measurements.size()) {
272 histograms_.fill(title + "match_measurement_dedx",
273 dedx_measurements[i]);
274 }
275 }
276
277 // For some particles
278
279 if (pidmap_.count(truth_trk->getPdgID()) != 0) {
280 histograms_.fill(title + "match_PID", pidmap_[truth_trk->getPdgID()]);
281
282 // TODO do this properly.
283
284 if (pidmap_[truth_trk->getPdgID()] == PIDBins::kminus) {
285 histograms_.fill(title + "match_kminus_p", truth_p);
286 }
287
288 if (pidmap_[truth_trk->getPdgID()] == PIDBins::kplus) {
289 histograms_.fill(title + "match_kplus_p", truth_p);
290 }
291
292 if (pidmap_[truth_trk->getPdgID()] == PIDBins::piminus) {
293 histograms_.fill(title + "match_piminus_p", truth_p);
294 }
295
296 if (pidmap_[truth_trk->getPdgID()] == PIDBins::piplus) {
297 histograms_.fill(title + "match_piplus_p", truth_p);
298 }
299
300 if (pidmap_[truth_trk->getPdgID()] == PIDBins::electron) {
301 histograms_.fill(title + "match_electron_p", truth_p);
302 }
303
304 if (pidmap_[truth_trk->getPdgID()] == PIDBins::positron) {
305 histograms_.fill(title + "match_positron_p", truth_p);
306 }
307
308 if (pidmap_[truth_trk->getPdgID()] == PIDBins::proton) {
309 histograms_.fill(title + "match_proton_p", truth_p);
310 }
311 }
312 } // Loop on tracks
313
314} // Efficiency plots
315
316void TrackingRecoDQM::trackMonitoring(
317 const std::vector<ldmx::Track>& tracks,
318 const std::vector<ldmx::Measurement>& measurements, const std::string title,
319 const bool& doDetail, const bool& doTruth) {
320 for (auto& track : tracks) {
321 // Perigee track parameters
322 auto trk_d0 = track.getD0();
323 auto trk_z0 = track.getZ0();
324 auto trk_qop = track.getQoP();
325 auto trk_theta = track.getTheta();
326 auto trk_phi = track.getPhi();
327 auto trk_p = 1000. / abs(trk_qop); // MeV
328 auto dedx_measurements = track.getDedxMeasurements();
329 auto measurement_idxs = track.getMeasurementsIdxs();
330 for (size_t i = 0; i < measurement_idxs.size(); ++i) {
331 histograms_.fill(title + "layers_hit",
332 measurements.at(measurement_idxs[i]).getLayer());
333 // Add histogram of the measurement dE/dx
334 if (i < dedx_measurements.size()) {
335 histograms_.fill(title + "measurement_dedx", dedx_measurements[i]);
336 }
337 }
338
339 // Per-layer unbiased U-residuals — only filled for the main unique-track
340 // call (title == title_). Fakes/duplicates use a different title prefix
341 // and do not have these histograms declared.
342 //
343 // Algebraic leave-one-out unbiasing (NIM A 262, 444, 1987):
344 // r_smooth = m - x_smooth
345 // denom = V - C (V = cov_uu, C = smoothed cov[loc0,loc0])
346 // r_ubs = V / denom * r_smooth
347 // pull = r_smooth / sqrt(denom)
348 if (title == title_) {
349 const auto& sm_loc0 = track.getSmoothedLoc0();
350 const auto& sm_cov = track.getSmoothedCovLoc0();
351 for (size_t i = 0; i < measurement_idxs.size(); ++i) {
352 if (i >= sm_loc0.size()) break;
353 const auto& meas = measurements.at(measurement_idxs[i]);
354 int layer = meas.getLayer();
355 float meas_u = meas.getLocalPosition()[0];
356 float V = meas.getLocalCovariance()[0]; // cov_uu
357 float C = sm_cov[i]; // smoothed cov[loc0, loc0]
358 float denom = V - C;
359 if (denom <= 0.f) continue; // degenerate — skip
360 float r_smooth = meas_u - sm_loc0[i];
361 float res_ubs = r_smooth * V / denom;
362 float pull_ubs = r_smooth / std::sqrt(denom);
363 histograms_.fill(title_ + "unbiased_res_u_l" + std::to_string(layer),
364 res_ubs);
365 histograms_.fill(title_ + "unbiased_pull_u_l" + std::to_string(layer),
366 pull_ubs);
367 }
368 }
369
370 auto trk_mom = track.getMomentumAtTarget();
371 // getMomentumAtTarget() returns MeV (LDMX convention)
372 double px_ldmx{0.}, py_ldmx{0.}, pz_ldmx{0.};
373 double pt_bending{0.}, pt_beam{0.};
374 if (trk_mom.size() == 3) {
375 px_ldmx = trk_mom[0]; // MeV
376 py_ldmx = trk_mom[1];
377 pz_ldmx = trk_mom[2];
378 // Bending-plane pT: horizontal (x_ldmx) + downstream (z_ldmx)
379 pt_bending = std::sqrt(px_ldmx * px_ldmx + pz_ldmx * pz_ldmx);
380 // Transverse pT perpendicular to beam: horizontal + vertical
381 pt_beam = std::sqrt(px_ldmx * px_ldmx + py_ldmx * py_ldmx);
382 }
383
384 // Covariance matrix
385 Acts::BoundSquareMatrix cov =
386 tracking::sim::utils::unpackCov(track.getPerigeeCov());
387
388 double sigmad0 = sqrt(
389 cov(Acts::BoundIndices::eBoundLoc0, Acts::BoundIndices::eBoundLoc0));
390 double sigmaz0 = sqrt(
391 cov(Acts::BoundIndices::eBoundLoc1, Acts::BoundIndices::eBoundLoc1));
392 double sigmaphi =
393 sqrt(cov(Acts::BoundIndices::eBoundPhi, Acts::BoundIndices::eBoundPhi));
394 double sigmatheta = sqrt(
395 cov(Acts::BoundIndices::eBoundTheta, Acts::BoundIndices::eBoundTheta));
396 double sigmaqop = sqrt(cov(Acts::BoundIndices::eBoundQOverP,
397 Acts::BoundIndices::eBoundQOverP));
398 double sigmap =
399 (1000. / trk_qop) * (1000. / trk_qop) * sigmaqop / 1000.; // MeV
400
401 histograms_.fill(title + "d0", trk_d0);
402 histograms_.fill(title + "z0", trk_z0);
403 histograms_.fill(title + "qop", trk_qop);
404 histograms_.fill(title + "phi", trk_phi);
405 histograms_.fill(title + "theta", trk_theta);
406 histograms_.fill(title + "p", trk_p);
407
408 if (doDetail) {
409 histograms_.fill(title + "px", px_ldmx);
410 histograms_.fill(title + "py", py_ldmx);
411 histograms_.fill(title + "pz", pz_ldmx);
412
413 histograms_.fill(title + "pt_bending", pt_bending);
414 histograms_.fill(title + "pt_beam", pt_beam);
415
416 histograms_.fill(title + "nHits", track.getNhits());
417 histograms_.fill(title + "Chi2", track.getChi2());
418 histograms_.fill(title + "ndf", track.getNdf());
419 histograms_.fill(title + "Chi2_per_ndf",
420 track.getChi2() / track.getNdf());
421 histograms_.fill(title + "nShared", track.getNsharedHits());
422
423 histograms_.fill(title + "d0_err", sigmad0);
424 histograms_.fill(title + "z0_err", sigmaz0);
425 histograms_.fill(title + "phi_err", sigmaphi);
426 histograms_.fill(title + "theta_err", sigmatheta);
427 histograms_.fill(title + "qop_err", sigmaqop);
428 histograms_.fill(title + "p_err", sigmap);
429
430 // 2D Error plots (p in MeV)
431 histograms_.fill(title + "d0_err_vs_p", trk_p, sigmad0);
432 histograms_.fill(title + "z0_err_vs_p", trk_p, sigmaz0);
433 histograms_.fill(title + "p_err_vs_p", trk_p, sigmap);
434
435 if (track.getNhits() == 8)
436 histograms_.fill(title + "p_err_vs_p_8hits", trk_p, sigmap);
437 else if (track.getNhits() == 9)
438 histograms_.fill(title + "p_err_vs_p_9hits", trk_p, sigmap);
439 else if (track.getNhits() == 10)
440 histograms_.fill(title + "p_err_vs_p_10hits", trk_p, sigmap);
441 }
442
443 if (doTruth) {
444 // Match to the truth track
445 ldmx::Track* truth_trk = nullptr;
446
447 auto it = std::find_if(truth_track_collection_->begin(),
448 truth_track_collection_->end(),
449 [&](const ldmx::Track& tt) {
450 return tt.getTrackID() == track.getTrackID();
451 });
452
453 double track_truth_prob = track.getTruthProb();
454
455 if (it != truth_track_collection_->end() &&
456 track_truth_prob >= track_prob_cut_)
457 truth_trk = &(*it);
458
459 // Found matched track
460 if (truth_trk) {
461 auto truth_d0 = truth_trk->getD0();
462 auto truth_z0 = truth_trk->getZ0();
463 auto truth_phi = truth_trk->getPhi();
464 auto truth_theta = truth_trk->getTheta();
465 auto truth_qop = truth_trk->getQoP();
466 auto truth_p = 1000. / abs(truth_trk->getQoP()); // MeV
467 auto truth_mom = truth_trk->getMomentumAtTarget();
468 // getMomentumAtTarget() returns MeV (LDMX convention)
469 double truth_pt_beam{0.};
470 if (truth_mom.size() == 3) {
471 truth_pt_beam = std::sqrt(truth_mom[0] * truth_mom[0] +
472 truth_mom[1] * truth_mom[1]);
473 }
474
475 // histograms_.fill(title+"truth_d0", truth_d0);
476 // histograms_.fill(title+"truth_z0", truth_z0);
477 // histograms_.fill(title+"truth_phi", truth_phi);
478 // histograms_.fill(title+"truth_theta",truth_theta);
479 // histograms_.fill(title+"truth_qop", truth_qop);
480 // histograms_.fill(title+"truth_p", truth_p);
481
482 double res_d0 = trk_d0 - truth_d0;
483 double res_z0 = trk_z0 - truth_z0;
484 double res_phi = trk_phi - truth_phi;
485 double res_theta = trk_theta - truth_theta;
486 double res_qop = trk_qop - truth_qop;
487 double res_p = trk_p - truth_p;
488 double res_pt_beam = pt_beam - truth_pt_beam;
489
490 histograms_.fill(title + "res_d0", res_d0);
491 histograms_.fill(title + "res_z0", res_z0);
492 histograms_.fill(title + "res_phi", res_phi);
493 histograms_.fill(title + "res_theta", res_theta);
494 histograms_.fill(title + "res_qop", res_qop);
495 histograms_.fill(title + "res_p", res_p);
496 histograms_.fill(title + "res_pt_beam", res_pt_beam);
497
498 double pull_d0 = res_d0 / sigmad0;
499 double pull_z0 = res_z0 / sigmaz0;
500 double pull_phi = res_phi / sigmaphi;
501 double pull_theta = res_theta / sigmatheta;
502 double pull_qop = res_qop / sigmaqop;
503 double pull_p = res_p / sigmap;
504
505 histograms_.fill(title + "pull_d0", pull_d0);
506 histograms_.fill(title + "pull_z0", pull_z0);
507 histograms_.fill(title + "pull_phi", pull_phi);
508 histograms_.fill(title + "pull_theta", pull_theta);
509 histograms_.fill(title + "pull_qop", pull_qop);
510 histograms_.fill(title + "pull_p", pull_p);
511
512 // Error plots from residuals
513
514 histograms_.fill(title + "res_p_vs_p", truth_p, res_p);
515
516 histograms_.fill(title + "res_qop_vs_p", truth_p, res_qop);
517 histograms_.fill(title + "res_d0_vs_p", truth_p, res_d0);
518 histograms_.fill(title + "res_z0_vs_p", truth_p, res_z0);
519 histograms_.fill(title + "res_phi_vs_p", truth_p, res_phi);
520 histograms_.fill(title + "res_theta_vs_p", truth_p, res_theta);
521
522 histograms_.fill(title + "pull_qop_vs_p", truth_p, pull_qop);
523 histograms_.fill(title + "pull_d0_vs_p", truth_p, pull_d0);
524 histograms_.fill(title + "pull_z0_vs_p", truth_p, pull_z0);
525 histograms_.fill(title + "pull_phi_vs_p", truth_p, pull_phi);
526 histograms_.fill(title + "pull_theta_vs_p", truth_p, pull_theta);
527
528 if (track.getNhits() == 8)
529 histograms_.fill(title + "res_p_vs_p_8hits", truth_p, res_p);
530 else if (track.getNhits() == 9)
531 histograms_.fill(title + "res_p_vs_p_9hits", truth_p, res_p);
532 else if (track.getNhits() == 10)
533 histograms_.fill(title + "res_p_vs_p_10hits", truth_p, res_p);
534
535 histograms_.fill(title + "res_pt_beam_vs_p", truth_pt_beam,
536 res_pt_beam);
537
538 } // found matched track
539 } // do TruthComparison
540 } // loop on tracks
541
542} // Track Monitoring
543
544void TrackingRecoDQM::trackStateMonitoring(
545 const std::vector<ldmx::Track>& tracks, ldmx::TrackStateType ts_type,
546 const std::string& ts_title) {
547 for (auto& track : tracks) {
548 // Match the tracks to truth
549 ldmx::Track* truth_trk = nullptr;
550
551 auto it = std::find_if(truth_track_collection_->begin(),
552 truth_track_collection_->end(),
553 [&](const ldmx::Track& tt) {
554 return tt.getTrackID() == track.getTrackID();
555 });
556
557 double track_truth_prob = track.getTruthProb();
558
559 if (it != truth_track_collection_->end() &&
560 track_truth_prob >= track_prob_cut_)
561 truth_trk = &(*it);
562
563 // Match not found, skip track
564 if (!truth_trk) continue;
565
566 auto trk_ts = track.getTrackState(ts_type);
567 auto truth_ts = truth_trk->getTrackState(ts_type);
568
569 if (!trk_ts.has_value()) continue;
570 if (!truth_ts.has_value()) continue;
571
572 const ldmx::Track::TrackState& target_state = trk_ts.value();
573 const ldmx::Track::TrackState& truth_target_state = truth_ts.value();
574
575 // Check that the covariance is filled
576 if (target_state.pos_mom_cov_.size() < 21) continue;
577
578 // Sigma from Cartesian position covariance (LDMX frame):
579 // pos_mom_cov_ upper-triangle layout: xx=0, yy=6
580 double sigmaloc0 = std::sqrt(target_state.pos_mom_cov_[0]); // sigma_x
581 double sigmaloc1 = std::sqrt(target_state.pos_mom_cov_[6]); // sigma_y
582
583 double trk_qop = track.getQoP();
584 double trk_p = 1000. / abs(trk_qop); // MeV
585
586 // loc0/loc1 = x/y position at the surface in LDMX global frame
587 double track_state_loc0 = target_state.pos_[0];
588 double track_state_loc1 = target_state.pos_[1];
589
590 double truth_state_loc0 = truth_target_state.pos_[0];
591 double truth_state_loc1 = truth_target_state.pos_[1];
592
593 histograms_.fill(title_ + "trk_" + ts_title + "_loc0", track_state_loc0);
594 histograms_.fill(title_ + "trk_" + ts_title + "_loc1", track_state_loc1);
595 histograms_.fill(title_ + ts_title + "_truth_loc0", truth_state_loc0);
596 histograms_.fill(title_ + ts_title + "_truth_loc1", truth_state_loc1);
597
598 // TH1F residuals
599 histograms_.fill(
600 title_ + "trk_" + ts_title + "_loc0-truth_" + ts_title + "_loc0",
601 track_state_loc0 - truth_state_loc0);
602 histograms_.fill(
603 title_ + "trk_" + ts_title + "_loc1-truth_" + ts_title + "_loc1",
604 track_state_loc1 - truth_state_loc1);
605
606 // TH1F The pulls of loc0 and loc1
607 histograms_.fill(title_ + ts_title + "_Pulls_of_loc0",
608 (track_state_loc0 - truth_state_loc0) / sigmaloc0);
609 histograms_.fill(title_ + ts_title + "_Pulls_of_loc1",
610 (track_state_loc1 - truth_state_loc1) / sigmaloc1);
611
612 // TODO:: TH1F The pulls of phi, theta, qop
613
614 // TH2F residual vs Nhits
615 histograms_.fill(title_ + ts_title + "_res_loc0-vs-N_hits",
616 track.getNhits(), track_state_loc0 - truth_state_loc0);
617 histograms_.fill(title_ + ts_title + "_res_loc1-vs-N_hits",
618 track.getNhits(), track_state_loc1 - truth_state_loc1);
619
620 // TH2F pulls vs Nhits
621 histograms_.fill(title_ + ts_title + "_pulls_loc0-vs-N_hits",
622 track.getNhits(),
623 (track_state_loc0 - truth_state_loc0) / sigmaloc0);
624 histograms_.fill(title_ + ts_title + "_pulls_loc1-vs-N_hits",
625 track.getNhits(),
626 (track_state_loc1 - truth_state_loc1) / sigmaloc1);
627
628 // TH2F residual vs trk_p
629 histograms_.fill(title_ + ts_title + "_res_loc0-vs-trk_p", trk_p,
630 track_state_loc0 - truth_state_loc0);
631 histograms_.fill(title_ + ts_title + "_res_loc1-vs-trk_p", trk_p,
632 track_state_loc1 - truth_state_loc1);
633
634 // TH2F pulls vs trk_p
635 histograms_.fill(title_ + ts_title + "_pulls_loc0-vs-trk_p", trk_p,
636 (track_state_loc0 - truth_state_loc0) / sigmaloc0);
637 histograms_.fill(title_ + ts_title + "_pulls_loc1-vs-trk_p", trk_p,
638 (track_state_loc1 - truth_state_loc1) / sigmaloc1);
639
640 } // loop on tracks
641}
642
643void TrackingRecoDQM::sortTracks(const std::vector<ldmx::Track>& tracks,
644 std::vector<ldmx::Track>& uniqueTracks,
645 std::vector<ldmx::Track>& duplicateTracks,
646 std::vector<ldmx::Track>& fakeTracks) {
647 // Create a copy of the const vector so we can sort it
648 std::vector<ldmx::Track> sorted_tracks = tracks;
649
650 // Sort the vector of Track objects based on their trackID member
651 std::sort(sorted_tracks.begin(), sorted_tracks.end(),
652 [](ldmx::Track& t1, ldmx::Track& t2) {
653 return t1.getTrackID() < t2.getTrackID();
654 });
655
656 // Loop over the sorted vector of Track objects
657 for (size_t i = 0; i < sorted_tracks.size(); i++) {
658 if (sorted_tracks[i].getTruthProb() < track_prob_cut_)
659 fakeTracks.push_back(sorted_tracks[i]);
660 else { // not a fake track
661 // If this is the first Track object with this trackID, add it to the
662 // uniqueTracks vector directly
663 if (uniqueTracks.size() == 0 ||
664 sorted_tracks[i].getTrackID() != sorted_tracks[i - 1].getTrackID()) {
665 uniqueTracks.push_back(sorted_tracks[i]);
666 }
667 // Otherwise, add it to the duplicateTracks vector if its truthProb is
668 // lower than the existing Track object Otherwise, if the truthProbability
669 // is higher than the track stored in uniqueTracks, put it in uniqueTracks
670 // and move the uniqueTracks.back to duplicateTracks.
671 else if (sorted_tracks[i].getTruthProb() >
672 uniqueTracks.back().getTruthProb()) {
673 duplicateTracks.push_back(uniqueTracks.back());
674 uniqueTracks.back() = sorted_tracks[i];
675 }
676 // Otherwise, add it to the duplicateTracks vector
677 else {
678 duplicateTracks.push_back(sorted_tracks[i]);
679 }
680 } // a real track
681 } // loop on sorted tracks
682 // The total number of elements in the uniqueTracks and duplicateTracks
683 // vectors should be equal to the number of elements in the original tracks
684 // vector
685 if (uniqueTracks.size() + duplicateTracks.size() + fakeTracks.size() !=
686 tracks.size()) {
687 std::cerr << "Error: unique and duplicate tracks vectors do not add up to "
688 "original tracks vector";
689 return;
690 }
691
692 // Iterate through the uniqueTracks vector and duplicateTracks vector
693 ldmx_log(trace) << "Unique tracks:";
694 for (const ldmx::Track& track : uniqueTracks) {
695 ldmx_log(trace) << "\tTrack ID: " << track.getTrackID()
696 << ", Truth Prob: " << track.getTruthProb();
697 }
698 ldmx_log(trace) << "Duplicate tracks:";
699 for (const ldmx::Track& track : duplicateTracks) {
700 ldmx_log(trace) << "\tTrack ID: " << track.getTrackID()
701 << ", Truth Prob: " << track.getTruthProb();
702 }
703 ldmx_log(trace) << "Fake tracks:";
704 for (const ldmx::Track& track : fakeTracks) {
705 ldmx_log(trace) << "\tTrack ID: " << track.getTrackID()
706 << ", Truth Prob: " << track.getTruthProb();
707 }
708}
709} // namespace tracking::dqm
710
711DECLARE_ANALYZER(tracking::dqm::TrackingRecoDQM)
DECLARE_ANALYZER
#define DECLARE_ANALYZER(CLASS)
Macro which allows the framework to construct an analyzer given its name during configuration.
Definition EventProcessor.h:363
framework::EventProcessor::histograms_
HistogramPool histograms_
helper object for making and filling histograms
Definition EventProcessor.h:221
framework::Event
Implements an event buffer system for storing event data.
Definition Event.h:42
framework::Event::exists
bool exists(const std::string &name, const std::string &passName, bool unique=true) const
Check for the existence of an object or collection with the given name and pass name in the event.
Definition Event.cxx:105
framework::HistogramPool::fill
void fill(const std::string &name, const T &val)
Fill a 1D histogram.
Definition HistogramPool.h:367
framework::config::Parameters
Class encapsulating parameters for configuring a processor.
Definition Parameters.h:29
framework::config::Parameters::get
const T & get(const std::string &name) const
Retrieve the parameter of the given name.
Definition Parameters.h:78
ldmx::SimTrackerHit
Represents a simulated tracker hit in the simulation.
Definition SimTrackerHit.h:24
ldmx::Track
Implementation of a track object.
Definition Track.h:53
ldmx::Track::getMomentumAtTarget
std::vector< double > getMomentumAtTarget() const
Returns the momentum (px, py, pz) in MeV in the LDMX global frame from the AtTarget TrackState.
Definition Track.h:218
tracking::dqm::TrackingRecoDQM::trackStateMonitoring
void trackStateMonitoring(const std::vector< ldmx::Track > &tracks, ldmx::TrackStateType ts_type, const std::string &ts_title)
Monitoring plots for tracks extrapolated to the ECAL Scoring plane.
Definition TrackingRecoDQM.cxx:544
tracking::dqm::TrackingRecoDQM::onProcessEnd
void onProcessEnd() override
Callback for the EventProcessor to take any necessary action when the processing of events finishes,...
Definition TrackingRecoDQM.cxx:148
tracking::dqm::TrackingRecoDQM::configure
void configure(framework::config::Parameters &parameters) override
Configure the analyzer using the given user specified parameters.
Definition TrackingRecoDQM.cxx:5
tracking::dqm::TrackingRecoDQM::analyze
void analyze(const framework::Event &event) override
Process the event and make histograms or summaries.
Definition TrackingRecoDQM.cxx:37
Generated by doxygen 1.12.0