LDMX Software
TrackingRecoDQM.cxx
1#include "Tracking/dqm/TrackingRecoDQM.h"
2
3#include <algorithm>
4#include <iostream>
5
6#include "Tracking/Sim/TrackingUtils.h"
7
8namespace tracking::dqm {
9
10void TrackingRecoDQM::configure(framework::config::Parameters& parameters) {
11 trackCollection_ =
12 parameters.getParameter<std::string>("track_collection", "TaggerTracks");
13 truthCollection_ = parameters.getParameter<std::string>("truth_collection",
14 "TaggerTruthTracks");
15 title_ = parameters.getParameter<std::string>("title", "tagger_trk_");
16 trackProb_cut_ = parameters.getParameter<double>("trackProb_cut", 0.5);
17 subdetector_ = parameters.getParameter<std::string>("subdetector", "Tagger");
18 trackStates_ =
19 parameters.getParameter<std::vector<std::string>>("trackStates", {});
20 measurementCollection_ = parameters.getParameter<std::string>(
21 "measurement_collection", "DigiTaggerSimHits");
22
23 ldmx_log(info) << "Track Collection " << trackCollection_ << std::endl;
24 ldmx_log(info) << "Truth Collection " << truthCollection_ << std::endl;
25
26 pidmap[-321] = PIDBins::kminus;
27 pidmap[321] = PIDBins::kplus;
28 pidmap[-211] = PIDBins::piminus;
29 pidmap[211] = PIDBins::piplus;
30 pidmap[11] = PIDBins::electron;
31 pidmap[-11] = PIDBins::positron;
32 pidmap[2212] = PIDBins::proton;
33 pidmap[-2212] = PIDBins::antiproton;
34}
35
36void TrackingRecoDQM::analyze(const framework::Event& event) {
37 ldmx_log(debug) << "DQM Reading in::" << trackCollection_ << std::endl;
38
39 if (!event.exists(trackCollection_)) {
40 ldmx_log(error) << "ERROR:: trackCollection " << trackCollection_
41 << " not in event" << std::endl;
42 return;
43 }
44 auto tracks{event.getCollection<ldmx::Track>(trackCollection_)};
45 auto measurements{
46 event.getCollection<ldmx::Measurement>(measurementCollection_)};
47 // The truth track collection
48 if (event.exists(truthCollection_)) {
49 truthTrackCollection_ = std::make_shared<ldmx::Tracks>(
50 event.getCollection<ldmx::Track>(truthCollection_));
51 doTruthComparison = true;
52 }
53
54 // The scoring plane hits
55 if (event.exists("EcalScoringPlaneHits")) {
56 ecal_scoring_hits_ = std::make_shared<std::vector<ldmx::SimTrackerHit>>(
57 event.getCollection<ldmx::SimTrackerHit>("EcalScoringPlaneHits"));
58 }
59
60 if (event.exists("TargetScoringPlaneHits")) {
61 target_scoring_hits_ = std::make_shared<std::vector<ldmx::SimTrackerHit>>(
62 event.getCollection<ldmx::SimTrackerHit>("TargetScoringPlaneHits"));
63 }
64
65 ldmx_log(debug) << "Do truth comparison::" << doTruthComparison << std::endl;
66
67 if (doTruthComparison) {
68 sortTracks(tracks, uniqueTracks_, duplicateTracks_, fakeTracks_);
69 } else {
70 uniqueTracks_ = tracks;
71 }
72
73 ldmx_log(debug) << "Filling histograms " << std::endl;
74
75 // General Plots
76 histograms_.fill(title_ + "N_tracks", tracks.size());
77
78 ldmx_log(debug) << "Track Monitoring on Unique Tracks" << std::endl;
79
80 TrackMonitoring(uniqueTracks_, measurements, title_, true, true);
81
82 ldmx_log(debug) << "Track Monitoring on duplicates and fakes" << std::endl;
83 // Fakes and duplicates
84 TrackMonitoring(duplicateTracks_, measurements, title_ + "dup_", false,
85 false);
86 TrackMonitoring(fakeTracks_, measurements, title_ + "fake_", false, false);
87
88 // Track Extrapolation to Ecal Monitoring
89
90 if (std::find(trackStates_.begin(), trackStates_.end(), "target") !=
91 trackStates_.end()) {
92 TrackStateMonitoring(tracks, ldmx::TrackStateType::AtTarget, "target");
93 }
94
95 if (std::find(trackStates_.begin(), trackStates_.end(), "ecal") !=
96 trackStates_.end()) {
97 TrackStateMonitoring(tracks, ldmx::TrackStateType::AtECAL, "ecal");
98 }
99
100 if (std::find(trackStates_.begin(), trackStates_.end(), "beamOrigin") !=
101 trackStates_.end()) {
102 TrackStateMonitoring(tracks, ldmx::TrackStateType::AtBeamOrigin,
103 "beamOrigin");
104 }
105
106 // Technical Efficiency plots
107 EfficiencyPlots(tracks, measurements, title_);
108
109 // Tagger Recoil Matching
110
111 // Clear the vectors
112 uniqueTracks_.clear();
113 duplicateTracks_.clear();
114 fakeTracks_.clear();
115}
116
117void TrackingRecoDQM::onProcessEnd() {
118 // Produce the efficiency plots. (TODO::Switch to TEfficiency instead)
119}
120
121void TrackingRecoDQM::EfficiencyPlots(
122 const std::vector<ldmx::Track>& tracks,
123 const std::vector<ldmx::Measurement>& measurements,
124 const std::string& title) {
125 // Do all truth track plots - denominator
126
127 histograms_.fill(title + "truth_N_tracks", truthTrackCollection_->size());
128 for (auto& truth_trk : *(truthTrackCollection_)) {
129 double truth_phi = truth_trk.getPhi();
130 double truth_d0 = truth_trk.getD0();
131 double truth_z0 = truth_trk.getZ0();
132 double truth_theta = truth_trk.getTheta();
133 double truth_qop = truth_trk.getQoP();
134 double truth_p = 1. / abs(truth_trk.getQoP());
135 int truth_nHits = truth_trk.getNhits();
136
137 std::vector<double> truth_mom = truth_trk.getMomentum();
138
139 // Polar angle
140 // The momentum in the plane transverse wrt the beam axis
141 double truth_pt_beam =
142 std::sqrt(truth_mom[1] * truth_mom[1] + truth_mom[2] * truth_mom[2]);
143
144 double truth_beam_angle = std::atan2(truth_pt_beam, truth_mom[0]);
145
146 histograms_.fill(title + "truth_nHits", truth_nHits);
147 histograms_.fill(title + "truth_d0", truth_d0);
148 histograms_.fill(title + "truth_z0", truth_z0);
149 histograms_.fill(title + "truth_phi", truth_phi);
150 histograms_.fill(title + "truth_theta", truth_theta);
151 histograms_.fill(title + "truth_qop", truth_qop);
152 histograms_.fill(title + "truth_p", truth_p);
153 histograms_.fill(title + "truth_beam_angle", truth_beam_angle);
154
155 if (pidmap.count(truth_trk.getPdgID()) != 0) {
156 histograms_.fill(title + "truth_PID", pidmap[truth_trk.getPdgID()]);
157
158 // TODO do this properly.
159
160 if (pidmap[truth_trk.getPdgID()] == PIDBins::kminus) {
161 histograms_.fill(title + "truth_kminus_p", truth_p);
162 }
163
164 if (pidmap[truth_trk.getPdgID()] == PIDBins::kplus) {
165 histograms_.fill(title + "truth_kplus_p", truth_p);
166 }
167
168 if (pidmap[truth_trk.getPdgID()] == PIDBins::piminus) {
169 histograms_.fill(title + "truth_piminus_p", truth_p);
170 }
171
172 if (pidmap[truth_trk.getPdgID()] == PIDBins::piplus) {
173 histograms_.fill(title + "truth_piplus_p", truth_p);
174 }
175
176 if (pidmap[truth_trk.getPdgID()] == PIDBins::electron) {
177 histograms_.fill(title + "truth_electron_p", truth_p);
178 }
179
180 if (pidmap[truth_trk.getPdgID()] == PIDBins::positron) {
181 histograms_.fill(title + "truth_positron_p", truth_p);
182 }
183
184 if (pidmap[truth_trk.getPdgID()] == PIDBins::proton) {
185 histograms_.fill(title + "truth_proton_p", truth_p);
186 }
187 }
188
189 } // loop on truth tracks
190
191 for (auto& track : tracks) {
192 // Match the tracks to truth
193 ldmx::Track* truth_trk = nullptr;
194
195 auto it =
196 std::find_if(truthTrackCollection_->begin(),
197 truthTrackCollection_->end(), [&](const ldmx::Track& tt) {
198 return tt.getTrackID() == track.getTrackID();
199 });
200
201 double trackTruthProb = track.getTruthProb();
202
203 if (it != truthTrackCollection_->end() && trackTruthProb >= trackProb_cut_)
204 truth_trk = &(*it);
205
206 // Match not found
207 if (!truth_trk) return;
208
209 double truth_phi = truth_trk->getPhi();
210 double truth_d0 = truth_trk->getD0();
211 double truth_z0 = truth_trk->getZ0();
212 double truth_theta = truth_trk->getTheta();
213 double truth_qop = truth_trk->getQoP();
214 double truth_p = 1. / abs(truth_trk->getQoP());
215 std::vector<double> truth_mom = truth_trk->getMomentum();
216
217 // Polar angle
218 // The momentum in the plane transverse wrt the beam axis
219 double truth_pt_beam =
220 std::sqrt(truth_mom[1] * truth_mom[1] + truth_mom[2] * truth_mom[2]);
221
222 double truth_beam_angle = std::atan2(truth_pt_beam, truth_mom[0]);
223
224 // Fill reco plots for efficiencies - numerator. The quantities are truth
225 histograms_.fill(title + "match_prob", trackTruthProb);
226 histograms_.fill(title + "match_d0", truth_d0);
227 histograms_.fill(title + "match_z0", truth_z0);
228 histograms_.fill(title + "match_phi", truth_phi);
229 histograms_.fill(title + "match_theta", truth_theta);
230 histograms_.fill(title + "match_p", truth_p);
231 histograms_.fill(title + "match_qop", truth_qop);
232 histograms_.fill(title + "match_beam_angle", truth_beam_angle);
233 histograms_.fill(title + "match_nHits", measurements.size());
234 for (auto trackHit : track.getMeasurementsIdxs()) {
235 histograms_.fill(title + "match_LayersHit",
236 measurements.at(trackHit).getLayer());
237 }
238
239 // For some particles
240
241 if (pidmap.count(truth_trk->getPdgID()) != 0) {
242 histograms_.fill(title + "match_PID", pidmap[truth_trk->getPdgID()]);
243
244 // TODO do this properly.
245
246 if (pidmap[truth_trk->getPdgID()] == PIDBins::kminus) {
247 histograms_.fill(title + "match_kminus_p", truth_p);
248 }
249
250 if (pidmap[truth_trk->getPdgID()] == PIDBins::kplus) {
251 histograms_.fill(title + "match_kplus_p", truth_p);
252 }
253
254 if (pidmap[truth_trk->getPdgID()] == PIDBins::piminus) {
255 histograms_.fill(title + "match_piminus_p", truth_p);
256 }
257
258 if (pidmap[truth_trk->getPdgID()] == PIDBins::piplus) {
259 histograms_.fill(title + "match_piplus_p", truth_p);
260 }
261
262 if (pidmap[truth_trk->getPdgID()] == PIDBins::electron) {
263 histograms_.fill(title + "match_electron_p", truth_p);
264 }
265
266 if (pidmap[truth_trk->getPdgID()] == PIDBins::positron) {
267 histograms_.fill(title + "match_positron_p", truth_p);
268 }
269
270 if (pidmap[truth_trk->getPdgID()] == PIDBins::proton) {
271 histograms_.fill(title + "match_proton_p", truth_p);
272 }
273 }
274 } // Loop on tracks
275
276} // Efficiency plots
277
278void TrackingRecoDQM::TrackMonitoring(
279 const std::vector<ldmx::Track>& tracks,
280 const std::vector<ldmx::Measurement>& measurements, const std::string title,
281 const bool& doDetail, const bool& doTruth) {
282 for (auto& track : tracks) {
283 // Perigee track parameters
284 double trk_d0 = track.getD0();
285 double trk_z0 = track.getZ0();
286 double trk_qop = track.getQoP();
287 double trk_theta = track.getTheta();
288 double trk_phi = track.getPhi();
289 double trk_p = 1. / abs(trk_qop);
290 for (auto trackHit : track.getMeasurementsIdxs()) {
291 histograms_.fill(title + "LayersHit",
292 measurements.at(trackHit).getLayer());
293 }
294
295 std::vector<double> trk_mom = track.getMomentum();
296
297 // The transverse momentum in the bending plane
298 double pt_bending =
299 std::sqrt(trk_mom[0] * trk_mom[0] + trk_mom[1] * trk_mom[1]);
300
301 // The momentum in the plane transverse wrt the beam axis
302 double pt_beam =
303 std::sqrt(trk_mom[1] * trk_mom[1] + trk_mom[2] * trk_mom[2]);
304
305 // Covariance matrix
306 Acts::BoundSquareMatrix cov =
307 tracking::sim::utils::unpackCov(track.getPerigeeCov());
308
309 double sigmad0 = sqrt(
310 cov(Acts::BoundIndices::eBoundLoc0, Acts::BoundIndices::eBoundLoc0));
311 double sigmaz0 = sqrt(
312 cov(Acts::BoundIndices::eBoundLoc1, Acts::BoundIndices::eBoundLoc1));
313 double sigmaphi =
314 sqrt(cov(Acts::BoundIndices::eBoundPhi, Acts::BoundIndices::eBoundPhi));
315 double sigmatheta = sqrt(
316 cov(Acts::BoundIndices::eBoundTheta, Acts::BoundIndices::eBoundTheta));
317 double sigmaqop = sqrt(cov(Acts::BoundIndices::eBoundQOverP,
318 Acts::BoundIndices::eBoundQOverP));
319 double sigmap = (1. / trk_qop) * (1. / trk_qop) * sigmaqop;
320
321 histograms_.fill(title + "d0", trk_d0);
322 histograms_.fill(title + "z0", trk_z0);
323 histograms_.fill(title + "qop", trk_qop);
324 histograms_.fill(title + "phi", trk_phi);
325 histograms_.fill(title + "theta", trk_theta);
326 histograms_.fill(title + "p", std::abs(1. / trk_qop));
327
328 if (doDetail) {
329 histograms_.fill(title + "px", trk_mom[0]);
330 histograms_.fill(title + "py", trk_mom[1]);
331 histograms_.fill(title + "pz", trk_mom[2]);
332
333 histograms_.fill(title + "pt_bending", pt_bending);
334 histograms_.fill(title + "pt_beam", pt_beam);
335
336 histograms_.fill(title + "nHits", track.getNhits());
337 histograms_.fill(title + "Chi2", track.getChi2());
338 histograms_.fill(title + "ndf", track.getNdf());
339 histograms_.fill(title + "Chi2_per_ndf",
340 track.getChi2() / track.getNdf());
341 histograms_.fill(title + "nShared", track.getNsharedHits());
342
343 histograms_.fill(title + "d0_err", sigmad0);
344 histograms_.fill(title + "z0_err", sigmaz0);
345 histograms_.fill(title + "phi_err", sigmaphi);
346 histograms_.fill(title + "theta_err", sigmatheta);
347 histograms_.fill(title + "qop_err", sigmaqop);
348 histograms_.fill(title + "p_err", sigmap);
349
350 // 2D Error plots
351
352 double p = std::abs(1. / trk_qop);
353 histograms_.fill(title + "d0_err_vs_p", p, sigmad0);
354 histograms_.fill(title + "z0_err_vs_p", std::abs(1. / trk_qop), sigmaz0);
355 histograms_.fill(title + "p_err_vs_p", std::abs(1. / trk_qop), sigmap);
356
357 if (track.getNhits() == 8)
358 histograms_.fill(title + "p_err_vs_p_8hits", p, sigmap);
359 else if (track.getNhits() == 9)
360 histograms_.fill(title + "p_err_vs_p_9hits", p, sigmap);
361 else if (track.getNhits() == 10)
362 histograms_.fill(title + "p_err_vs_p_10hits", p, sigmap);
363 }
364
365 if (doTruth) {
366 // Match to the truth track
367 ldmx::Track* truth_trk = nullptr;
368
369 auto it = std::find_if(truthTrackCollection_->begin(),
370 truthTrackCollection_->end(),
371 [&](const ldmx::Track& tt) {
372 return tt.getTrackID() == track.getTrackID();
373 });
374
375 double trackTruthProb = track.getTruthProb();
376
377 if (it != truthTrackCollection_->end() &&
378 trackTruthProb >= trackProb_cut_)
379 truth_trk = &(*it);
380
381 // Found matched track
382 if (truth_trk) {
383 double truth_d0 = truth_trk->getD0();
384 double truth_z0 = truth_trk->getZ0();
385 double truth_phi = truth_trk->getPhi();
386 double truth_theta = truth_trk->getTheta();
387 double truth_qop = truth_trk->getQoP();
388 double truth_p = 1. / abs(truth_trk->getQoP());
389 std::vector<double> truth_mom = truth_trk->getMomentum();
390 // Polar angle
391 // The momentum in the plane transverse wrt the beam axis
392 double truth_pt_beam = std::sqrt(truth_mom[1] * truth_mom[1] +
393 truth_mom[2] * truth_mom[2]);
394
395 // histograms_.fill(title+"truth_d0", truth_d0);
396 // histograms_.fill(title+"truth_z0", truth_z0);
397 // histograms_.fill(title+"truth_phi", truth_phi);
398 // histograms_.fill(title+"truth_theta",truth_theta);
399 // histograms_.fill(title+"truth_qop", truth_qop);
400 // histograms_.fill(title+"truth_p", truth_p);
401
402 double res_d0 = trk_d0 - truth_d0;
403 double res_z0 = trk_z0 - truth_z0;
404 double res_phi = trk_phi - truth_phi;
405 double res_theta = trk_theta - truth_theta;
406 double res_qop = trk_qop - truth_qop;
407 double res_p = trk_p - truth_p;
408 double res_pt_beam = pt_beam - truth_pt_beam;
409
410 histograms_.fill(title + "res_d0", res_d0);
411 histograms_.fill(title + "res_z0", res_z0);
412 histograms_.fill(title + "res_phi", res_phi);
413 histograms_.fill(title + "res_theta", res_theta);
414 histograms_.fill(title + "res_qop", res_qop);
415 histograms_.fill(title + "res_p", res_p);
416 histograms_.fill(title + "res_pt_beam", res_pt_beam);
417
418 double pull_d0 = res_d0 / sigmad0;
419 double pull_z0 = res_z0 / sigmaz0;
420 double pull_phi = res_phi / sigmaphi;
421 double pull_theta = res_theta / sigmatheta;
422 double pull_qop = res_qop / sigmaqop;
423 double pull_p = res_p / sigmap;
424
425 histograms_.fill(title + "pull_d0", pull_d0);
426 histograms_.fill(title + "pull_z0", pull_z0);
427 histograms_.fill(title + "pull_phi", pull_phi);
428 histograms_.fill(title + "pull_theta", pull_theta);
429 histograms_.fill(title + "pull_qop", pull_qop);
430 histograms_.fill(title + "pull_p", pull_p);
431
432 // Error plots from residuals
433
434 histograms_.fill(title + "res_p_vs_p", truth_p, res_p);
435
436 histograms_.fill(title + "res_qop_vs_p", truth_p, res_qop);
437 histograms_.fill(title + "res_d0_vs_p", truth_p, res_d0);
438 histograms_.fill(title + "res_z0_vs_p", truth_p, res_z0);
439 histograms_.fill(title + "res_phi_vs_p", truth_p, res_phi);
440 histograms_.fill(title + "res_theta_vs_p", truth_p, res_theta);
441
442 histograms_.fill(title + "pull_qop_vs_p", truth_p, pull_qop);
443 histograms_.fill(title + "pull_d0_vs_p", truth_p, pull_d0);
444 histograms_.fill(title + "pull_z0_vs_p", truth_p, pull_z0);
445 histograms_.fill(title + "pull_phi_vs_p", truth_p, pull_phi);
446 histograms_.fill(title + "pull_theta_vs_p", truth_p, pull_theta);
447
448 if (track.getNhits() == 8)
449 histograms_.fill(title + "res_p_vs_p_8hits", truth_p, res_p);
450 else if (track.getNhits() == 9)
451 histograms_.fill(title + "res_p_vs_p_9hits", truth_p, res_p);
452 else if (track.getNhits() == 10)
453 histograms_.fill(title + "res_p_vs_p_10hits", truth_p, res_p);
454
455 histograms_.fill(title + "res_pt_beam_vs_p", truth_pt_beam,
456 res_pt_beam);
457
458 } // found matched track
459 } // do TruthComparison
460 } // loop on tracks
461
462} // Track Monitoring
463
464void TrackingRecoDQM::TrackStateMonitoring(const ldmx::Tracks& tracks,
465 ldmx::TrackStateType ts_type,
466 const std::string& ts_title) {
467 for (auto& track : tracks) {
468 // Match the tracks to truth
469 ldmx::Track* truth_trk = nullptr;
470
471 auto it =
472 std::find_if(truthTrackCollection_->begin(),
473 truthTrackCollection_->end(), [&](const ldmx::Track& tt) {
474 return tt.getTrackID() == track.getTrackID();
475 });
476
477 double trackTruthProb = track.getTruthProb();
478
479 if (it != truthTrackCollection_->end() && trackTruthProb >= trackProb_cut_)
480 truth_trk = &(*it);
481
482 // Match not found, skip track
483 if (!truth_trk) continue;
484
485 // TruthTrack doesn't have the right amount of states
486
487 auto trk_ts = track.getTrackState(ts_type);
488 auto truth_ts = truth_trk->getTrackState(ts_type);
489
490 if (!trk_ts.has_value()) continue;
491
492 if (!truth_ts.has_value()) continue;
493
494 ldmx::Track::TrackState& truthTargetState = truth_ts.value();
495 ldmx::Track::TrackState& TargetState = trk_ts.value();
496
497 ldmx_log(debug) << "Unpacking covariance matrix" << std::endl;
498 Acts::BoundSquareMatrix cov =
499 tracking::sim::utils::unpackCov(TargetState.cov);
500
501 [[maybe_unused]] double sigmaloc0 = sqrt(
502 cov(Acts::BoundIndices::eBoundLoc0, Acts::BoundIndices::eBoundLoc0));
503 [[maybe_unused]] double sigmaloc1 = sqrt(
504 cov(Acts::BoundIndices::eBoundLoc1, Acts::BoundIndices::eBoundLoc1));
505 [[maybe_unused]] double sigmaphi =
506 sqrt(cov(Acts::BoundIndices::eBoundPhi, Acts::BoundIndices::eBoundPhi));
507 [[maybe_unused]] double sigmatheta = sqrt(
508 cov(Acts::BoundIndices::eBoundTheta, Acts::BoundIndices::eBoundTheta));
509 [[maybe_unused]] double sigmaqop = sqrt(cov(
510 Acts::BoundIndices::eBoundQOverP, Acts::BoundIndices::eBoundQOverP));
511
512 double trk_qop = track.getQoP();
513 double trk_p = 1. / abs(trk_qop);
514
515 double track_state_loc0 = TargetState.params[0];
516 double track_state_loc1 = TargetState.params[1];
517 [[maybe_unused]] double track_state_phi = TargetState.params[2];
518 [[maybe_unused]] double track_state_theta = TargetState.params[3];
519 [[maybe_unused]] double track_state_p = TargetState.params[4];
520
521 double truth_state_loc0 = truthTargetState.params[0];
522 double truth_state_loc1 = truthTargetState.params[1];
523 [[maybe_unused]] double truth_state_phi = truthTargetState.params[2];
524 [[maybe_unused]] double truth_state_theta = truthTargetState.params[3];
525 [[maybe_unused]] double truth_state_p = truthTargetState.params[4];
526
527 // Check that the track state is filled
528 if (TargetState.params.size() < 5) continue;
529
530 histograms_.fill(title_ + "trk_" + ts_title + "_loc0", track_state_loc0);
531 histograms_.fill(title_ + "trk_" + ts_title + "_loc1", track_state_loc1);
532 histograms_.fill(title_ + ts_title + "_sp_hit_X", truth_state_loc0);
533 histograms_.fill(title_ + ts_title + "_sp_hit_Y", truth_state_loc1);
534
535 // TH1F The difference(residual) between end_loc0 and sp_hit_X
536 histograms_.fill(title_ + "trk_" + ts_title + "_loc0-sp_hit_X",
537 track_state_loc0 - truth_state_loc0);
538 histograms_.fill(title_ + "trk_" + ts_title + "_loc1-sp_hit_Y",
539 track_state_loc1 - truth_state_loc1);
540
541 // TH1F The pulls of loc0 and loc1
542 histograms_.fill(title_ + ts_title + "_Pulls_of_loc0",
543 (track_state_loc0 - truth_state_loc0) / sigmaloc0);
544 histograms_.fill(title_ + ts_title + "_Pulls_of_loc1",
545 (track_state_loc1 - truth_state_loc1) / sigmaloc1);
546
547 // TODO:: TH1F The pulls of phi, theta, qop
548
549 // TH2F residual vs Nhits
550 histograms_.fill(title_ + ts_title + "_res_loc0-vs-N_hits",
551 track.getNhits(), track_state_loc0 - truth_state_loc0);
552 histograms_.fill(title_ + ts_title + "_res_loc1-vs-N_hits",
553 track.getNhits(), track_state_loc1 - truth_state_loc1);
554
555 // TH2F pulls vs Nhits
556 histograms_.fill(title_ + ts_title + "_pulls_loc0-vs-N_hits",
557 track.getNhits(),
558 (track_state_loc0 - truth_state_loc0) / sigmaloc0);
559 histograms_.fill(title_ + ts_title + "_pulls_loc1-vs-N_hits",
560 track.getNhits(),
561 (track_state_loc1 - truth_state_loc1) / sigmaloc1);
562
563 // TH2F residual vs trk_p
564 histograms_.fill(title_ + ts_title + "_res_loc0-vs-trk_p", trk_p,
565 track_state_loc0 - truth_state_loc0);
566 histograms_.fill(title_ + ts_title + "_res_loc1-vs-trk_p", trk_p,
567 track_state_loc1 - truth_state_loc1);
568
569 // TH2F pulls vs trk_p
570 histograms_.fill(title_ + ts_title + "_pulls_loc0-vs-trk_p", trk_p,
571 (track_state_loc0 - truth_state_loc0) / sigmaloc0);
572 histograms_.fill(title_ + ts_title + "_pulls_loc1-vs-trk_p", trk_p,
573 (track_state_loc1 - truth_state_loc1) / sigmaloc1);
574
575 } // loop on tracks
576}
577
578void TrackingRecoDQM::sortTracks(const std::vector<ldmx::Track>& tracks,
579 std::vector<ldmx::Track>& uniqueTracks,
580 std::vector<ldmx::Track>& duplicateTracks,
581 std::vector<ldmx::Track>& fakeTracks) {
582 // Create a copy of the const vector so we can sort it
583 std::vector<ldmx::Track> sortedTracks = tracks;
584
585 // Sort the vector of Track objects based on their trackID member
586 std::sort(sortedTracks.begin(), sortedTracks.end(),
587 [](ldmx::Track& t1, ldmx::Track& t2) {
588 return t1.getTrackID() < t2.getTrackID();
589 });
590
591 // Loop over the sorted vector of Track objects
592 for (size_t i = 0; i < sortedTracks.size(); i++) {
593 if (sortedTracks[i].getTruthProb() < trackProb_cut_)
594 fakeTracks.push_back(sortedTracks[i]);
595 else { // not a fake track
596 // If this is the first Track object with this trackID, add it to the
597 // uniqueTracks vector directly
598 if (uniqueTracks.size() == 0 ||
599 sortedTracks[i].getTrackID() != sortedTracks[i - 1].getTrackID()) {
600 uniqueTracks.push_back(sortedTracks[i]);
601 }
602 // Otherwise, add it to the duplicateTracks vector if its truthProb is
603 // lower than the existing Track object Otherwise, if the truthProbability
604 // is higher than the track stored in uniqueTracks, put it in uniqueTracks
605 // and move the uniqueTracks.back to duplicateTracks.
606 else if (sortedTracks[i].getTruthProb() >
607 uniqueTracks.back().getTruthProb()) {
608 duplicateTracks.push_back(uniqueTracks.back());
609 uniqueTracks.back() = sortedTracks[i];
610 }
611 // Otherwise, add it to the duplicateTracks vector
612 else {
613 duplicateTracks.push_back(sortedTracks[i]);
614 }
615 } // a real track
616 } // loop on sorted tracks
617 // The total number of elements in the uniqueTracks and duplicateTracks
618 // vectors should be equal to the number of elements in the original tracks
619 // vector
620 if (uniqueTracks.size() + duplicateTracks.size() + fakeTracks.size() !=
621 tracks.size()) {
622 std::cerr << "Error: unique and duplicate tracks vectors do not add up to "
623 "original tracks vector"
624 << std::endl;
625 return;
626 }
627
628 if (debug_) {
629 // Iterate through the uniqueTracks vector and duplicateTracks vector
630 std::cout << "Unique tracks:" << std::endl;
631 for (const ldmx::Track& track : uniqueTracks) {
632 std::cout << "Track ID: " << track.getTrackID()
633 << ", Truth Prob: " << track.getTruthProb() << std::endl;
634 }
635 std::cout << "Duplicate tracks:" << std::endl;
636 for (const ldmx::Track& track : duplicateTracks) {
637 std::cout << "Track ID: " << track.getTrackID()
638 << ", Truth Prob: " << track.getTruthProb() << std::endl;
639 }
640 std::cout << "Fake tracks:" << std::endl;
641 for (const ldmx::Track& track : fakeTracks) {
642 std::cout << "Track ID: " << track.getTrackID()
643 << ", Truth Prob: " << track.getTruthProb() << std::endl;
644 }
645 }
646}
647} // namespace tracking::dqm
648
649DECLARE_ANALYZER_NS(tracking::dqm, TrackingRecoDQM)
DECLARE_ANALYZER_NS
#define DECLARE_ANALYZER_NS(NS, CLASS)
Macro which allows the framework to construct an analyzer given its name during configuration.
Definition EventProcessor.h:391
framework::EventProcessor::histograms_
HistogramHelper histograms_
Interface class for making and filling histograms.
Definition EventProcessor.h:216
framework::Event
Implements an event buffer system for storing event data.
Definition Event.h:41
framework::Event::getCollection
const std::vector< ContentType > & getCollection(const std::string &collectionName, const std::string &passName="") const
Get a collection (std::vector) of objects from the event bus.
Definition Event.h:386
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:92
framework::HistogramHelper::fill
void fill(const std::string &name, const double &val)
Fill a 1D histogram.
Definition Histograms.h:166
framework::config::Parameters
Class encapsulating parameters for configuring a processor.
Definition Parameters.h:27
framework::config::Parameters::getParameter
T getParameter(const std::string &name) const
Retrieve the parameter of the given name.
Definition Parameters.h:89
ldmx::SimTrackerHit
Represents a simulated tracker hit in the simulation.
Definition SimTrackerHit.h:24
ldmx::Track
Implementation of a track object.
Definition Track.h:52
tracking::dqm::TrackingRecoDQM::TrackStateMonitoring
void TrackStateMonitoring(const ldmx::Tracks &tracks, ldmx::TrackStateType ts_type, const std::string &ts_title)
Monitoring plots for tracks extrapolated to the ECAL Scoring plane.
Definition TrackingRecoDQM.cxx:464
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:117
tracking::dqm::TrackingRecoDQM::configure
void configure(framework::config::Parameters &parameters) override
Configure the analyzer using the given user specified parameters.
Definition TrackingRecoDQM.cxx:10
tracking::dqm::TrackingRecoDQM::analyze
void analyze(const framework::Event &event) override
Process the event and make histograms or summaries.
Definition TrackingRecoDQM.cxx:36
Generated by doxygen 1.9.8