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 const std::vector<double> cov = track.getCov();
139 double sigma_phi = phiAngleError(track.getSlopeX(), cov);
140 double sigma_theta =
141 thetaAngleError(track.getSlopeX(), track.getSlopeY(), cov);
142 double sigma_loc0_target = std::sqrt(cov[4]);
143 double sigma_loc1_target = std::sqrt(cov[9]);
144 double sigma_loc0_ecal =
145 locError(cov.at(0), cov.at(4), cov.at(1), track.getEcalLayer1Z());
146 double sigma_loc1_ecal =
147 locError(cov.at(7), cov.at(9), cov.at(8), track.getEcalLayer1Z());
148
149 histograms_.fill(title + "phi", trk_phi);
150 histograms_.fill(title + "theta", trk_theta);
151 histograms_.fill(title_ + "trk_target_loc0", track_state_loc0_target);
152 histograms_.fill(title_ + "trk_target_loc1", track_state_loc1_target);
153 histograms_.fill(title_ + "trk_ecal_loc0", track_state_loc0_ecal);
154 histograms_.fill(title_ + "trk_ecal_loc1", track_state_loc1_ecal);
155 histograms_.fill(title + "trk_PID", track_pdg_id);
156
157 if (do_detail) {
158 histograms_.fill(title + "nHits", track.getNhits());
159 histograms_.fill(title + "Chi2", track.getChi2());
160 histograms_.fill(title + "ndf", track.getNdf());
161 histograms_.fill(title + "Chi2_per_ndf",
162 track.getChi2() / track.getNdf());
163 histograms_.fill(title + "dRecHit", track.getDistanceToRecHit());
164
165 histograms_.fill(title + "phi_err", sigma_phi);
166 histograms_.fill(title + "theta_err", sigma_theta);
167
168 if (do_truth) {
169 // Match to the truth track
170 // TODO: Currently, we assume the truth tracks have no uncertainty, but
171 // the parameters for the truth tracks
172 // TODO: are found from a fitting algorithm, and so have some
173 // uncertainty. Is there a better way to represent
174 // TODO: the truth tracks for a more accurate comparison?
175 ldmx::StraightTrack* truth_trk = nullptr;
176
177 // Only compare truth and reco tracks with same ID
178 auto it = std::find_if(truth_track_collection_->begin(),
179 truth_track_collection_->end(),
180 [&](const ldmx::StraightTrack& tt) {
181 return tt.getTrackID() == track.getTrackID();
182 });
183
184 double track_truth_prob = track.getTruthProb();
185
186 // make sure truth track has good enough truth prob.
187 if (it != truth_track_collection_->end() &&
188 track_truth_prob >= track_prob_cut_) {
189 truth_trk = &(*it);
190 }
191
192 // Found matched track
193 if (truth_trk) {
194 double truth_theta = truth_trk->getTheta();
195 double truth_phi = truth_trk->getPhi();
196 double truth_state_loc0_target = truth_trk->getTargetX();
197 double truth_state_loc1_target = truth_trk->getTargetY();
198 double truth_state_loc0_ecal = truth_trk->getEcalLayer1X();
199 double truth_state_loc1_ecal = truth_trk->getEcalLayer1Y();
200 int truth_pdg_id = truth_trk->getPdgID();
201
202 histograms_.fill(title + "truth_phi", truth_phi);
203 histograms_.fill(title + "truth_theta", truth_theta);
204 histograms_.fill(title + "truth_PID", truth_pdg_id);
205
206 double res_phi = trk_phi - truth_phi;
207 double res_theta = trk_theta - truth_theta;
208
209 histograms_.fill(title + "res_phi", res_phi);
210 histograms_.fill(title + "res_theta", res_theta);
211
212 double pull_phi = res_phi / sigma_phi;
213 double pull_theta = res_theta / sigma_theta;
214 histograms_.fill(title + "pull_phi", pull_phi);
215 histograms_.fill(title + "pull_theta", pull_theta);
216
217 // TH1F The difference(residual) between end_loc and truth_endloc
218 histograms_.fill(title_ + "trk_target_loc0-truth_target_loc0",
219 track_state_loc0_target - truth_state_loc0_target);
220 histograms_.fill(title_ + "trk_target_loc1-truth_target_loc1",
221 track_state_loc1_target - truth_state_loc1_target);
222 histograms_.fill(title_ + "trk_ecal_loc0-truth_ecal_loc0",
223 track_state_loc0_ecal - truth_state_loc0_ecal);
224 histograms_.fill(title_ + "trk_ecal_loc1-truth_ecal_loc1",
225 track_state_loc1_ecal - truth_state_loc1_ecal);
226
227 // TH1F The pulls of loc0 and loc1
228 histograms_.fill(title_ + "target_Pulls_of_loc0",
229 (track_state_loc0_target - truth_state_loc0_target) /
230 sigma_loc0_target);
231 histograms_.fill(title_ + "target_Pulls_of_loc1",
232 (track_state_loc1_target - truth_state_loc1_target) /
233 sigma_loc1_target);
234 histograms_.fill(title_ + "ecal_Pulls_of_loc0",
235 (track_state_loc0_ecal - truth_state_loc0_ecal) /
236 sigma_loc0_ecal);
237 histograms_.fill(title_ + "ecal_Pulls_of_loc1",
238 (track_state_loc1_ecal - truth_state_loc1_ecal) /
239 sigma_loc1_ecal);
240
241 // TH2F residual vs Nhits
242 histograms_.fill(title_ + "target_res_loc0-vs-N_hits",
243 track.getNhits(),
244 track_state_loc0_target - truth_state_loc0_target);
245 histograms_.fill(title_ + "target_res_loc1-vs-N_hits",
246 track.getNhits(),
247 track_state_loc1_target - truth_state_loc1_target);
248 histograms_.fill(title_ + "ecal_res_loc0-vs-N_hits", track.getNhits(),
249 track_state_loc0_ecal - truth_state_loc0_ecal);
250 histograms_.fill(title_ + "ecal_res_loc1-vs-N_hits", track.getNhits(),
251 track_state_loc1_ecal - truth_state_loc1_ecal);
252
253 // TH2F pulls vs Nhits
254 histograms_.fill(title_ + "target_pulls_loc0-vs-N_hits",
255 track.getNhits(),
256 (track_state_loc0_target - truth_state_loc0_target) /
257 sigma_loc0_target);
258 histograms_.fill(title_ + "target_pulls_loc1-vs-N_hits",
259 track.getNhits(),
260 (track_state_loc1_target - truth_state_loc1_target) /
261 sigma_loc1_target);
262 histograms_.fill(title_ + "ecal_pulls_loc0-vs-N_hits",
263 track.getNhits(),
264 (track_state_loc0_ecal - truth_state_loc0_ecal) /
265 sigma_loc0_ecal);
266 histograms_.fill(title_ + "ecal_pulls_loc1-vs-N_hits",
267 track.getNhits(),
268 (track_state_loc1_ecal - truth_state_loc1_ecal) /
269 sigma_loc1_ecal);
270
271 } // loop on tracks
272
273 } // do truth
274 } // do detail
275
276 } // for tracks
277} // TrackMonitoringUnique
278
279void StraightTracksDQM::sortTracks(
280 const std::vector<ldmx::StraightTrack>& tracks,
281 std::vector<ldmx::StraightTrack>& unique_tracks,
282 std::vector<ldmx::StraightTrack>& duplicate_tracks,
283 std::vector<ldmx::StraightTrack>& fake_tracks) {
284 std::vector<ldmx::StraightTrack> sorted_tracks = tracks;
285
286 // Sort the vector of Track objects based on their trackID member
287 std::sort(sorted_tracks.begin(), sorted_tracks.end(),
288 [](ldmx::StraightTrack& trk1, ldmx::StraightTrack& trk2) {
289 return trk1.getTrackID() < trk2.getTrackID();
290 });
291
292 // Loop over the sorted vector of Track objects
293 for (size_t i = 0; i < sorted_tracks.size(); i++) {
294 if (sorted_tracks[i].getTruthProb() < track_prob_cut_)
295 fake_tracks.push_back(sorted_tracks[i]);
296 else {
297 // If this is the first Track object with this trackID, add it to the
298 // uniqueTracks vector directly
299 if (unique_tracks.size() == 0 ||
300 sorted_tracks[i].getTrackID() != sorted_tracks[i - 1].getTrackID()) {
301 unique_tracks.push_back(sorted_tracks[i]);
302 }
303
304 // Otherwise, add it to the duplicateTracks vector if its truthProb is
305 // lower than the existing Track object Otherwise, if the truthProbability
306 // is higher than the track stored in uniqueTracks, put it in uniqueTracks
307 // and move the uniqueTracks.back to duplicateTracks.
308 else if (sorted_tracks[i].getTruthProb() >
309 unique_tracks.back().getTruthProb()) {
310 duplicate_tracks.push_back(unique_tracks.back());
311 unique_tracks.back() = sorted_tracks[i];
312 }
313
314 // Otherwise, add it to the duplicateTracks vector
315 else {
316 duplicate_tracks.push_back(sorted_tracks[i]);
317 }
318 } // else (a real track)
319 } // loop on sorted tracks
320
321 // The total number of elements in the uniqueTracks and duplicateTracks
322 // vectors should be equal to the number of elements in the original tracks
323 // vector
324 if ((unique_tracks.size() + duplicate_tracks.size() + fake_tracks.size()) !=
325 tracks.size()) {
326 ldmx_log(error) << "Unique and duplicate track vectors do not add up to "
327 "original tracks vector";
328 return;
329 } // if different tracks don't add up to correct total
330
331 // Iterate through the uniqueTracks vector and duplicateTracks vector
332 ldmx_log(trace) << "Unique tracks:";
333 for (const ldmx::StraightTrack& track : unique_tracks) {
334 ldmx_log(trace) << "Track ID: " << track.getTrackID()
335 << ", Truth Prob: " << track.getTruthProb();
336 }
337 ldmx_log(trace) << "Duplicate tracks:";
338 for (const ldmx::StraightTrack& track : duplicate_tracks) {
339 ldmx_log(trace) << "Track ID: " << track.getTrackID()
340 << ", Truth Prob: " << track.getTruthProb();
341 }
342 ldmx_log(trace) << "Fake tracks:";
343 for (const ldmx::StraightTrack& track : fake_tracks) {
344 ldmx_log(trace) << "Track ID: " << track.getTrackID()
345 << ", Truth Prob: " << track.getTruthProb();
346 }
347
348} // sortTracks
349
350double StraightTracksDQM::thetaAngleError(
351 double m_x, double m_y, const std::vector<double>& covariance_vector) {
352 double sqrt_term = std::sqrt(1 + (m_x * m_x));
353 double sum_term = (1 + (m_x * m_x) + (m_y * m_y));
354
355 double dtheta_dmx = (-m_x * m_y) / (sqrt_term * sum_term);
356 double dtheta_dmy = (sqrt_term / sum_term);
357
358 double sigma_mx2 = covariance_vector[0];
359 double sigma_my2 = covariance_vector[7];
360 double cov_mx_my = covariance_vector[2];
361
362 double sigma_theta2 = (dtheta_dmx * dtheta_dmx * sigma_mx2) +
363 (dtheta_dmy * dtheta_dmy * sigma_my2) +
364 (2 * dtheta_dmx * dtheta_dmy * cov_mx_my);
365
366 return std::sqrt(sigma_theta2);
367} // thetaAngleError
368
369double StraightTracksDQM::phiAngleError(
370 double m_x, const std::vector<double>& covariance_vector) {
371 double sum_term = (1 + (m_x * m_x));
372
373 double sigma_mx = std::sqrt(covariance_vector[0]);
374
375 double sigma_phi = (sigma_mx / sum_term);
376
377 return sigma_phi;
378} // phiAngleError
379
380double StraightTracksDQM::locError(double var_slope, double var_intercept,
381 double cov_slope_intercept, double z_pos) {
382 return std::sqrt((z_pos * z_pos * var_slope) + var_intercept +
383 (2 * z_pos * cov_slope_intercept));
384} // locAngleError
385
386} // namespace tracking::dqm
387
388DECLARE_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: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
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