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