LDMX Software
TrackComparisonProcessor.cxx
1#include "Tracking/Reco/TrackComparisonProcessor.h"
2
3#include <cmath>
4#include <map>
5#include <vector>
6
7#include "Framework/Logger.h"
8#include "SimCore/Event/SimParticle.h"
9#include "Tracking/Event/Track.h"
10
11namespace tracking::reco {
12
13TrackComparisonProcessor::TrackComparisonProcessor(const std::string& name,
14 framework::Process& process)
15 : framework::Analyzer(name, process) {}
16
17void TrackComparisonProcessor::configure(
18 framework::config::Parameters& parameters) {
19 trk_collection_smear_ =
20 parameters.get<std::string>("trk_collection_smear", "TaggerTracks");
21 trk_collection_digi_ =
22 parameters.get<std::string>("trk_collection_digi", "TaggerDigiTracks");
23 pass_name_smear_ = parameters.get<std::string>("pass_name_smear", "");
24 pass_name_digi_ = parameters.get<std::string>("pass_name_digi", "");
25 do_tagger_ = parameters.get<bool>("do_tagger", true);
26 do_recoil_ = parameters.get<bool>("do_recoil", false);
27 recoil_collection_smear_ =
28 parameters.get<std::string>("recoil_collection_smear", "RecoilTracks");
29 recoil_collection_digi_ =
30 parameters.get<std::string>("recoil_collection_digi", "RecoilDigiTracks");
31 recoil_pass_smear_ = parameters.get<std::string>("recoil_pass_smear", "");
32 recoil_pass_digi_ = parameters.get<std::string>("recoil_pass_digi", "");
33 min_truth_prob_ = parameters.get<double>("min_truth_prob", 0.5);
34 sim_particles_pass_ = parameters.get<std::string>("sim_particles_pass", "");
35 output_file_ =
36 parameters.get<std::string>("output_file", "track_comparison.root");
37}
38
39void TrackComparisonProcessor::setupTree(TTree* tree, PairVars& v) {
40 tree->Branch("track_id", &v.track_id);
41 tree->Branch("truth_prob_s", &v.truth_prob_s);
42 tree->Branch("truth_prob_d", &v.truth_prob_d);
43 tree->Branch("nhits_s", &v.nhits_s);
44 tree->Branch("nhits_d", &v.nhits_d);
45 tree->Branch("chi2ndf_s", &v.chi2ndf_s);
46 tree->Branch("chi2ndf_d", &v.chi2ndf_d);
47 tree->Branch("d0_s", &v.d0_s);
48 tree->Branch("d0_d", &v.d0_d);
49 tree->Branch("z0_s", &v.z0_s);
50 tree->Branch("z0_d", &v.z0_d);
51 tree->Branch("phi_s", &v.phi_s);
52 tree->Branch("phi_d", &v.phi_d);
53 tree->Branch("theta_s", &v.theta_s);
54 tree->Branch("theta_d", &v.theta_d);
55 tree->Branch("qop_s", &v.qop_s);
56 tree->Branch("qop_d", &v.qop_d);
57 tree->Branch("p_s", &v.p_s);
58 tree->Branch("p_d", &v.p_d);
59 tree->Branch("delta_d0", &v.delta_d0);
60 tree->Branch("delta_z0", &v.delta_z0);
61 tree->Branch("delta_phi", &v.delta_phi);
62 tree->Branch("delta_theta", &v.delta_theta);
63 tree->Branch("delta_p_over_p", &v.delta_p_over_p);
64 tree->Branch("px_s", &v.px_s);
65 tree->Branch("py_s", &v.py_s);
66 tree->Branch("pz_s", &v.pz_s);
67 tree->Branch("px_d", &v.px_d);
68 tree->Branch("py_d", &v.py_d);
69 tree->Branch("pz_d", &v.pz_d);
70 tree->Branch("px_t", &v.px_t);
71 tree->Branch("py_t", &v.py_t);
72 tree->Branch("pz_t", &v.pz_t);
73 tree->Branch("p_t", &v.p_t);
74 tree->Branch("vx_t", &v.vx_t);
75 tree->Branch("vy_t", &v.vy_t);
76 tree->Branch("vz_t", &v.vz_t);
77 tree->Branch("delta_p_over_p_s", &v.delta_p_over_p_s);
78 tree->Branch("delta_p_over_p_d", &v.delta_p_over_p_d);
79}
80
81void TrackComparisonProcessor::onProcessStart() {
82 file_ = new TFile(output_file_.c_str(), "RECREATE");
83 file_->cd();
84
85 if (do_tagger_) {
86 tagger_tree_ =
87 new TTree("tagger_pairs", "Tagger smear-vs-digi track pairs");
88 tagger_tree_->SetDirectory(file_);
89 setupTree(tagger_tree_, tagger_vars_);
90
91 histograms_.create("tagger_delta_d0", "#Delta d_{0} (digi-smear) [mm]", 200,
92 -0.5, 0.5);
93 histograms_.create("tagger_delta_z0", "#Delta z_{0} (digi-smear) [mm]", 200,
94 -2.0, 2.0);
95 histograms_.create("tagger_delta_phi", "#Delta #phi (digi-smear) [rad]",
96 200, -0.02, 0.02);
97 histograms_.create("tagger_delta_theta", "#Delta #theta (digi-smear) [rad]",
98 200, -0.02, 0.02);
99 histograms_.create("tagger_delta_p_over_p", "#Delta p/p (digi-smear)/smear",
100 200, -0.1, 0.1);
101 histograms_.create("tagger_nhits_s", "N hits (smear)", 15, 0, 15);
102 histograms_.create("tagger_nhits_d", "N hits (digi)", 15, 0, 15);
103 histograms_.create("tagger_chi2ndf_s", "#chi^{2}/ndf (smear)", 100, 0, 10);
104 histograms_.create("tagger_chi2ndf_d", "#chi^{2}/ndf (digi)", 100, 0, 10);
105 histograms_.create("tagger_p_s", "p (smear) [MeV]", 200, 0, 8000);
106 histograms_.create("tagger_p_d", "p (digi) [MeV]", 200, 0, 8000);
107 histograms_.create("tagger_p_t", "p (truth) [MeV]", 200, 0, 8000);
108 histograms_.create("tagger_delta_p_over_p_s",
109 "#Delta p/p (smear-truth)/truth", 200, -0.1, 0.1);
110 histograms_.create("tagger_delta_p_over_p_d",
111 "#Delta p/p (digi-truth)/truth", 200, -0.1, 0.1);
112 }
113
114 if (do_recoil_) {
115 recoil_tree_ =
116 new TTree("recoil_pairs", "Recoil smear-vs-digi track pairs");
117 recoil_tree_->SetDirectory(file_);
118 setupTree(recoil_tree_, recoil_vars_);
119
120 histograms_.create("recoil_delta_d0", "#Delta d_{0} (digi-smear) [mm]", 200,
121 -2.0, 2.0);
122 histograms_.create("recoil_delta_z0", "#Delta z_{0} (digi-smear) [mm]", 200,
123 -5.0, 5.0);
124 histograms_.create("recoil_delta_phi", "#Delta #phi (digi-smear) [rad]",
125 200, -0.02, 0.02);
126 histograms_.create("recoil_delta_theta", "#Delta #theta (digi-smear) [rad]",
127 200, -0.02, 0.02);
128 histograms_.create("recoil_delta_p_over_p", "#Delta p/p (digi-smear)/smear",
129 200, -0.4, 0.4);
130 histograms_.create("recoil_nhits_s", "N hits (smear)", 15, 0, 15);
131 histograms_.create("recoil_nhits_d", "N hits (digi)", 15, 0, 15);
132 histograms_.create("recoil_chi2ndf_s", "#chi^{2}/ndf (smear)", 100, 0, 10);
133 histograms_.create("recoil_chi2ndf_d", "#chi^{2}/ndf (digi)", 100, 0, 10);
134 histograms_.create("recoil_p_s", "p (smear) [MeV]", 200, 0, 8000);
135 histograms_.create("recoil_p_d", "p (digi) [MeV]", 200, 0, 8000);
136 histograms_.create("recoil_p_t", "p (truth) [MeV]", 200, 0, 8000);
137 histograms_.create("recoil_delta_p_over_p_s",
138 "#Delta p/p (smear-truth)/truth", 200, -0.4, 0.4);
139 histograms_.create("recoil_delta_p_over_p_d",
140 "#Delta p/p (digi-truth)/truth", 200, -0.4, 0.4);
141 }
142}
143
144void TrackComparisonProcessor::fillPair(const ldmx::Track& smear,
145 const ldmx::Track& digi,
146 const ldmx::SimParticle& truth,
147 PairVars& v,
148 const std::string& prefix) {
149 v.track_id = smear.getTrackID();
150 v.truth_prob_s = smear.getTruthProb();
151 v.truth_prob_d = digi.getTruthProb();
152 v.nhits_s = smear.getNhits();
153 v.nhits_d = digi.getNhits();
154 v.chi2ndf_s = (smear.getNdf() > 0) ? smear.getChi2() / smear.getNdf() : -1;
155 v.chi2ndf_d = (digi.getNdf() > 0) ? digi.getChi2() / digi.getNdf() : -1;
156 v.d0_s = smear.getD0();
157 v.d0_d = digi.getD0();
158 v.z0_s = smear.getZ0();
159 v.z0_d = digi.getZ0();
160 v.phi_s = smear.getPhi();
161 v.phi_d = digi.getPhi();
162 v.theta_s = smear.getTheta();
163 v.theta_d = digi.getTheta();
164 v.qop_s = smear.getQoP();
165 v.qop_d = digi.getQoP();
166 // QoP is in e/GeV (ACTS units); convert to MeV so p_s/p_d match SimParticle
167 // units.
168 v.p_s = (v.qop_s != 0) ? std::abs(1000.0 / v.qop_s) : -1;
169 v.p_d = (v.qop_d != 0) ? std::abs(1000.0 / v.qop_d) : -1;
170 v.delta_d0 = v.d0_d - v.d0_s;
171 v.delta_z0 = v.z0_d - v.z0_s;
172 v.delta_phi = v.phi_d - v.phi_s;
173 v.delta_theta = v.theta_d - v.theta_s;
174 v.delta_p_over_p = (v.p_s > 0) ? (v.p_d - v.p_s) / v.p_s : -999;
175
176 auto mom_s = smear.getMomentumAtTarget();
177 if (mom_s.size() == 3) {
178 v.px_s = mom_s[0];
179 v.py_s = mom_s[1];
180 v.pz_s = mom_s[2];
181 }
182 auto mom_d = digi.getMomentumAtTarget();
183 if (mom_d.size() == 3) {
184 v.px_d = mom_d[0];
185 v.py_d = mom_d[1];
186 v.pz_d = mom_d[2];
187 }
188
189 auto mom_t = truth.getMomentum();
190 v.px_t = mom_t[0];
191 v.py_t = mom_t[1];
192 v.pz_t = mom_t[2];
193 v.p_t = std::sqrt(v.px_t * v.px_t + v.py_t * v.py_t + v.pz_t * v.pz_t);
194 auto vtx = truth.getVertex();
195 v.vx_t = vtx[0];
196 v.vy_t = vtx[1];
197 v.vz_t = vtx[2];
198
199 v.delta_p_over_p_s = (v.p_t > 0) ? (v.p_s - v.p_t) / v.p_t : -999;
200 v.delta_p_over_p_d = (v.p_t > 0) ? (v.p_d - v.p_t) / v.p_t : -999;
201
202 histograms_.fill(prefix + "delta_d0", v.delta_d0);
203 histograms_.fill(prefix + "delta_z0", v.delta_z0);
204 histograms_.fill(prefix + "delta_phi", v.delta_phi);
205 histograms_.fill(prefix + "delta_theta", v.delta_theta);
206 histograms_.fill(prefix + "delta_p_over_p", v.delta_p_over_p);
207 histograms_.fill(prefix + "nhits_s", v.nhits_s);
208 histograms_.fill(prefix + "nhits_d", v.nhits_d);
209 histograms_.fill(prefix + "chi2ndf_s", v.chi2ndf_s);
210 histograms_.fill(prefix + "chi2ndf_d", v.chi2ndf_d);
211 histograms_.fill(prefix + "p_s", v.p_s);
212 histograms_.fill(prefix + "p_d", v.p_d);
213 histograms_.fill(prefix + "p_t", v.p_t);
214 histograms_.fill(prefix + "delta_p_over_p_s", v.delta_p_over_p_s);
215 histograms_.fill(prefix + "delta_p_over_p_d", v.delta_p_over_p_d);
216}
217
218void TrackComparisonProcessor::processTracker(const framework::Event& event,
219 const std::string& coll_smear,
220 const std::string& pass_smear,
221 const std::string& coll_digi,
222 const std::string& pass_digi,
223 TTree* tree, PairVars& vars,
224 const std::string& histo_prefix) {
225 if (!event.exists(coll_smear, pass_smear)) {
226 ldmx_log(warn) << "Smear collection " << coll_smear << " not found";
227 return;
228 }
229 if (!event.exists(coll_digi, pass_digi)) {
230 ldmx_log(warn) << "Digi collection " << coll_digi << " not found";
231 return;
232 }
233
234 const auto& tracks_smear =
235 event.getCollection<ldmx::Track>(coll_smear, pass_smear);
236 const auto& tracks_digi =
237 event.getCollection<ldmx::Track>(coll_digi, pass_digi);
238
239 const auto& particle_map =
240 event.getMap<int, ldmx::SimParticle>("SimParticles", sim_particles_pass_);
241
242 std::map<int, const ldmx::Track*> smear_by_id;
243 for (const auto& t : tracks_smear) {
244 if (t.getTruthProb() >= min_truth_prob_) smear_by_id[t.getTrackID()] = &t;
245 }
246
247 for (const auto& t : tracks_digi) {
248 if (t.getTruthProb() < min_truth_prob_) continue;
249 auto it = smear_by_id.find(t.getTrackID());
250 if (it == smear_by_id.end()) continue;
251 auto pit = particle_map.find(t.getTrackID());
252 if (pit == particle_map.end()) continue;
253
254 fillPair(*it->second, t, pit->second, vars, histo_prefix);
255 tree->Fill();
256 }
257}
258
259void TrackComparisonProcessor::analyze(const framework::Event& event) {
260 if (do_tagger_) {
261 processTracker(event, trk_collection_smear_, pass_name_smear_,
262 trk_collection_digi_, pass_name_digi_, tagger_tree_,
263 tagger_vars_, "tagger_");
264 }
265 if (do_recoil_) {
266 processTracker(event, recoil_collection_smear_, recoil_pass_smear_,
267 recoil_collection_digi_, recoil_pass_digi_, recoil_tree_,
268 recoil_vars_, "recoil_");
269 }
270}
271
272void TrackComparisonProcessor::onProcessEnd() {
273 file_->cd();
274 if (tagger_tree_) tagger_tree_->Write("", TObject::kOverwrite);
275 if (recoil_tree_) recoil_tree_->Write("", TObject::kOverwrite);
276 file_->Close();
277}
278
279} // namespace tracking::reco
280
281DECLARE_ANALYZER(tracking::reco::TrackComparisonProcessor)
DECLARE_ANALYZER
#define DECLARE_ANALYZER(CLASS)
Macro which allows the framework to construct an analyzer given its name during configuration.
Definition EventProcessor.h:361
framework::Event
Implements an event buffer system for storing event data.
Definition Event.h:40
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:107
framework::Process
Class which represents the process under execution.
Definition Process.h:34
framework::config::Parameters
Class encapsulating parameters for configuring a processor.
Definition Parameters.h:26
framework::config::Parameters::get
const T & get(const std::string &name) const
Retrieve the parameter of the given name.
Definition Parameters.h:75
ldmx::SimParticle
Class representing a simulated particle.
Definition SimParticle.h:25
ldmx::SimParticle::getVertex
std::vector< double > getVertex() const
Get a vector containing the vertex of this particle in mm.
Definition SimParticle.h:115
ldmx::SimParticle::getMomentum
std::vector< double > getMomentum() const
Get a vector containing the momentum of this particle [MeV].
Definition SimParticle.h:150
ldmx::Track
Implementation of a track object.
Definition Track.h:54
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:220
tracking::reco::TrackComparisonProcessor
Compares tracking performance between a truth-smeared hit chain and a charge-digitized hit chain on a...
Definition TrackComparisonProcessor.h:26
framework
All classes in the ldmx-sw project use this namespace.
Generated by doxygen 1.12.0