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