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