LDMX Software
StraightTracksDQM.cxx
1#include "Tracking/dqm/StraightTracksDQM.h"
2
3#include <algorithm>
4#include <iostream>
5
6#include "Tracking/Sim/TrackingUtils.h"
7
8namespace tracking::dqm {
9
10void StraightTracksDQM::configure(framework::config::Parameters& parameters) {
11 track_collection_ =
12 parameters.get<std::string>("track_collection", "LinearRecoilTracks");
13 truth_collection_ = parameters.get<std::string>("truth_collection",
14 "LinearRecoilTruthTracks");
15 title_ = parameters.get<std::string>("title", "recoil_lin_trk_");
16 track_prob_cut_ = parameters.get<double>("trackProb_cut", 0.5);
17 subdetector_ = parameters.get<std::string>("subdetector", "Recoil");
18 measurement_collection_ = parameters.get<std::string>(
19 "measurement_collection", "DigiRecoilSimHits");
20
21 track_collection_events_passname_ =
22 parameters.get<std::string>("track_collection_events_passname");
23 truth_collection_events_passname_ =
24 parameters.get<std::string>("truth_collection_events_passname");
25 input_pass_name_ = parameters.get<std::string>("input_pass_name");
26
27 ldmx_log(info) << "Track Collection " << track_collection_;
28 ldmx_log(info) << "Truth Collection " << truth_collection_;
29} // configure
30
31void StraightTracksDQM::analyze(const framework::Event& event) {
32 ldmx_log(debug) << "DQM Reading in::" << track_collection_;
33
34 if (!event.exists(track_collection_, track_collection_events_passname_)) {
35 ldmx_log(error) << "trackCollection " << track_collection_
36 << " not in event";
37 return;
38 }
39
40 const std::vector<ldmx::StraightTrack> tracks =
41 event.getCollection<ldmx::StraightTrack>(track_collection_,
42 input_pass_name_);
43 const std::vector<ldmx::Measurement> measurements =
44 event.getCollection<ldmx::Measurement>(measurement_collection_,
45 input_pass_name_);
46
47 // Get the truth track collection
48 if (event.exists(truth_collection_, truth_collection_events_passname_)) {
49 truth_track_collection_ =
50 std::make_shared<std::vector<ldmx::StraightTrack>>(
51 event.getCollection<ldmx::StraightTrack>(truth_collection_,
52 input_pass_name_));
53 do_truth_comparison_ = true;
54 }
55
56 ldmx_log(debug) << "Do truth comparison::" << do_truth_comparison_;
57
58 if (do_truth_comparison_) {
59 sortTracks(tracks, unique_tracks_, duplicate_tracks_, fake_tracks_);
60 } else {
61 unique_tracks_ = tracks;
62 }
63
64 ldmx_log(debug) << "Filling histograms ";
65
66 // General Plots
67 histograms_.fill(title_ + "N_tracks", tracks.size());
68
69 ldmx_log(debug) << "Track Monitoring on Unique Tracks";
70
71 trackMonitoringUnique(unique_tracks_, measurements, title_, true, true);
72
73 ldmx_log(debug) << "Track Monitoring on duplicates and fakes";
74
75 // Fakes and duplicates
76 trackMonitoring(duplicate_tracks_, measurements, title_ + "dup_", false);
77 trackMonitoring(fake_tracks_, measurements, title_ + "fake_", false);
78
79 // Clear the vectors
80 unique_tracks_.clear();
81 duplicate_tracks_.clear();
82 fake_tracks_.clear();
83} // analyze
84
85void StraightTracksDQM::trackMonitoring(
86 const std::vector<ldmx::StraightTrack>& tracks,
87 const std::vector<ldmx::Measurement>& measurements, const std::string title,
88 const bool& do_detail) {
89 for (auto& track : tracks) {
90 double trk_theta = track.getTheta();
91 double trk_phi = track.getPhi();
92 double track_state_loc0_target = track.getTargetX();
93 double track_state_loc1_target = track.getTargetY();
94 double track_state_loc0_ecal = track.getEcalLayer1X();
95 double track_state_loc1_ecal = track.getEcalLayer1Y();
96 int track_pdg_id = track.getPdgID();
97
98 double sigma_phi = phiAngleError(track.getSlopeX(), track.getCov());
99 double sigma_theta =
100 thetaAngleError(track.getSlopeX(), track.getSlopeY(), track.getCov());
101
102 histograms_.fill(title + "phi", trk_phi);
103 histograms_.fill(title + "theta", trk_theta);
104 histograms_.fill(title_ + "trk_target_loc0", track_state_loc0_target);
105 histograms_.fill(title_ + "trk_target_loc1", track_state_loc1_target);
106 histograms_.fill(title_ + "trk_ecal_loc0", track_state_loc0_ecal);
107 histograms_.fill(title_ + "trk_ecal_loc1", track_state_loc1_ecal);
108 histograms_.fill(title + "trk_PID", track_pdg_id);
109
110 if (do_detail) {
111 histograms_.fill(title + "nHits", track.getNhits());
112 histograms_.fill(title + "Chi2", track.getChi2());
113 histograms_.fill(title + "ndf", track.getNdf());
114 histograms_.fill(title + "Chi2_per_ndf",
115 track.getChi2() / track.getNdf());
116 histograms_.fill(title + "dRecHit", track.getDistanceToRecHit());
117
118 histograms_.fill(title + "phi_err", sigma_phi);
119 histograms_.fill(title + "theta_err", sigma_theta);
120 } // do detail
121
122 } // for tracks
123} // TrackMonitoring
124
125void StraightTracksDQM::trackMonitoringUnique(
126 const std::vector<ldmx::StraightTrack>& tracks,
127 const std::vector<ldmx::Measurement>& measurements, const std::string title,
128 const bool& do_detail, const bool& do_truth) {
129 for (auto& track : tracks) {
130 double trk_theta = track.getTheta();
131 double trk_phi = track.getPhi();
132 double track_state_loc0_target = track.getTargetX();
133 double track_state_loc1_target = track.getTargetY();
134 double track_state_loc0_ecal = track.getEcalLayer1X();
135 double track_state_loc1_ecal = track.getEcalLayer1Y();
136 int track_pdg_id = track.getPdgID();
137
138 double sigma_phi = phiAngleError(track.getSlopeX(), track.getCov());
139 double sigma_theta =
140 thetaAngleError(track.getSlopeX(), track.getSlopeY(), track.getCov());
141 double sigma_loc0_target = std::sqrt(track.getCov()[4]);
142 double sigma_loc1_target = std::sqrt(track.getCov()[9]);
143 double sigma_loc0_ecal =
144 locError(track.getCov()[0], track.getCov()[4], track.getCov()[1],
145 track.getEcalLayer1Z());
146 double sigma_loc1_ecal =
147 locError(track.getCov()[7], track.getCov()[9], track.getCov()[8],
148 track.getEcalLayer1Z());
149
150 histograms_.fill(title + "phi", trk_phi);
151 histograms_.fill(title + "theta", trk_theta);
152 histograms_.fill(title_ + "trk_target_loc0", track_state_loc0_target);
153 histograms_.fill(title_ + "trk_target_loc1", track_state_loc1_target);
154 histograms_.fill(title_ + "trk_ecal_loc0", track_state_loc0_ecal);
155 histograms_.fill(title_ + "trk_ecal_loc1", track_state_loc1_ecal);
156 histograms_.fill(title + "trk_PID", track_pdg_id);
157
158 if (do_detail) {
159 histograms_.fill(title + "nHits", track.getNhits());
160 histograms_.fill(title + "Chi2", track.getChi2());
161 histograms_.fill(title + "ndf", track.getNdf());
162 histograms_.fill(title + "Chi2_per_ndf",
163 track.getChi2() / track.getNdf());
164 histograms_.fill(title + "dRecHit", track.getDistanceToRecHit());
165
166 histograms_.fill(title + "phi_err", sigma_phi);
167 histograms_.fill(title + "theta_err", sigma_theta);
168
169 if (do_truth) {
170 // Match to the truth track
171 // TODO: Currently, we assume the truth tracks have no uncertainty, but
172 // the parameters for the truth tracks
173 // TODO: are found from a fitting algorithm, and so have some
174 // uncertainty. Is there a better way to represent
175 // TODO: the truth tracks for a more accurate comparison?
176 ldmx::StraightTrack* truth_trk = nullptr;
177
178 // Only compare truth and reco tracks with same ID
179 auto it = std::find_if(truth_track_collection_->begin(),
180 truth_track_collection_->end(),
181 [&](const ldmx::StraightTrack& tt) {
182 return tt.getTrackID() == track.getTrackID();
183 });
184
185 double track_truth_prob = track.getTruthProb();
186
187 // make sure truth track has good enough truth prob.
188 if (it != truth_track_collection_->end() &&
189 track_truth_prob >= track_prob_cut_) {
190 truth_trk = &(*it);
191 }
192
193 // Found matched track
194 if (truth_trk) {
195 double truth_theta = truth_trk->getTheta();
196 double truth_phi = truth_trk->getPhi();
197 double truth_state_loc0_target = truth_trk->getTargetX();
198 double truth_state_loc1_target = truth_trk->getTargetY();
199 double truth_state_loc0_ecal = truth_trk->getEcalLayer1X();
200 double truth_state_loc1_ecal = truth_trk->getEcalLayer1Y();
201 int truth_pdg_id = truth_trk->getPdgID();
202
203 histograms_.fill(title + "truth_phi", truth_phi);
204 histograms_.fill(title + "truth_theta", truth_theta);
205 histograms_.fill(title + "truth_PID", truth_pdg_id);
206
207 double res_phi = trk_phi - truth_phi;
208 double res_theta = trk_theta - truth_theta;
209
210 histograms_.fill(title + "res_phi", res_phi);
211 histograms_.fill(title + "res_theta", res_theta);
212
213 double pull_phi = res_phi / sigma_phi;
214 double pull_theta = res_theta / sigma_theta;
215 histograms_.fill(title + "pull_phi", pull_phi);
216 histograms_.fill(title + "pull_theta", pull_theta);
217
218 // TH1F The difference(residual) between end_loc and truth_endloc
219 histograms_.fill(title_ + "trk_target_loc0-truth_target_loc0",
220 track_state_loc0_target - truth_state_loc0_target);
221 histograms_.fill(title_ + "trk_target_loc1-truth_target_loc1",
222 track_state_loc1_target - truth_state_loc1_target);
223 histograms_.fill(title_ + "trk_ecal_loc0-truth_ecal_loc0",
224 track_state_loc0_ecal - truth_state_loc0_ecal);
225 histograms_.fill(title_ + "trk_ecal_loc1-truth_ecal_loc1",
226 track_state_loc1_ecal - truth_state_loc1_ecal);
227
228 // TH1F The pulls of loc0 and loc1
229 histograms_.fill(title_ + "target_Pulls_of_loc0",
230 (track_state_loc0_target - truth_state_loc0_target) /
231 sigma_loc0_target);
232 histograms_.fill(title_ + "target_Pulls_of_loc1",
233 (track_state_loc1_target - truth_state_loc1_target) /
234 sigma_loc1_target);
235 histograms_.fill(title_ + "ecal_Pulls_of_loc0",
236 (track_state_loc0_ecal - truth_state_loc0_ecal) /
237 sigma_loc0_ecal);
238 histograms_.fill(title_ + "ecal_Pulls_of_loc1",
239 (track_state_loc1_ecal - truth_state_loc1_ecal) /
240 sigma_loc1_ecal);
241
242 // TH2F residual vs Nhits
243 histograms_.fill(title_ + "target_res_loc0-vs-N_hits",
244 track.getNhits(),
245 track_state_loc0_target - truth_state_loc0_target);
246 histograms_.fill(title_ + "target_res_loc1-vs-N_hits",
247 track.getNhits(),
248 track_state_loc1_target - truth_state_loc1_target);
249 histograms_.fill(title_ + "ecal_res_loc0-vs-N_hits", track.getNhits(),
250 track_state_loc0_ecal - truth_state_loc0_ecal);
251 histograms_.fill(title_ + "ecal_res_loc1-vs-N_hits", track.getNhits(),
252 track_state_loc1_ecal - truth_state_loc1_ecal);
253
254 // TH2F pulls vs Nhits
255 histograms_.fill(title_ + "target_pulls_loc0-vs-N_hits",
256 track.getNhits(),
257 (track_state_loc0_target - truth_state_loc0_target) /
258 sigma_loc0_target);
259 histograms_.fill(title_ + "target_pulls_loc1-vs-N_hits",
260 track.getNhits(),
261 (track_state_loc1_target - truth_state_loc1_target) /
262 sigma_loc1_target);
263 histograms_.fill(title_ + "ecal_pulls_loc0-vs-N_hits",
264 track.getNhits(),
265 (track_state_loc0_ecal - truth_state_loc0_ecal) /
266 sigma_loc0_ecal);
267 histograms_.fill(title_ + "ecal_pulls_loc1-vs-N_hits",
268 track.getNhits(),
269 (track_state_loc1_ecal - truth_state_loc1_ecal) /
270 sigma_loc1_ecal);
271
272 } // loop on tracks
273
274 } // do truth
275 } // do detail
276
277 } // for tracks
278} // TrackMonitoringUnique
279
280void StraightTracksDQM::sortTracks(
281 const std::vector<ldmx::StraightTrack>& tracks,
282 std::vector<ldmx::StraightTrack>& unique_tracks,
283 std::vector<ldmx::StraightTrack>& duplicate_tracks,
284 std::vector<ldmx::StraightTrack>& fake_tracks) {
285 std::vector<ldmx::StraightTrack> sorted_tracks = tracks;
286
287 // Sort the vector of Track objects based on their trackID member
288 std::sort(sorted_tracks.begin(), sorted_tracks.end(),
289 [](ldmx::StraightTrack& trk1, ldmx::StraightTrack& trk2) {
290 return trk1.getTrackID() < trk2.getTrackID();
291 });
292
293 // Loop over the sorted vector of Track objects
294 for (size_t i = 0; i < sorted_tracks.size(); i++) {
295 if (sorted_tracks[i].getTruthProb() < track_prob_cut_)
296 fake_tracks.push_back(sorted_tracks[i]);
297 else {
298 // If this is the first Track object with this trackID, add it to the
299 // uniqueTracks vector directly
300 if (unique_tracks.size() == 0 ||
301 sorted_tracks[i].getTrackID() != sorted_tracks[i - 1].getTrackID()) {
302 unique_tracks.push_back(sorted_tracks[i]);
303 }
304
305 // Otherwise, add it to the duplicateTracks vector if its truthProb is
306 // lower than the existing Track object Otherwise, if the truthProbability
307 // is higher than the track stored in uniqueTracks, put it in uniqueTracks
308 // and move the uniqueTracks.back to duplicateTracks.
309 else if (sorted_tracks[i].getTruthProb() >
310 unique_tracks.back().getTruthProb()) {
311 duplicate_tracks.push_back(unique_tracks.back());
312 unique_tracks.back() = sorted_tracks[i];
313 }
314
315 // Otherwise, add it to the duplicateTracks vector
316 else {
317 duplicate_tracks.push_back(sorted_tracks[i]);
318 }
319 } // else (a real track)
320 } // loop on sorted tracks
321
322 // The total number of elements in the uniqueTracks and duplicateTracks
323 // vectors should be equal to the number of elements in the original tracks
324 // vector
325 if ((unique_tracks.size() + duplicate_tracks.size() + fake_tracks.size()) !=
326 tracks.size()) {
327 ldmx_log(error) << "Unique and duplicate track vectors do not add up to "
328 "original tracks vector";
329 return;
330 } // if different tracks don't add up to correct total
331
332 // Iterate through the uniqueTracks vector and duplicateTracks vector
333 ldmx_log(trace) << "Unique tracks:";
334 for (const ldmx::StraightTrack& track : unique_tracks) {
335 ldmx_log(trace) << "Track ID: " << track.getTrackID()
336 << ", Truth Prob: " << track.getTruthProb();
337 }
338 ldmx_log(trace) << "Duplicate tracks:";
339 for (const ldmx::StraightTrack& track : duplicate_tracks) {
340 ldmx_log(trace) << "Track ID: " << track.getTrackID()
341 << ", Truth Prob: " << track.getTruthProb();
342 }
343 ldmx_log(trace) << "Fake tracks:";
344 for (const ldmx::StraightTrack& track : fake_tracks) {
345 ldmx_log(trace) << "Track ID: " << track.getTrackID()
346 << ", Truth Prob: " << track.getTruthProb();
347 }
348
349} // sortTracks
350
351double StraightTracksDQM::thetaAngleError(
352 double m_x, double m_y, const std::vector<double>& covariance_vector) {
353 double sqrt_term = std::sqrt(1 + (m_x * m_x));
354 double sum_term = (1 + (m_x * m_x) + (m_y * m_y));
355
356 double dtheta_dmx = (-m_x * m_y) / (sqrt_term * sum_term);
357 double dtheta_dmy = (sqrt_term / sum_term);
358
359 double sigma_mx2 = covariance_vector[0];
360 double sigma_my2 = covariance_vector[7];
361 double cov_mx_my = covariance_vector[2];
362
363 double sigma_theta2 = (dtheta_dmx * dtheta_dmx * sigma_mx2) +
364 (dtheta_dmy * dtheta_dmy * sigma_my2) +
365 (2 * dtheta_dmx * dtheta_dmy * cov_mx_my);
366
367 return std::sqrt(sigma_theta2);
368} // thetaAngleError
369
370double StraightTracksDQM::phiAngleError(
371 double m_x, const std::vector<double>& covariance_vector) {
372 double sum_term = (1 + (m_x * m_x));
373
374 double sigma_mx = std::sqrt(covariance_vector[0]);
375
376 double sigma_phi = (sigma_mx / sum_term);
377
378 return sigma_phi;
379} // phiAngleError
380
381double StraightTracksDQM::locError(double var_slope, double var_intercept,
382 double cov_slope_intercept, double z_pos) {
383 return std::sqrt((z_pos * z_pos * var_slope) + var_intercept +
384 (2 * z_pos * cov_slope_intercept));
385} // locAngleError
386
387} // namespace tracking::dqm
388
389DECLARE_ANALYZER(tracking::dqm::StraightTracksDQM)
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:92
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
tracking::dqm::StraightTracksDQM::configure
void configure(framework::config::Parameters &parameters) override
Callback for the EventProcessor to configure itself from the given set of parameters.
Definition StraightTracksDQM.cxx:10
tracking::dqm::StraightTracksDQM::analyze
void analyze(const framework::Event &event) override
Process the event and make histograms or summaries.
Definition StraightTracksDQM.cxx:31
Generated by doxygen 1.12.0