EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
RootTrajectoryWriter.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file RootTrajectoryWriter.cpp
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2019-2020 CERN for the benefit of the Acts project
4 //
5 // This Source Code Form is subject to the terms of the Mozilla Public
6 // License, v. 2.0. If a copy of the MPL was not distributed with this
7 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
8 
10 
19 
20 #include <ios>
21 #include <stdexcept>
22 
23 #include <TFile.h>
24 #include <TTree.h>
25 
30 using Measurement =
33 
37  : WriterT(cfg.inputTrajectories, "RootTrajectoryWriter", lvl),
38  m_cfg(cfg),
39  m_outputFile(cfg.rootFile) {
40  // An input collection name and tree name must be specified
41  if (m_cfg.inputTrajectories.empty()) {
42  throw std::invalid_argument("Missing input trajectory collection");
43  } else if (m_cfg.inputParticles.empty()) {
44  throw std::invalid_argument("Missing input particle collection");
45  } else if (cfg.outputFilename.empty()) {
46  throw std::invalid_argument("Missing output filename");
47  } else if (m_cfg.outputTreename.empty()) {
48  throw std::invalid_argument("Missing tree name");
49  }
50 
51  // Setup ROOT I/O
52  if (m_outputFile == nullptr) {
54  m_outputFile = TFile::Open(path.c_str(), m_cfg.fileMode.c_str());
55  if (m_outputFile == nullptr) {
56  throw std::ios_base::failure("Could not open '" + path);
57  }
58  }
59  m_outputFile->cd();
60  m_outputTree =
61  new TTree(m_cfg.outputTreename.c_str(), m_cfg.outputTreename.c_str());
62  if (m_outputTree == nullptr)
63  throw std::bad_alloc();
64  else {
65  // I/O parameters
66  m_outputTree->Branch("event_nr", &m_eventNr);
67  m_outputTree->Branch("traj_nr", &m_trajNr);
68  m_outputTree->Branch("t_barcode", &m_t_barcode, "t_barcode/l");
69  m_outputTree->Branch("t_charge", &m_t_charge);
70  m_outputTree->Branch("t_time", &m_t_time);
71  m_outputTree->Branch("t_vx", &m_t_vx);
72  m_outputTree->Branch("t_vy", &m_t_vy);
73  m_outputTree->Branch("t_vz", &m_t_vz);
74  m_outputTree->Branch("t_px", &m_t_px);
75  m_outputTree->Branch("t_py", &m_t_py);
76  m_outputTree->Branch("t_pz", &m_t_pz);
77  m_outputTree->Branch("t_theta", &m_t_theta);
78  m_outputTree->Branch("t_phi", &m_t_phi);
79  m_outputTree->Branch("t_eta", &m_t_eta);
80  m_outputTree->Branch("t_pT", &m_t_pT);
81 
82  m_outputTree->Branch("t_x", &m_t_x);
83  m_outputTree->Branch("t_y", &m_t_y);
84  m_outputTree->Branch("t_z", &m_t_z);
85  m_outputTree->Branch("t_r", &m_t_r);
86  m_outputTree->Branch("t_dx", &m_t_dx);
87  m_outputTree->Branch("t_dy", &m_t_dy);
88  m_outputTree->Branch("t_dz", &m_t_dz);
89  m_outputTree->Branch("t_eLOC0", &m_t_eLOC0);
90  m_outputTree->Branch("t_eLOC1", &m_t_eLOC1);
91  m_outputTree->Branch("t_ePHI", &m_t_ePHI);
92  m_outputTree->Branch("t_eTHETA", &m_t_eTHETA);
93  m_outputTree->Branch("t_eQOP", &m_t_eQOP);
94  m_outputTree->Branch("t_eT", &m_t_eT);
95 
96  m_outputTree->Branch("nStates", &m_nStates);
97  m_outputTree->Branch("nMeasurements", &m_nMeasurements);
98  m_outputTree->Branch("volume_id", &m_volumeID);
99  m_outputTree->Branch("layer_id", &m_layerID);
100  m_outputTree->Branch("module_id", &m_moduleID);
101  m_outputTree->Branch("l_x_hit", &m_lx_hit);
102  m_outputTree->Branch("l_y_hit", &m_ly_hit);
103  m_outputTree->Branch("g_x_hit", &m_x_hit);
104  m_outputTree->Branch("g_y_hit", &m_y_hit);
105  m_outputTree->Branch("g_z_hit", &m_z_hit);
106  m_outputTree->Branch("res_x_hit", &m_res_x_hit);
107  m_outputTree->Branch("res_y_hit", &m_res_y_hit);
108  m_outputTree->Branch("err_x_hit", &m_err_x_hit);
109  m_outputTree->Branch("err_y_hit", &m_err_y_hit);
110  m_outputTree->Branch("pull_x_hit", &m_pull_x_hit);
111  m_outputTree->Branch("pull_y_hit", &m_pull_y_hit);
112  m_outputTree->Branch("dim_hit", &m_dim_hit);
113 
114  m_outputTree->Branch("hasFittedParams", &m_hasFittedParams);
115  m_outputTree->Branch("eLOC0_fit", &m_eLOC0_fit);
116  m_outputTree->Branch("eLOC1_fit", &m_eLOC1_fit);
117  m_outputTree->Branch("ePHI_fit", &m_ePHI_fit);
118  m_outputTree->Branch("eTHETA_fit", &m_eTHETA_fit);
119  m_outputTree->Branch("eQOP_fit", &m_eQOP_fit);
120  m_outputTree->Branch("eT_fit", &m_eT_fit);
121  m_outputTree->Branch("err_eLOC0_fit", &m_err_eLOC0_fit);
122  m_outputTree->Branch("err_eLOC1_fit", &m_err_eLOC1_fit);
123  m_outputTree->Branch("err_ePHI_fit", &m_err_ePHI_fit);
124  m_outputTree->Branch("err_eTHETA_fit", &m_err_eTHETA_fit);
125  m_outputTree->Branch("err_eQOP_fit", &m_err_eQOP_fit);
126  m_outputTree->Branch("err_eT_fit", &m_err_eT_fit);
127 
128  m_outputTree->Branch("nPredicted", &m_nPredicted);
129  m_outputTree->Branch("predicted", &m_prt);
130  m_outputTree->Branch("eLOC0_prt", &m_eLOC0_prt);
131  m_outputTree->Branch("eLOC1_prt", &m_eLOC1_prt);
132  m_outputTree->Branch("ePHI_prt", &m_ePHI_prt);
133  m_outputTree->Branch("eTHETA_prt", &m_eTHETA_prt);
134  m_outputTree->Branch("eQOP_prt", &m_eQOP_prt);
135  m_outputTree->Branch("eT_prt", &m_eT_prt);
136  m_outputTree->Branch("res_eLOC0_prt", &m_res_eLOC0_prt);
137  m_outputTree->Branch("res_eLOC1_prt", &m_res_eLOC1_prt);
138  m_outputTree->Branch("res_ePHI_prt", &m_res_ePHI_prt);
139  m_outputTree->Branch("res_eTHETA_prt", &m_res_eTHETA_prt);
140  m_outputTree->Branch("res_eQOP_prt", &m_res_eQOP_prt);
141  m_outputTree->Branch("res_eT_prt", &m_res_eT_prt);
142  m_outputTree->Branch("err_eLOC0_prt", &m_err_eLOC0_prt);
143  m_outputTree->Branch("err_eLOC1_prt", &m_err_eLOC1_prt);
144  m_outputTree->Branch("err_ePHI_prt", &m_err_ePHI_prt);
145  m_outputTree->Branch("err_eTHETA_prt", &m_err_eTHETA_prt);
146  m_outputTree->Branch("err_eQOP_prt", &m_err_eQOP_prt);
147  m_outputTree->Branch("err_eT_prt", &m_err_eT_prt);
148  m_outputTree->Branch("pull_eLOC0_prt", &m_pull_eLOC0_prt);
149  m_outputTree->Branch("pull_eLOC1_prt", &m_pull_eLOC1_prt);
150  m_outputTree->Branch("pull_ePHI_prt", &m_pull_ePHI_prt);
151  m_outputTree->Branch("pull_eTHETA_prt", &m_pull_eTHETA_prt);
152  m_outputTree->Branch("pull_eQOP_prt", &m_pull_eQOP_prt);
153  m_outputTree->Branch("pull_eT_prt", &m_pull_eT_prt);
154  m_outputTree->Branch("g_x_prt", &m_x_prt);
155  m_outputTree->Branch("g_y_prt", &m_y_prt);
156  m_outputTree->Branch("g_z_prt", &m_z_prt);
157  m_outputTree->Branch("px_prt", &m_px_prt);
158  m_outputTree->Branch("py_prt", &m_py_prt);
159  m_outputTree->Branch("pz_prt", &m_pz_prt);
160  m_outputTree->Branch("eta_prt", &m_eta_prt);
161  m_outputTree->Branch("pT_prt", &m_pT_prt);
162 
163  m_outputTree->Branch("nFiltered", &m_nFiltered);
164  m_outputTree->Branch("filtered", &m_flt);
165  m_outputTree->Branch("eLOC0_flt", &m_eLOC0_flt);
166  m_outputTree->Branch("eLOC1_flt", &m_eLOC1_flt);
167  m_outputTree->Branch("ePHI_flt", &m_ePHI_flt);
168  m_outputTree->Branch("eTHETA_flt", &m_eTHETA_flt);
169  m_outputTree->Branch("eQOP_flt", &m_eQOP_flt);
170  m_outputTree->Branch("eT_flt", &m_eT_flt);
171  m_outputTree->Branch("res_eLOC0_flt", &m_res_eLOC0_flt);
172  m_outputTree->Branch("res_eLOC1_flt", &m_res_eLOC1_flt);
173  m_outputTree->Branch("res_ePHI_flt", &m_res_ePHI_flt);
174  m_outputTree->Branch("res_eTHETA_flt", &m_res_eTHETA_flt);
175  m_outputTree->Branch("res_eQOP_flt", &m_res_eQOP_flt);
176  m_outputTree->Branch("res_eT_flt", &m_res_eT_flt);
177  m_outputTree->Branch("err_eLOC0_flt", &m_err_eLOC0_flt);
178  m_outputTree->Branch("err_eLOC1_flt", &m_err_eLOC1_flt);
179  m_outputTree->Branch("err_ePHI_flt", &m_err_ePHI_flt);
180  m_outputTree->Branch("err_eTHETA_flt", &m_err_eTHETA_flt);
181  m_outputTree->Branch("err_eQOP_flt", &m_err_eQOP_flt);
182  m_outputTree->Branch("err_eT_flt", &m_err_eT_flt);
183  m_outputTree->Branch("pull_eLOC0_flt", &m_pull_eLOC0_flt);
184  m_outputTree->Branch("pull_eLOC1_flt", &m_pull_eLOC1_flt);
185  m_outputTree->Branch("pull_ePHI_flt", &m_pull_ePHI_flt);
186  m_outputTree->Branch("pull_eTHETA_flt", &m_pull_eTHETA_flt);
187  m_outputTree->Branch("pull_eQOP_flt", &m_pull_eQOP_flt);
188  m_outputTree->Branch("pull_eT_flt", &m_pull_eT_flt);
189  m_outputTree->Branch("g_x_flt", &m_x_flt);
190  m_outputTree->Branch("g_y_flt", &m_y_flt);
191  m_outputTree->Branch("g_z_flt", &m_z_flt);
192  m_outputTree->Branch("px_flt", &m_px_flt);
193  m_outputTree->Branch("py_flt", &m_py_flt);
194  m_outputTree->Branch("pz_flt", &m_pz_flt);
195  m_outputTree->Branch("eta_flt", &m_eta_flt);
196  m_outputTree->Branch("pT_flt", &m_pT_flt);
197  m_outputTree->Branch("chi2", &m_chi2);
198 
199  m_outputTree->Branch("nSmoothed", &m_nSmoothed);
200  m_outputTree->Branch("smoothed", &m_smt);
201  m_outputTree->Branch("eLOC0_smt", &m_eLOC0_smt);
202  m_outputTree->Branch("eLOC1_smt", &m_eLOC1_smt);
203  m_outputTree->Branch("ePHI_smt", &m_ePHI_smt);
204  m_outputTree->Branch("eTHETA_smt", &m_eTHETA_smt);
205  m_outputTree->Branch("eQOP_smt", &m_eQOP_smt);
206  m_outputTree->Branch("eT_smt", &m_eT_smt);
207  m_outputTree->Branch("res_eLOC0_smt", &m_res_eLOC0_smt);
208  m_outputTree->Branch("res_eLOC1_smt", &m_res_eLOC1_smt);
209  m_outputTree->Branch("res_ePHI_smt", &m_res_ePHI_smt);
210  m_outputTree->Branch("res_eTHETA_smt", &m_res_eTHETA_smt);
211  m_outputTree->Branch("res_eQOP_smt", &m_res_eQOP_smt);
212  m_outputTree->Branch("res_eT_smt", &m_res_eT_smt);
213  m_outputTree->Branch("err_eLOC0_smt", &m_err_eLOC0_smt);
214  m_outputTree->Branch("err_eLOC1_smt", &m_err_eLOC1_smt);
215  m_outputTree->Branch("err_ePHI_smt", &m_err_ePHI_smt);
216  m_outputTree->Branch("err_eTHETA_smt", &m_err_eTHETA_smt);
217  m_outputTree->Branch("err_eQOP_smt", &m_err_eQOP_smt);
218  m_outputTree->Branch("err_eT_smt", &m_err_eT_smt);
219  m_outputTree->Branch("pull_eLOC0_smt", &m_pull_eLOC0_smt);
220  m_outputTree->Branch("pull_eLOC1_smt", &m_pull_eLOC1_smt);
221  m_outputTree->Branch("pull_ePHI_smt", &m_pull_ePHI_smt);
222  m_outputTree->Branch("pull_eTHETA_smt", &m_pull_eTHETA_smt);
223  m_outputTree->Branch("pull_eQOP_smt", &m_pull_eQOP_smt);
224  m_outputTree->Branch("pull_eT_smt", &m_pull_eT_smt);
225  m_outputTree->Branch("g_x_smt", &m_x_smt);
226  m_outputTree->Branch("g_y_smt", &m_y_smt);
227  m_outputTree->Branch("g_z_smt", &m_z_smt);
228  m_outputTree->Branch("px_smt", &m_px_smt);
229  m_outputTree->Branch("py_smt", &m_py_smt);
230  m_outputTree->Branch("pz_smt", &m_pz_smt);
231  m_outputTree->Branch("eta_smt", &m_eta_smt);
232  m_outputTree->Branch("pT_smt", &m_pT_smt);
233  }
234 }
235 
237  if (m_outputFile) {
238  m_outputFile->Close();
239  }
240 }
241 
243  if (m_outputFile) {
244  m_outputFile->cd();
245  m_outputTree->Write();
246  ACTS_INFO("Write trajectories to tree '"
247  << m_cfg.outputTreename << "' in '"
248  << joinPaths(m_cfg.outputDir, m_cfg.outputFilename) << "'");
249  }
250  return ProcessCode::SUCCESS;
251 }
252 
254  const AlgorithmContext& ctx, const TrajectoryContainer& trajectories) {
255  if (m_outputFile == nullptr)
256  return ProcessCode::SUCCESS;
257 
258  auto& gctx = ctx.geoContext;
259 
260  // read truth particles from input collection
261  const auto& particles =
262  ctx.eventStore.get<SimParticleContainer>(m_cfg.inputParticles);
263 
264  // Exclusive access to the tree while writing
265  std::lock_guard<std::mutex> lock(m_writeMutex);
266 
267  // Get the event number
268  m_eventNr = ctx.eventNumber;
269 
270  // Loop over the trajectories
271  int iTraj = 0;
272  for (const auto& traj : trajectories) {
273  m_trajNr = iTraj;
274 
275  // The trajectory entry indices and the multiTrajectory
276  const auto& [trackTips, mj] = traj.trajectory();
277  if (trackTips.empty()) {
278  ACTS_WARNING("Empty multiTrajectory.");
279  continue;
280  }
281  // Check the size of the trajectory entry indices. For track fitting, there
282  // should be at most one trajectory
283  if (trackTips.size() > 1) {
284  ACTS_ERROR("Track fitting should not result in multiple trajectories.");
285  return ProcessCode::ABORT;
286  }
287  // Get the entry index for the single trajectory
288  auto& trackTip = trackTips.front();
289 
290  // Collect the trajectory summary info
291  auto trajState =
293  m_nMeasurements = trajState.nMeasurements;
294  m_nStates = trajState.nStates;
295 
296  // Get the majority truth particle to this track
297  const auto particleHitCount = traj.identifyMajorityParticle(trackTip);
298  if (not particleHitCount.empty()) {
299  // Get the barcode of the majority truth particle
300  m_t_barcode = particleHitCount.front().particleId.value();
301  // Find the truth particle via the barcode
302  auto ip = particles.find(m_t_barcode);
303  if (ip != particles.end()) {
304  const auto& particle = *ip;
305  ACTS_DEBUG("Find the truth particle with barcode = " << m_t_barcode);
306  // Get the truth particle info at vertex
307  const auto p = particle.absMomentum();
308  m_t_charge = particle.charge();
309  m_t_time = particle.time();
310  m_t_vx = particle.position().x();
311  m_t_vy = particle.position().y();
312  m_t_vz = particle.position().z();
313  m_t_px = p * particle.unitDirection().x();
314  m_t_py = p * particle.unitDirection().y();
315  m_t_pz = p * particle.unitDirection().z();
316  m_t_theta = theta(particle.unitDirection());
317  m_t_phi = phi(particle.unitDirection());
318  m_t_eta = eta(particle.unitDirection());
319  m_t_pT = p * perp(particle.unitDirection());
320  } else {
321  ACTS_WARNING("Truth particle with barcode = " << m_t_barcode
322  << " not found!");
323  }
324  }
325 
326  // Get the fitted track parameter
327  m_hasFittedParams = false;
328  if (traj.hasTrackParameters(trackTip)) {
329  m_hasFittedParams = true;
330  const auto& boundParam = traj.trackParameters(trackTip);
331  const auto& parameter = boundParam.parameters();
332  const auto& covariance = *boundParam.covariance();
333  m_eLOC0_fit = parameter[Acts::eBoundLoc0];
334  m_eLOC1_fit = parameter[Acts::eBoundLoc1];
335  m_ePHI_fit = parameter[Acts::eBoundPhi];
336  m_eTHETA_fit = parameter[Acts::eBoundTheta];
337  m_eQOP_fit = parameter[Acts::eBoundQOverP];
338  m_eT_fit = parameter[Acts::eBoundTime];
339  m_err_eLOC0_fit = sqrt(covariance(Acts::eBoundLoc0, Acts::eBoundLoc0));
340  m_err_eLOC1_fit = sqrt(covariance(Acts::eBoundLoc1, Acts::eBoundLoc1));
341  m_err_ePHI_fit = sqrt(covariance(Acts::eBoundPhi, Acts::eBoundPhi));
342  m_err_eTHETA_fit = sqrt(covariance(Acts::eBoundTheta, Acts::eBoundTheta));
343  m_err_eQOP_fit = sqrt(covariance(Acts::eBoundQOverP, Acts::eBoundQOverP));
344  m_err_eT_fit = sqrt(covariance(Acts::eBoundTime, Acts::eBoundTime));
345  }
346 
347  // Get the trackStates on the trajectory
348  m_nPredicted = 0;
349  m_nFiltered = 0;
350  m_nSmoothed = 0;
351  mj.visitBackwards(trackTip, [&](const auto& state) {
352  // we only fill the track states with non-outlier measurement
353  auto typeFlags = state.typeFlags();
354  if (not typeFlags.test(Acts::TrackStateFlag::MeasurementFlag)) {
355  return true;
356  }
357 
358  auto meas = std::get<Measurement>(*state.uncalibrated());
359  auto& surface = meas.referenceObject();
360 
361  // get the geometry ID
362  auto geoID = surface.geometryId();
363  m_volumeID.push_back(geoID.volume());
364  m_layerID.push_back(geoID.layer());
365  m_moduleID.push_back(geoID.sensitive());
366 
367  // get local position
368  Acts::Vector2D local(meas.parameters()[Acts::eBoundLoc0],
369  meas.parameters()[Acts::eBoundLoc1]);
370  // get global position
371  Acts::Vector3D mom(1, 1, 1);
372  Acts::Vector3D global = surface.localToGlobal(ctx.geoContext, local, mom);
373 
374  // get measurement covariance
375  auto cov = meas.covariance();
376  // float resX = sqrt(cov(Acts::eBoundLoc0, Acts::eBoundLoc0));
377  // float resY = sqrt(cov(Acts::eBoundLoc1, Acts::eBoundLoc1));
378 
379  // push the measurement info
380  m_lx_hit.push_back(local.x());
381  m_ly_hit.push_back(local.y());
382  m_x_hit.push_back(global.x());
383  m_y_hit.push_back(global.y());
384  m_z_hit.push_back(global.z());
385 
386  // get the truth hit corresponding to this trackState
387  const auto& truthHit = state.uncalibrated().truthHit();
388  // get local truth position
389  Acts::Vector2D truthlocal{0., 0.};
390  auto lpResult = surface.globalToLocal(gctx, truthHit.position(),
391  truthHit.unitDirection());
392  if (not lpResult.ok()) {
393  ACTS_FATAL("Global to local transformation did not succeed.");
394  }
395  truthlocal = lpResult.value();
396 
397  // push the truth hit info
398  m_t_x.push_back(truthHit.position().x());
399  m_t_y.push_back(truthHit.position().y());
400  m_t_z.push_back(truthHit.position().z());
401  m_t_r.push_back(perp(truthHit.position()));
402  m_t_dx.push_back(truthHit.unitDirection().x());
403  m_t_dy.push_back(truthHit.unitDirection().y());
404  m_t_dz.push_back(truthHit.unitDirection().z());
405 
406  // get the truth track parameter at this track State
407  float truthLOC0 = 0, truthLOC1 = 0, truthPHI = 0, truthTHETA = 0,
408  truthQOP = 0, truthTIME = 0;
409  truthLOC0 = truthlocal.x();
410  truthLOC1 = truthlocal.y();
411  truthPHI = phi(truthHit.unitDirection());
412  truthTHETA = theta(truthHit.unitDirection());
413  truthQOP =
414  m_t_charge / truthHit.momentum4Before().template head<3>().norm();
415  truthTIME = truthHit.time();
416 
417  // push the truth track parameter at this track State
418  m_t_eLOC0.push_back(truthLOC0);
419  m_t_eLOC1.push_back(truthLOC1);
420  m_t_ePHI.push_back(truthPHI);
421  m_t_eTHETA.push_back(truthTHETA);
422  m_t_eQOP.push_back(truthQOP);
423  m_t_eT.push_back(truthTIME);
424 
425  // get the predicted parameter
426  bool predicted = false;
427  if (state.hasPredicted()) {
428  predicted = true;
429  m_nPredicted++;
430  auto parameters = state.predicted();
431  auto covariance = state.predictedCovariance();
432  // local hit residual info
433  auto H = meas.projector();
434  auto resCov = cov + H * covariance * H.transpose();
435  auto residual = meas.residual(parameters);
436  m_res_x_hit.push_back(residual(Acts::eBoundLoc0));
437  m_res_y_hit.push_back(residual(Acts::eBoundLoc1));
438  m_err_x_hit.push_back(sqrt(resCov(Acts::eBoundLoc0, Acts::eBoundLoc0)));
439  m_err_y_hit.push_back(sqrt(resCov(Acts::eBoundLoc1, Acts::eBoundLoc1)));
440  m_pull_x_hit.push_back(
441  residual(Acts::eBoundLoc0) /
442  sqrt(resCov(Acts::eBoundLoc0, Acts::eBoundLoc0)));
443  m_pull_y_hit.push_back(
444  residual(Acts::eBoundLoc1) /
445  sqrt(resCov(Acts::eBoundLoc1, Acts::eBoundLoc1)));
446  m_dim_hit.push_back(state.calibratedSize());
447 
448  // predicted parameter
449  m_eLOC0_prt.push_back(parameters[Acts::eBoundLoc0]);
450  m_eLOC1_prt.push_back(parameters[Acts::eBoundLoc1]);
451  m_ePHI_prt.push_back(parameters[Acts::eBoundPhi]);
452  m_eTHETA_prt.push_back(parameters[Acts::eBoundTheta]);
453  m_eQOP_prt.push_back(parameters[Acts::eBoundQOverP]);
454  m_eT_prt.push_back(parameters[Acts::eBoundTime]);
455 
456  // predicted residual
457  m_res_eLOC0_prt.push_back(parameters[Acts::eBoundLoc0] - truthLOC0);
458  m_res_eLOC1_prt.push_back(parameters[Acts::eBoundLoc1] - truthLOC1);
459  m_res_ePHI_prt.push_back(parameters[Acts::eBoundPhi] - truthPHI);
460  m_res_eTHETA_prt.push_back(parameters[Acts::eBoundTheta] - truthTHETA);
461  m_res_eQOP_prt.push_back(parameters[Acts::eBoundQOverP] - truthQOP);
462  m_res_eT_prt.push_back(parameters[Acts::eBoundTime] - truthTIME);
463 
464  // predicted parameter error
465  m_err_eLOC0_prt.push_back(
466  sqrt(covariance(Acts::eBoundLoc0, Acts::eBoundLoc0)));
467  m_err_eLOC1_prt.push_back(
468  sqrt(covariance(Acts::eBoundLoc1, Acts::eBoundLoc1)));
469  m_err_ePHI_prt.push_back(
470  sqrt(covariance(Acts::eBoundPhi, Acts::eBoundPhi)));
471  m_err_eTHETA_prt.push_back(
472  sqrt(covariance(Acts::eBoundTheta, Acts::eBoundTheta)));
473  m_err_eQOP_prt.push_back(
474  sqrt(covariance(Acts::eBoundQOverP, Acts::eBoundQOverP)));
475  m_err_eT_prt.push_back(
476  sqrt(covariance(Acts::eBoundTime, Acts::eBoundTime)));
477 
478  // predicted parameter pull
479  m_pull_eLOC0_prt.push_back(
480  (parameters[Acts::eBoundLoc0] - truthLOC0) /
481  sqrt(covariance(Acts::eBoundLoc0, Acts::eBoundLoc0)));
482  m_pull_eLOC1_prt.push_back(
483  (parameters[Acts::eBoundLoc1] - truthLOC1) /
484  sqrt(covariance(Acts::eBoundLoc1, Acts::eBoundLoc1)));
485  m_pull_ePHI_prt.push_back(
486  (parameters[Acts::eBoundPhi] - truthPHI) /
487  sqrt(covariance(Acts::eBoundPhi, Acts::eBoundPhi)));
488  m_pull_eTHETA_prt.push_back(
489  (parameters[Acts::eBoundTheta] - truthTHETA) /
490  sqrt(covariance(Acts::eBoundTheta, Acts::eBoundTheta)));
491  m_pull_eQOP_prt.push_back(
492  (parameters[Acts::eBoundQOverP] - truthQOP) /
493  sqrt(covariance(Acts::eBoundQOverP, Acts::eBoundQOverP)));
494  m_pull_eT_prt.push_back(
495  (parameters[Acts::eBoundTime] - truthTIME) /
496  sqrt(covariance(Acts::eBoundTime, Acts::eBoundTime)));
497 
498  // further predicted parameter info
499  Acts::FreeVector freeParams =
501  parameters);
502  m_x_prt.push_back(freeParams[Acts::eFreePos0]);
503  m_y_prt.push_back(freeParams[Acts::eFreePos1]);
504  m_z_prt.push_back(freeParams[Acts::eFreePos2]);
505  auto p = std::abs(1 / freeParams[Acts::eFreeQOverP]);
506  m_px_prt.push_back(p * freeParams[Acts::eFreeDir0]);
507  m_py_prt.push_back(p * freeParams[Acts::eFreeDir1]);
508  m_pz_prt.push_back(p * freeParams[Acts::eFreeDir2]);
509  m_pT_prt.push_back(p * std::hypot(freeParams[Acts::eFreeDir0],
510  freeParams[Acts::eFreeDir1]));
511  m_eta_prt.push_back(
512  Acts::VectorHelpers::eta(freeParams.segment<3>(Acts::eFreeDir0)));
513  } else {
514  // push default values if no predicted parameter
515  m_res_x_hit.push_back(-99.);
516  m_res_y_hit.push_back(-99.);
517  m_err_x_hit.push_back(-99.);
518  m_err_y_hit.push_back(-99.);
519  m_pull_x_hit.push_back(-99.);
520  m_pull_y_hit.push_back(-99.);
521  m_dim_hit.push_back(-99.);
522  m_eLOC0_prt.push_back(-99.);
523  m_eLOC1_prt.push_back(-99.);
524  m_ePHI_prt.push_back(-99.);
525  m_eTHETA_prt.push_back(-99.);
526  m_eQOP_prt.push_back(-99.);
527  m_eT_prt.push_back(-99.);
528  m_res_eLOC0_prt.push_back(-99.);
529  m_res_eLOC1_prt.push_back(-99.);
530  m_res_ePHI_prt.push_back(-99.);
531  m_res_eTHETA_prt.push_back(-99.);
532  m_res_eQOP_prt.push_back(-99.);
533  m_res_eT_prt.push_back(-99.);
534  m_err_eLOC0_prt.push_back(-99);
535  m_err_eLOC1_prt.push_back(-99);
536  m_err_ePHI_prt.push_back(-99);
537  m_err_eTHETA_prt.push_back(-99);
538  m_err_eQOP_prt.push_back(-99);
539  m_err_eT_prt.push_back(-99);
540  m_pull_eLOC0_prt.push_back(-99.);
541  m_pull_eLOC1_prt.push_back(-99.);
542  m_pull_ePHI_prt.push_back(-99.);
543  m_pull_eTHETA_prt.push_back(-99.);
544  m_pull_eQOP_prt.push_back(-99.);
545  m_pull_eT_prt.push_back(-99.);
546  m_x_prt.push_back(-99.);
547  m_y_prt.push_back(-99.);
548  m_z_prt.push_back(-99.);
549  m_px_prt.push_back(-99.);
550  m_py_prt.push_back(-99.);
551  m_pz_prt.push_back(-99.);
552  m_pT_prt.push_back(-99.);
553  m_eta_prt.push_back(-99.);
554  }
555 
556  // get the filtered parameter
557  bool filtered = false;
558  if (state.hasFiltered()) {
559  filtered = true;
560  m_nFiltered++;
561  auto parameters = state.filtered();
562  auto covariance = state.filteredCovariance();
563  // filtered parameter
564  m_eLOC0_flt.push_back(parameters[Acts::eBoundLoc0]);
565  m_eLOC1_flt.push_back(parameters[Acts::eBoundLoc1]);
566  m_ePHI_flt.push_back(parameters[Acts::eBoundPhi]);
567  m_eTHETA_flt.push_back(parameters[Acts::eBoundTheta]);
568  m_eQOP_flt.push_back(parameters[Acts::eBoundQOverP]);
569  m_eT_flt.push_back(parameters[Acts::eBoundTime]);
570 
571  // filtered residual
572  m_res_eLOC0_flt.push_back(parameters[Acts::eBoundLoc0] - truthLOC0);
573  m_res_eLOC1_flt.push_back(parameters[Acts::eBoundLoc1] - truthLOC1);
574  m_res_ePHI_flt.push_back(parameters[Acts::eBoundPhi] - truthPHI);
575  m_res_eTHETA_flt.push_back(parameters[Acts::eBoundTheta] - truthTHETA);
576  m_res_eQOP_flt.push_back(parameters[Acts::eBoundQOverP] - truthQOP);
577  m_res_eT_flt.push_back(parameters[Acts::eBoundTime] - truthTIME);
578 
579  // filtered parameter error
580  m_err_eLOC0_flt.push_back(
581  sqrt(covariance(Acts::eBoundLoc0, Acts::eBoundLoc0)));
582  m_err_eLOC1_flt.push_back(
583  sqrt(covariance(Acts::eBoundLoc1, Acts::eBoundLoc1)));
584  m_err_ePHI_flt.push_back(
585  sqrt(covariance(Acts::eBoundPhi, Acts::eBoundPhi)));
586  m_err_eTHETA_flt.push_back(
587  sqrt(covariance(Acts::eBoundTheta, Acts::eBoundTheta)));
588  m_err_eQOP_flt.push_back(
589  sqrt(covariance(Acts::eBoundQOverP, Acts::eBoundQOverP)));
590  m_err_eT_flt.push_back(
591  sqrt(covariance(Acts::eBoundTime, Acts::eBoundTime)));
592 
593  // filtered parameter pull
594  m_pull_eLOC0_flt.push_back(
595  (parameters[Acts::eBoundLoc0] - truthLOC0) /
596  sqrt(covariance(Acts::eBoundLoc0, Acts::eBoundLoc0)));
597  m_pull_eLOC1_flt.push_back(
598  (parameters[Acts::eBoundLoc1] - truthLOC1) /
599  sqrt(covariance(Acts::eBoundLoc1, Acts::eBoundLoc1)));
600  m_pull_ePHI_flt.push_back(
601  (parameters[Acts::eBoundPhi] - truthPHI) /
602  sqrt(covariance(Acts::eBoundPhi, Acts::eBoundPhi)));
603  m_pull_eTHETA_flt.push_back(
604  (parameters[Acts::eBoundTheta] - truthTHETA) /
605  sqrt(covariance(Acts::eBoundTheta, Acts::eBoundTheta)));
606  m_pull_eQOP_flt.push_back(
607  (parameters[Acts::eBoundQOverP] - truthQOP) /
608  sqrt(covariance(Acts::eBoundQOverP, Acts::eBoundQOverP)));
609  m_pull_eT_flt.push_back(
610  (parameters[Acts::eBoundTime] - truthTIME) /
611  sqrt(covariance(Acts::eBoundTime, Acts::eBoundTime)));
612 
613  // more filtered parameter info
614  const Acts::FreeVector freeParams =
616  parameters);
617  m_x_flt.push_back(freeParams[Acts::eFreePos0]);
618  m_y_flt.push_back(freeParams[Acts::eFreePos1]);
619  m_z_flt.push_back(freeParams[Acts::eFreePos2]);
620  const auto p = std::abs(1 / freeParams[Acts::eFreeQOverP]);
621  m_px_flt.push_back(p * freeParams[Acts::eFreeDir0]);
622  m_py_flt.push_back(p * freeParams[Acts::eFreeDir1]);
623  m_pz_flt.push_back(p * freeParams[Acts::eFreeDir2]);
624  m_pT_flt.push_back(p * std::hypot(freeParams[Acts::eFreeDir0],
625  freeParams[Acts::eFreeDir1]));
626  m_eta_flt.push_back(
627  Acts::VectorHelpers::eta(freeParams.segment<3>(Acts::eFreeDir0)));
628  m_chi2.push_back(state.chi2());
629  } else {
630  // push default values if no filtered parameter
631  m_eLOC0_flt.push_back(-99.);
632  m_eLOC1_flt.push_back(-99.);
633  m_ePHI_flt.push_back(-99.);
634  m_eTHETA_flt.push_back(-99.);
635  m_eQOP_flt.push_back(-99.);
636  m_eT_flt.push_back(-99.);
637  m_res_eLOC0_flt.push_back(-99.);
638  m_res_eLOC1_flt.push_back(-99.);
639  m_res_ePHI_flt.push_back(-99.);
640  m_res_eTHETA_flt.push_back(-99.);
641  m_res_eQOP_flt.push_back(-99.);
642  m_res_eT_flt.push_back(-99.);
643  m_err_eLOC0_flt.push_back(-99);
644  m_err_eLOC1_flt.push_back(-99);
645  m_err_ePHI_flt.push_back(-99);
646  m_err_eTHETA_flt.push_back(-99);
647  m_err_eQOP_flt.push_back(-99);
648  m_err_eT_flt.push_back(-99);
649  m_pull_eLOC0_flt.push_back(-99.);
650  m_pull_eLOC1_flt.push_back(-99.);
651  m_pull_ePHI_flt.push_back(-99.);
652  m_pull_eTHETA_flt.push_back(-99.);
653  m_pull_eQOP_flt.push_back(-99.);
654  m_pull_eT_flt.push_back(-99.);
655  m_x_flt.push_back(-99.);
656  m_y_flt.push_back(-99.);
657  m_z_flt.push_back(-99.);
658  m_py_flt.push_back(-99.);
659  m_pz_flt.push_back(-99.);
660  m_pT_flt.push_back(-99.);
661  m_eta_flt.push_back(-99.);
662  m_chi2.push_back(-99.0);
663  }
664 
665  // get the smoothed parameter
666  bool smoothed = false;
667  if (state.hasSmoothed()) {
668  smoothed = true;
669  m_nSmoothed++;
670  auto parameters = state.smoothed();
671  auto covariance = state.smoothedCovariance();
672 
673  // smoothed parameter
674  m_eLOC0_smt.push_back(parameters[Acts::eBoundLoc0]);
675  m_eLOC1_smt.push_back(parameters[Acts::eBoundLoc1]);
676  m_ePHI_smt.push_back(parameters[Acts::eBoundPhi]);
677  m_eTHETA_smt.push_back(parameters[Acts::eBoundTheta]);
678  m_eQOP_smt.push_back(parameters[Acts::eBoundQOverP]);
679  m_eT_smt.push_back(parameters[Acts::eBoundTime]);
680 
681  // smoothed residual
682  m_res_eLOC0_smt.push_back(parameters[Acts::eBoundLoc0] - truthLOC0);
683  m_res_eLOC1_smt.push_back(parameters[Acts::eBoundLoc1] - truthLOC1);
684  m_res_ePHI_smt.push_back(parameters[Acts::eBoundPhi] - truthPHI);
685  m_res_eTHETA_smt.push_back(parameters[Acts::eBoundTheta] - truthTHETA);
686  m_res_eQOP_smt.push_back(parameters[Acts::eBoundQOverP] - truthQOP);
687  m_res_eT_smt.push_back(parameters[Acts::eBoundTime] - truthTIME);
688 
689  // smoothed parameter error
690  m_err_eLOC0_smt.push_back(
691  sqrt(covariance(Acts::eBoundLoc0, Acts::eBoundLoc0)));
692  m_err_eLOC1_smt.push_back(
693  sqrt(covariance(Acts::eBoundLoc1, Acts::eBoundLoc1)));
694  m_err_ePHI_smt.push_back(
695  sqrt(covariance(Acts::eBoundPhi, Acts::eBoundPhi)));
696  m_err_eTHETA_smt.push_back(
697  sqrt(covariance(Acts::eBoundTheta, Acts::eBoundTheta)));
698  m_err_eQOP_smt.push_back(
699  sqrt(covariance(Acts::eBoundQOverP, Acts::eBoundQOverP)));
700  m_err_eT_smt.push_back(
701  sqrt(covariance(Acts::eBoundTime, Acts::eBoundTime)));
702 
703  // smoothed parameter pull
704  m_pull_eLOC0_smt.push_back(
705  (parameters[Acts::eBoundLoc0] - truthLOC0) /
706  sqrt(covariance(Acts::eBoundLoc0, Acts::eBoundLoc0)));
707  m_pull_eLOC1_smt.push_back(
708  (parameters[Acts::eBoundLoc1] - truthLOC1) /
709  sqrt(covariance(Acts::eBoundLoc1, Acts::eBoundLoc1)));
710  m_pull_ePHI_smt.push_back(
711  (parameters[Acts::eBoundPhi] - truthPHI) /
712  sqrt(covariance(Acts::eBoundPhi, Acts::eBoundPhi)));
713  m_pull_eTHETA_smt.push_back(
714  (parameters[Acts::eBoundTheta] - truthTHETA) /
715  sqrt(covariance(Acts::eBoundTheta, Acts::eBoundTheta)));
716  m_pull_eQOP_smt.push_back(
717  (parameters[Acts::eBoundQOverP] - truthQOP) /
718  sqrt(covariance(Acts::eBoundQOverP, Acts::eBoundQOverP)));
719  m_pull_eT_smt.push_back(
720  (parameters[Acts::eBoundTime] - truthTIME) /
721  sqrt(covariance(Acts::eBoundTime, Acts::eBoundTime)));
722 
723  // further smoothed parameter info
724  const Acts::FreeVector freeParams =
726  parameters);
727  m_x_smt.push_back(freeParams[Acts::eFreePos0]);
728  m_y_smt.push_back(freeParams[Acts::eFreePos1]);
729  m_z_smt.push_back(freeParams[Acts::eFreePos2]);
730  const auto p = std::abs(1 / freeParams[Acts::eFreeQOverP]);
731  m_px_smt.push_back(p * freeParams[Acts::eFreeDir0]);
732  m_py_smt.push_back(p * freeParams[Acts::eFreeDir1]);
733  m_pz_smt.push_back(p * freeParams[Acts::eFreeDir2]);
734  m_pT_smt.push_back(p * std::hypot(freeParams[Acts::eFreeDir0],
735  freeParams[Acts::eFreeDir1]));
736  m_eta_smt.push_back(
737  Acts::VectorHelpers::eta(freeParams.segment<3>(Acts::eFreeDir0)));
738  } else {
739  // push default values if no smoothed parameter
740  m_eLOC0_smt.push_back(-99.);
741  m_eLOC1_smt.push_back(-99.);
742  m_ePHI_smt.push_back(-99.);
743  m_eTHETA_smt.push_back(-99.);
744  m_eQOP_smt.push_back(-99.);
745  m_eT_smt.push_back(-99.);
746  m_res_eLOC0_smt.push_back(-99.);
747  m_res_eLOC1_smt.push_back(-99.);
748  m_res_ePHI_smt.push_back(-99.);
749  m_res_eTHETA_smt.push_back(-99.);
750  m_res_eQOP_smt.push_back(-99.);
751  m_res_eT_smt.push_back(-99.);
752  m_err_eLOC0_smt.push_back(-99);
753  m_err_eLOC1_smt.push_back(-99);
754  m_err_ePHI_smt.push_back(-99);
755  m_err_eTHETA_smt.push_back(-99);
756  m_err_eQOP_smt.push_back(-99);
757  m_err_eT_smt.push_back(-99);
758  m_pull_eLOC0_smt.push_back(-99.);
759  m_pull_eLOC1_smt.push_back(-99.);
760  m_pull_ePHI_smt.push_back(-99.);
761  m_pull_eTHETA_smt.push_back(-99.);
762  m_pull_eQOP_smt.push_back(-99.);
763  m_pull_eT_smt.push_back(-99.);
764  m_x_smt.push_back(-99.);
765  m_y_smt.push_back(-99.);
766  m_z_smt.push_back(-99.);
767  m_px_smt.push_back(-99.);
768  m_py_smt.push_back(-99.);
769  m_pz_smt.push_back(-99.);
770  m_pT_smt.push_back(-99.);
771  m_eta_smt.push_back(-99.);
772  }
773 
774  m_prt.push_back(predicted);
775  m_flt.push_back(filtered);
776  m_smt.push_back(smoothed);
777  return true;
778  }); // all states
779 
780  // fill the variables for one track to tree
781  m_outputTree->Fill();
782 
783  // now reset
784  m_t_x.clear();
785  m_t_y.clear();
786  m_t_z.clear();
787  m_t_r.clear();
788  m_t_dx.clear();
789  m_t_dy.clear();
790  m_t_dz.clear();
791  m_t_eLOC0.clear();
792  m_t_eLOC1.clear();
793  m_t_ePHI.clear();
794  m_t_eTHETA.clear();
795  m_t_eQOP.clear();
796  m_t_eT.clear();
797 
798  m_volumeID.clear();
799  m_layerID.clear();
800  m_moduleID.clear();
801  m_lx_hit.clear();
802  m_ly_hit.clear();
803  m_x_hit.clear();
804  m_y_hit.clear();
805  m_z_hit.clear();
806  m_res_x_hit.clear();
807  m_res_y_hit.clear();
808  m_err_x_hit.clear();
809  m_err_y_hit.clear();
810  m_pull_x_hit.clear();
811  m_pull_y_hit.clear();
812  m_dim_hit.clear();
813 
814  m_prt.clear();
815  m_eLOC0_prt.clear();
816  m_eLOC1_prt.clear();
817  m_ePHI_prt.clear();
818  m_eTHETA_prt.clear();
819  m_eQOP_prt.clear();
820  m_eT_prt.clear();
821  m_res_eLOC0_prt.clear();
822  m_res_eLOC1_prt.clear();
823  m_res_ePHI_prt.clear();
824  m_res_eTHETA_prt.clear();
825  m_res_eQOP_prt.clear();
826  m_res_eT_prt.clear();
827  m_err_eLOC0_prt.clear();
828  m_err_eLOC1_prt.clear();
829  m_err_ePHI_prt.clear();
830  m_err_eTHETA_prt.clear();
831  m_err_eQOP_prt.clear();
832  m_err_eT_prt.clear();
833  m_pull_eLOC0_prt.clear();
834  m_pull_eLOC1_prt.clear();
835  m_pull_ePHI_prt.clear();
836  m_pull_eTHETA_prt.clear();
837  m_pull_eQOP_prt.clear();
838  m_pull_eT_prt.clear();
839  m_x_prt.clear();
840  m_y_prt.clear();
841  m_z_prt.clear();
842  m_px_prt.clear();
843  m_py_prt.clear();
844  m_pz_prt.clear();
845  m_eta_prt.clear();
846  m_pT_prt.clear();
847 
848  m_flt.clear();
849  m_eLOC0_flt.clear();
850  m_eLOC1_flt.clear();
851  m_ePHI_flt.clear();
852  m_eTHETA_flt.clear();
853  m_eQOP_flt.clear();
854  m_eT_flt.clear();
855  m_res_eLOC0_flt.clear();
856  m_res_eLOC1_flt.clear();
857  m_res_ePHI_flt.clear();
858  m_res_eTHETA_flt.clear();
859  m_res_eQOP_flt.clear();
860  m_res_eT_flt.clear();
861  m_err_eLOC0_flt.clear();
862  m_err_eLOC1_flt.clear();
863  m_err_ePHI_flt.clear();
864  m_err_eTHETA_flt.clear();
865  m_err_eQOP_flt.clear();
866  m_err_eT_flt.clear();
867  m_pull_eLOC0_flt.clear();
868  m_pull_eLOC1_flt.clear();
869  m_pull_ePHI_flt.clear();
870  m_pull_eTHETA_flt.clear();
871  m_pull_eQOP_flt.clear();
872  m_pull_eT_flt.clear();
873  m_x_flt.clear();
874  m_y_flt.clear();
875  m_z_flt.clear();
876  m_px_flt.clear();
877  m_py_flt.clear();
878  m_pz_flt.clear();
879  m_eta_flt.clear();
880  m_pT_flt.clear();
881  m_chi2.clear();
882 
883  m_smt.clear();
884  m_eLOC0_smt.clear();
885  m_eLOC1_smt.clear();
886  m_ePHI_smt.clear();
887  m_eTHETA_smt.clear();
888  m_eQOP_smt.clear();
889  m_eT_smt.clear();
890  m_res_eLOC0_smt.clear();
891  m_res_eLOC1_smt.clear();
892  m_res_ePHI_smt.clear();
893  m_res_eTHETA_smt.clear();
894  m_res_eQOP_smt.clear();
895  m_res_eT_smt.clear();
896  m_err_eLOC0_smt.clear();
897  m_err_eLOC1_smt.clear();
898  m_err_ePHI_smt.clear();
899  m_err_eTHETA_smt.clear();
900  m_err_eQOP_smt.clear();
901  m_err_eT_smt.clear();
902  m_pull_eLOC0_smt.clear();
903  m_pull_eLOC1_smt.clear();
904  m_pull_ePHI_smt.clear();
905  m_pull_eTHETA_smt.clear();
906  m_pull_eQOP_smt.clear();
907  m_pull_eT_smt.clear();
908  m_x_smt.clear();
909  m_y_smt.clear();
910  m_z_smt.clear();
911  m_px_smt.clear();
912  m_py_smt.clear();
913  m_pz_smt.clear();
914  m_eta_smt.clear();
915  m_pT_smt.clear();
916 
917  iTraj++;
918  } // all trajectories
919 
920  return ProcessCode::SUCCESS;
921 }