Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
RootTrajectoryStatesWriter.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file RootTrajectoryStatesWriter.cpp
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2019-2023 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 
9 #include "ActsExamples/Io/Root/RootTrajectoryStatesWriter.hpp"
10 
11 #include "Acts/Definitions/Algebra.hpp"
12 #include "Acts/Definitions/Common.hpp"
13 #include "Acts/Definitions/TrackParametrization.hpp"
14 #include "Acts/EventData/MultiTrajectory.hpp"
15 #include "Acts/EventData/MultiTrajectoryHelpers.hpp"
16 #include "Acts/EventData/VectorMultiTrajectory.hpp"
17 #include "Acts/EventData/detail/TransformationBoundToFree.hpp"
18 #include "Acts/Geometry/GeometryContext.hpp"
19 #include "Acts/Utilities/MultiIndex.hpp"
20 #include "Acts/Utilities/detail/periodic.hpp"
21 #include "ActsExamples/EventData/AverageSimHits.hpp"
22 #include "ActsExamples/Framework/AlgorithmContext.hpp"
23 #include "ActsExamples/Utilities/Range.hpp"
24 #include "ActsExamples/Validation/TrackClassification.hpp"
25 #include "ActsFatras/EventData/Barcode.hpp"
26 #include "ActsFatras/EventData/Particle.hpp"
27 
28 #include <cmath>
29 #include <cstddef>
30 #include <ios>
31 #include <limits>
32 #include <memory>
33 #include <optional>
34 #include <ostream>
35 #include <stdexcept>
36 #include <utility>
37 
38 #include <TFile.h>
39 #include <TTree.h>
40 
41 namespace ActsExamples {
42 class IndexSourceLink;
43 } // namespace ActsExamples
44 
45 using Acts::VectorHelpers::eta;
46 using Acts::VectorHelpers::perp;
47 using Acts::VectorHelpers::phi;
48 using Acts::VectorHelpers::theta;
49 
50 ActsExamples::RootTrajectoryStatesWriter::RootTrajectoryStatesWriter(
51  const ActsExamples::RootTrajectoryStatesWriter::Config& config,
52  Acts::Logging::Level level)
53  : WriterT(config.inputTrajectories, "RootTrajectoryStatesWriter", level),
54  m_cfg(config) {
55  // trajectories collection name is already checked by base ctor
56  if (m_cfg.inputParticles.empty()) {
57  throw std::invalid_argument("Missing particles input collection");
58  }
59  if (m_cfg.inputSimHits.empty()) {
60  throw std::invalid_argument("Missing simulated hits input collection");
61  }
62  if (m_cfg.inputMeasurementParticlesMap.empty()) {
63  throw std::invalid_argument("Missing hit-particles map input collection");
64  }
65  if (m_cfg.inputMeasurementSimHitsMap.empty()) {
66  throw std::invalid_argument(
67  "Missing hit-simulated-hits map input collection");
68  }
69  if (m_cfg.filePath.empty()) {
70  throw std::invalid_argument("Missing output filename");
71  }
72  if (m_cfg.treeName.empty()) {
73  throw std::invalid_argument("Missing tree name");
74  }
75 
76  m_inputParticles.initialize(m_cfg.inputParticles);
77  m_inputSimHits.initialize(m_cfg.inputSimHits);
78  m_inputMeasurementParticlesMap.initialize(m_cfg.inputMeasurementParticlesMap);
79  m_inputMeasurementSimHitsMap.initialize(m_cfg.inputMeasurementSimHitsMap);
80 
81  // Setup ROOT I/O
82  auto path = m_cfg.filePath;
83  m_outputFile = TFile::Open(path.c_str(), m_cfg.fileMode.c_str());
84  if (m_outputFile == nullptr) {
85  throw std::ios_base::failure("Could not open '" + path);
86  }
87  m_outputFile->cd();
88  m_outputTree = new TTree(m_cfg.treeName.c_str(), m_cfg.treeName.c_str());
89  if (m_outputTree == nullptr) {
90  throw std::bad_alloc();
91  } else {
92  // I/O parameters
93  m_outputTree->Branch("event_nr", &m_eventNr);
94  m_outputTree->Branch("multiTraj_nr", &m_multiTrajNr);
95  m_outputTree->Branch("subTraj_nr", &m_subTrajNr);
96 
97  m_outputTree->Branch("t_x", &m_t_x);
98  m_outputTree->Branch("t_y", &m_t_y);
99  m_outputTree->Branch("t_z", &m_t_z);
100  m_outputTree->Branch("t_r", &m_t_r);
101  m_outputTree->Branch("t_dx", &m_t_dx);
102  m_outputTree->Branch("t_dy", &m_t_dy);
103  m_outputTree->Branch("t_dz", &m_t_dz);
104  m_outputTree->Branch("t_eLOC0", &m_t_eLOC0);
105  m_outputTree->Branch("t_eLOC1", &m_t_eLOC1);
106  m_outputTree->Branch("t_ePHI", &m_t_ePHI);
107  m_outputTree->Branch("t_eTHETA", &m_t_eTHETA);
108  m_outputTree->Branch("t_eQOP", &m_t_eQOP);
109  m_outputTree->Branch("t_eT", &m_t_eT);
110 
111  m_outputTree->Branch("nStates", &m_nStates);
112  m_outputTree->Branch("nMeasurements", &m_nMeasurements);
113  m_outputTree->Branch("volume_id", &m_volumeID);
114  m_outputTree->Branch("layer_id", &m_layerID);
115  m_outputTree->Branch("module_id", &m_moduleID);
116  m_outputTree->Branch("pathLength", &m_pathLength);
117  m_outputTree->Branch("l_x_hit", &m_lx_hit);
118  m_outputTree->Branch("l_y_hit", &m_ly_hit);
119  m_outputTree->Branch("g_x_hit", &m_x_hit);
120  m_outputTree->Branch("g_y_hit", &m_y_hit);
121  m_outputTree->Branch("g_z_hit", &m_z_hit);
122  m_outputTree->Branch("res_x_hit", &m_res_x_hit);
123  m_outputTree->Branch("res_y_hit", &m_res_y_hit);
124  m_outputTree->Branch("err_x_hit", &m_err_x_hit);
125  m_outputTree->Branch("err_y_hit", &m_err_y_hit);
126  m_outputTree->Branch("pull_x_hit", &m_pull_x_hit);
127  m_outputTree->Branch("pull_y_hit", &m_pull_y_hit);
128  m_outputTree->Branch("dim_hit", &m_dim_hit);
129 
130  m_outputTree->Branch("nPredicted", &m_nParams[0]);
131  m_outputTree->Branch("predicted", &m_hasParams[0]);
132  m_outputTree->Branch("eLOC0_prt", &m_eLOC0[0]);
133  m_outputTree->Branch("eLOC1_prt", &m_eLOC1[0]);
134  m_outputTree->Branch("ePHI_prt", &m_ePHI[0]);
135  m_outputTree->Branch("eTHETA_prt", &m_eTHETA[0]);
136  m_outputTree->Branch("eQOP_prt", &m_eQOP[0]);
137  m_outputTree->Branch("eT_prt", &m_eT[0]);
138  m_outputTree->Branch("res_eLOC0_prt", &m_res_eLOC0[0]);
139  m_outputTree->Branch("res_eLOC1_prt", &m_res_eLOC1[0]);
140  m_outputTree->Branch("res_ePHI_prt", &m_res_ePHI[0]);
141  m_outputTree->Branch("res_eTHETA_prt", &m_res_eTHETA[0]);
142  m_outputTree->Branch("res_eQOP_prt", &m_res_eQOP[0]);
143  m_outputTree->Branch("res_eT_prt", &m_res_eT[0]);
144  m_outputTree->Branch("err_eLOC0_prt", &m_err_eLOC0[0]);
145  m_outputTree->Branch("err_eLOC1_prt", &m_err_eLOC1[0]);
146  m_outputTree->Branch("err_ePHI_prt", &m_err_ePHI[0]);
147  m_outputTree->Branch("err_eTHETA_prt", &m_err_eTHETA[0]);
148  m_outputTree->Branch("err_eQOP_prt", &m_err_eQOP[0]);
149  m_outputTree->Branch("err_eT_prt", &m_err_eT[0]);
150  m_outputTree->Branch("pull_eLOC0_prt", &m_pull_eLOC0[0]);
151  m_outputTree->Branch("pull_eLOC1_prt", &m_pull_eLOC1[0]);
152  m_outputTree->Branch("pull_ePHI_prt", &m_pull_ePHI[0]);
153  m_outputTree->Branch("pull_eTHETA_prt", &m_pull_eTHETA[0]);
154  m_outputTree->Branch("pull_eQOP_prt", &m_pull_eQOP[0]);
155  m_outputTree->Branch("pull_eT_prt", &m_pull_eT[0]);
156  m_outputTree->Branch("g_x_prt", &m_x[0]);
157  m_outputTree->Branch("g_y_prt", &m_y[0]);
158  m_outputTree->Branch("g_z_prt", &m_z[0]);
159  m_outputTree->Branch("px_prt", &m_px[0]);
160  m_outputTree->Branch("py_prt", &m_py[0]);
161  m_outputTree->Branch("pz_prt", &m_pz[0]);
162  m_outputTree->Branch("eta_prt", &m_eta[0]);
163  m_outputTree->Branch("pT_prt", &m_pT[0]);
164 
165  m_outputTree->Branch("nFiltered", &m_nParams[1]);
166  m_outputTree->Branch("filtered", &m_hasParams[1]);
167  m_outputTree->Branch("eLOC0_flt", &m_eLOC0[1]);
168  m_outputTree->Branch("eLOC1_flt", &m_eLOC1[1]);
169  m_outputTree->Branch("ePHI_flt", &m_ePHI[1]);
170  m_outputTree->Branch("eTHETA_flt", &m_eTHETA[1]);
171  m_outputTree->Branch("eQOP_flt", &m_eQOP[1]);
172  m_outputTree->Branch("eT_flt", &m_eT[1]);
173  m_outputTree->Branch("res_eLOC0_flt", &m_res_eLOC0[1]);
174  m_outputTree->Branch("res_eLOC1_flt", &m_res_eLOC1[1]);
175  m_outputTree->Branch("res_ePHI_flt", &m_res_ePHI[1]);
176  m_outputTree->Branch("res_eTHETA_flt", &m_res_eTHETA[1]);
177  m_outputTree->Branch("res_eQOP_flt", &m_res_eQOP[1]);
178  m_outputTree->Branch("res_eT_flt", &m_res_eT[1]);
179  m_outputTree->Branch("err_eLOC0_flt", &m_err_eLOC0[1]);
180  m_outputTree->Branch("err_eLOC1_flt", &m_err_eLOC1[1]);
181  m_outputTree->Branch("err_ePHI_flt", &m_err_ePHI[1]);
182  m_outputTree->Branch("err_eTHETA_flt", &m_err_eTHETA[1]);
183  m_outputTree->Branch("err_eQOP_flt", &m_err_eQOP[1]);
184  m_outputTree->Branch("err_eT_flt", &m_err_eT[1]);
185  m_outputTree->Branch("pull_eLOC0_flt", &m_pull_eLOC0[1]);
186  m_outputTree->Branch("pull_eLOC1_flt", &m_pull_eLOC1[1]);
187  m_outputTree->Branch("pull_ePHI_flt", &m_pull_ePHI[1]);
188  m_outputTree->Branch("pull_eTHETA_flt", &m_pull_eTHETA[1]);
189  m_outputTree->Branch("pull_eQOP_flt", &m_pull_eQOP[1]);
190  m_outputTree->Branch("pull_eT_flt", &m_pull_eT[1]);
191  m_outputTree->Branch("g_x_flt", &m_x[1]);
192  m_outputTree->Branch("g_y_flt", &m_y[1]);
193  m_outputTree->Branch("g_z_flt", &m_z[1]);
194  m_outputTree->Branch("px_flt", &m_px[1]);
195  m_outputTree->Branch("py_flt", &m_py[1]);
196  m_outputTree->Branch("pz_flt", &m_pz[1]);
197  m_outputTree->Branch("eta_flt", &m_eta[1]);
198  m_outputTree->Branch("pT_flt", &m_pT[1]);
199 
200  m_outputTree->Branch("nSmoothed", &m_nParams[2]);
201  m_outputTree->Branch("smoothed", &m_hasParams[2]);
202  m_outputTree->Branch("eLOC0_smt", &m_eLOC0[2]);
203  m_outputTree->Branch("eLOC1_smt", &m_eLOC1[2]);
204  m_outputTree->Branch("ePHI_smt", &m_ePHI[2]);
205  m_outputTree->Branch("eTHETA_smt", &m_eTHETA[2]);
206  m_outputTree->Branch("eQOP_smt", &m_eQOP[2]);
207  m_outputTree->Branch("eT_smt", &m_eT[2]);
208  m_outputTree->Branch("res_eLOC0_smt", &m_res_eLOC0[2]);
209  m_outputTree->Branch("res_eLOC1_smt", &m_res_eLOC1[2]);
210  m_outputTree->Branch("res_ePHI_smt", &m_res_ePHI[2]);
211  m_outputTree->Branch("res_eTHETA_smt", &m_res_eTHETA[2]);
212  m_outputTree->Branch("res_eQOP_smt", &m_res_eQOP[2]);
213  m_outputTree->Branch("res_eT_smt", &m_res_eT[2]);
214  m_outputTree->Branch("err_eLOC0_smt", &m_err_eLOC0[2]);
215  m_outputTree->Branch("err_eLOC1_smt", &m_err_eLOC1[2]);
216  m_outputTree->Branch("err_ePHI_smt", &m_err_ePHI[2]);
217  m_outputTree->Branch("err_eTHETA_smt", &m_err_eTHETA[2]);
218  m_outputTree->Branch("err_eQOP_smt", &m_err_eQOP[2]);
219  m_outputTree->Branch("err_eT_smt", &m_err_eT[2]);
220  m_outputTree->Branch("pull_eLOC0_smt", &m_pull_eLOC0[2]);
221  m_outputTree->Branch("pull_eLOC1_smt", &m_pull_eLOC1[2]);
222  m_outputTree->Branch("pull_ePHI_smt", &m_pull_ePHI[2]);
223  m_outputTree->Branch("pull_eTHETA_smt", &m_pull_eTHETA[2]);
224  m_outputTree->Branch("pull_eQOP_smt", &m_pull_eQOP[2]);
225  m_outputTree->Branch("pull_eT_smt", &m_pull_eT[2]);
226  m_outputTree->Branch("g_x_smt", &m_x[2]);
227  m_outputTree->Branch("g_y_smt", &m_y[2]);
228  m_outputTree->Branch("g_z_smt", &m_z[2]);
229  m_outputTree->Branch("px_smt", &m_px[2]);
230  m_outputTree->Branch("py_smt", &m_py[2]);
231  m_outputTree->Branch("pz_smt", &m_pz[2]);
232  m_outputTree->Branch("eta_smt", &m_eta[2]);
233  m_outputTree->Branch("pT_smt", &m_pT[2]);
234 
235  m_outputTree->Branch("nUnbiased", &m_nParams[3]);
236  m_outputTree->Branch("unbiased", &m_hasParams[3]);
237  m_outputTree->Branch("eLOC0_ubs", &m_eLOC0[3]);
238  m_outputTree->Branch("eLOC1_ubs", &m_eLOC1[3]);
239  m_outputTree->Branch("ePHI_ubs", &m_ePHI[3]);
240  m_outputTree->Branch("eTHETA_ubs", &m_eTHETA[3]);
241  m_outputTree->Branch("eQOP_ubs", &m_eQOP[3]);
242  m_outputTree->Branch("eT_ubs", &m_eT[3]);
243  m_outputTree->Branch("res_eLOC0_ubs", &m_res_eLOC0[3]);
244  m_outputTree->Branch("res_eLOC1_ubs", &m_res_eLOC1[3]);
245  m_outputTree->Branch("res_ePHI_ubs", &m_res_ePHI[3]);
246  m_outputTree->Branch("res_eTHETA_ubs", &m_res_eTHETA[3]);
247  m_outputTree->Branch("res_eQOP_ubs", &m_res_eQOP[3]);
248  m_outputTree->Branch("res_eT_ubs", &m_res_eT[3]);
249  m_outputTree->Branch("err_eLOC0_ubs", &m_err_eLOC0[3]);
250  m_outputTree->Branch("err_eLOC1_ubs", &m_err_eLOC1[3]);
251  m_outputTree->Branch("err_ePHI_ubs", &m_err_ePHI[3]);
252  m_outputTree->Branch("err_eTHETA_ubs", &m_err_eTHETA[3]);
253  m_outputTree->Branch("err_eQOP_ubs", &m_err_eQOP[3]);
254  m_outputTree->Branch("err_eT_ubs", &m_err_eT[3]);
255  m_outputTree->Branch("pull_eLOC0_ubs", &m_pull_eLOC0[3]);
256  m_outputTree->Branch("pull_eLOC1_ubs", &m_pull_eLOC1[3]);
257  m_outputTree->Branch("pull_ePHI_ubs", &m_pull_ePHI[3]);
258  m_outputTree->Branch("pull_eTHETA_ubs", &m_pull_eTHETA[3]);
259  m_outputTree->Branch("pull_eQOP_ubs", &m_pull_eQOP[3]);
260  m_outputTree->Branch("pull_eT_ubs", &m_pull_eT[3]);
261  m_outputTree->Branch("g_x_ubs", &m_x[3]);
262  m_outputTree->Branch("g_y_ubs", &m_y[3]);
263  m_outputTree->Branch("g_z_ubs", &m_z[3]);
264  m_outputTree->Branch("px_ubs", &m_px[3]);
265  m_outputTree->Branch("py_ubs", &m_py[3]);
266  m_outputTree->Branch("pz_ubs", &m_pz[3]);
267  m_outputTree->Branch("eta_ubs", &m_eta[3]);
268  m_outputTree->Branch("pT_ubs", &m_pT[3]);
269 
270  m_outputTree->Branch("chi2", &m_chi2);
271  }
272 }
273 
274 ActsExamples::RootTrajectoryStatesWriter::~RootTrajectoryStatesWriter() {
275  m_outputFile->Close();
276 }
277 
278 ActsExamples::ProcessCode ActsExamples::RootTrajectoryStatesWriter::finalize() {
279  m_outputFile->cd();
280  m_outputTree->Write();
281  m_outputFile->Close();
282 
283  ACTS_INFO("Wrote states of trajectories to tree '"
284  << m_cfg.treeName << "' in '" << m_cfg.treeName << "'");
285 
286  return ProcessCode::SUCCESS;
287 }
288 
289 ActsExamples::ProcessCode ActsExamples::RootTrajectoryStatesWriter::writeT(
290  const AlgorithmContext& ctx, const TrajectoriesContainer& trajectories) {
291  auto& gctx = ctx.geoContext;
292  // Read additional input collections
293  const auto& particles = m_inputParticles(ctx);
294  const auto& simHits = m_inputSimHits(ctx);
295  const auto& hitParticlesMap = m_inputMeasurementParticlesMap(ctx);
296  const auto& hitSimHitsMap = m_inputMeasurementSimHitsMap(ctx);
297 
298  // For each particle within a track, how many hits did it contribute
299  std::vector<ParticleHitCount> particleHitCounts;
300 
301  // Exclusive access to the tree while writing
302  std::lock_guard<std::mutex> lock(m_writeMutex);
303 
304  // Get the event number
305  m_eventNr = ctx.eventNumber;
306 
307  // Loop over the trajectories
308  for (size_t itraj = 0; itraj < trajectories.size(); ++itraj) {
309  const auto& traj = trajectories[itraj];
310 
311  // The trajectory index
312  m_multiTrajNr = itraj;
313 
314  // The trajectory entry indices
315  const auto& trackTips = traj.tips();
316 
317  // Don't write empty MultiTrajectory
318  if (trackTips.empty()) {
319  continue;
320  }
321 
322  // The MultiTrajectory
323  const auto& mj = traj.multiTrajectory();
324 
325  // Loop over the entry indices for the subtrajectories
326  for (unsigned int isubtraj = 0; isubtraj < trackTips.size(); ++isubtraj) {
327  // The subtrajectory index
328  m_subTrajNr = isubtraj;
329  // The entry index for this subtrajectory
330  const auto& trackTip = trackTips[isubtraj];
331  // Collect the trajectory summary info
332  auto trajState =
333  Acts::MultiTrajectoryHelpers::trajectoryState(mj, trackTip);
334  m_nMeasurements = trajState.nMeasurements;
335  m_nStates = trajState.nStates;
336 
337  // Get the majority truth particle to this track
338  int truthQ = 1.;
339  identifyContributingParticles(hitParticlesMap, traj, trackTip,
340  particleHitCounts);
341  if (not particleHitCounts.empty()) {
342  // Get the barcode of the majority truth particle
343  auto barcode = particleHitCounts.front().particleId;
344  // Find the truth particle via the barcode
345  auto ip = particles.find(barcode);
346  if (ip != particles.end()) {
347  const auto& particle = *ip;
348  ACTS_VERBOSE("Find the truth particle with barcode "
349  << barcode << "=" << barcode.value());
350  // Get the truth particle charge
351  truthQ = static_cast<int>(particle.charge());
352  } else {
353  ACTS_DEBUG("Truth particle with barcode "
354  << barcode << "=" << barcode.value() << " not found!");
355  }
356  }
357 
358  // Get the trackStates on the trajectory
359  m_nParams = {0, 0, 0};
360  mj.visitBackwards(trackTip, [&](const auto& state) {
361  // we only fill the track states with non-outlier measurement
362  auto typeFlags = state.typeFlags();
363  if (not typeFlags.test(Acts::TrackStateFlag::MeasurementFlag)) {
364  return true;
365  }
366 
367  const auto& surface = state.referenceSurface();
368 
369  // get the truth hits corresponding to this trackState
370  // Use average truth in the case of multiple contributing sim hits
371  auto sl =
372  state.getUncalibratedSourceLink().template get<IndexSourceLink>();
373  const auto hitIdx = sl.index();
374  auto indices = makeRange(hitSimHitsMap.equal_range(hitIdx));
375  auto [truthLocal, truthPos4, truthUnitDir] =
376  averageSimHits(ctx.geoContext, surface, simHits, indices, logger());
377  // momentum averaging makes even less sense than averaging position and
378  // direction. use the first momentum or set q/p to zero
379  float truthQOP = 0.0f;
380  if (not indices.empty()) {
381  // we assume that the indices are within valid ranges so we do not
382  // need to check their validity again.
383  const auto simHitIdx0 = indices.begin()->second;
384  const auto& simHit0 = *simHits.nth(simHitIdx0);
385  const auto p =
386  simHit0.momentum4Before().template segment<3>(Acts::eMom0).norm();
387  truthQOP = truthQ / p;
388  }
389 
390  // fill the truth hit info
391  m_t_x.push_back(truthPos4[Acts::ePos0]);
392  m_t_y.push_back(truthPos4[Acts::ePos1]);
393  m_t_z.push_back(truthPos4[Acts::ePos2]);
394  m_t_r.push_back(perp(truthPos4.template segment<3>(Acts::ePos0)));
395  m_t_dx.push_back(truthUnitDir[Acts::eMom0]);
396  m_t_dy.push_back(truthUnitDir[Acts::eMom1]);
397  m_t_dz.push_back(truthUnitDir[Acts::eMom2]);
398 
399  // get the truth track parameter at this track State
400  float truthLOC0 = truthLocal[Acts::ePos0];
401  float truthLOC1 = truthLocal[Acts::ePos1];
402  float truthTIME = truthPos4[Acts::eTime];
403  float truthPHI = phi(truthUnitDir);
404  float truthTHETA = theta(truthUnitDir);
405 
406  // fill the truth track parameter at this track State
407  m_t_eLOC0.push_back(truthLOC0);
408  m_t_eLOC1.push_back(truthLOC1);
409  m_t_ePHI.push_back(truthPHI);
410  m_t_eTHETA.push_back(truthTHETA);
411  m_t_eQOP.push_back(truthQOP);
412  m_t_eT.push_back(truthTIME);
413 
414  // get the geometry ID
415  auto geoID = surface.geometryId();
416  m_volumeID.push_back(geoID.volume());
417  m_layerID.push_back(geoID.layer());
418  m_moduleID.push_back(geoID.sensitive());
419 
420  // get the path length
421  m_pathLength.push_back(state.pathLength());
422 
423  // expand the local measurements into the full bound space
424  Acts::BoundVector meas = state.effectiveProjector().transpose() *
425  state.effectiveCalibrated();
426  // extract local and global position
427  Acts::Vector2 local(meas[Acts::eBoundLoc0], meas[Acts::eBoundLoc1]);
428  Acts::Vector3 mom(1, 1, 1);
429  Acts::Vector3 global =
430  surface.localToGlobal(ctx.geoContext, local, mom);
431 
432  // fill the measurement info
433  m_lx_hit.push_back(local[Acts::ePos0]);
434  m_ly_hit.push_back(local[Acts::ePos1]);
435  m_x_hit.push_back(global[Acts::ePos0]);
436  m_y_hit.push_back(global[Acts::ePos1]);
437  m_z_hit.push_back(global[Acts::ePos2]);
438 
439  // status of the fitted track parameters
440  std::array<bool, 4> hasParams = {false, false, false, false};
441  // optional fitted track parameters
442  std::optional<std::pair<Acts::BoundVector, Acts::BoundMatrix>>
443  trackParamsOpt = std::nullopt;
444  // lambda to get the fitted track parameters
445  auto getTrackParams = [&](unsigned int ipar) {
446  if (ipar == 0 && state.hasPredicted()) {
447  hasParams[0] = true;
448  m_nParams[0]++;
449  trackParamsOpt =
450  std::make_pair(state.predicted(), state.predictedCovariance());
451  } else if (ipar == 1 && state.hasFiltered()) {
452  hasParams[1] = true;
453  m_nParams[1]++;
454  trackParamsOpt =
455  std::make_pair(state.filtered(), state.filteredCovariance());
456  } else if (ipar == 2 && state.hasSmoothed()) {
457  hasParams[2] = true;
458  m_nParams[2]++;
459  trackParamsOpt =
460  std::make_pair(state.smoothed(), state.smoothedCovariance());
461  } else if (ipar == 3 && state.hasSmoothed()) {
462  hasParams[3] = true;
463  m_nParams[3]++;
464  // calculate the unbiased track parameters (i.e. fitted track
465  // parameters with this measurement removed) using Eq.(12a)-Eq.(12c)
466  // of NIMA 262, 444 (1987)
467  auto m = state.effectiveCalibrated();
468  auto H = state.effectiveProjector();
469  auto V = state.effectiveCalibratedCovariance();
470  auto K =
471  (state.smoothedCovariance() * H.transpose() *
472  (H * state.smoothedCovariance() * H.transpose() - V).inverse())
473  .eval();
474  auto unbiasedParamsVec =
475  state.smoothed() + K * (m - H * state.smoothed());
476  auto unbiasedParamsCov =
477  state.smoothedCovariance() - K * H * state.smoothedCovariance();
478  trackParamsOpt =
479  std::make_pair(unbiasedParamsVec, unbiasedParamsCov);
480  }
481  };
482 
483  // fill the fitted track parameters
484  for (unsigned int ipar = 0; ipar < 4; ++ipar) {
485  // get the fitted track parameters
486  getTrackParams(ipar);
487  if (trackParamsOpt) {
488  const auto& [parameters, covariance] = *trackParamsOpt;
489  if (ipar == 0) {
490  //
491  // local hit residual info
492  auto H = state.effectiveProjector();
493  auto V = state.effectiveCalibratedCovariance();
494  auto resCov = V + H * covariance * H.transpose();
495  Acts::ActsDynamicVector res(state.calibratedSize());
496  res.setZero();
497 
498  res = state.effectiveCalibrated() - H * parameters;
499 
500  m_res_x_hit.push_back(res[Acts::eBoundLoc0]);
501  m_err_x_hit.push_back(
502  std::sqrt(V(Acts::eBoundLoc0, Acts::eBoundLoc0)));
503  m_pull_x_hit.push_back(
504  res[Acts::eBoundLoc0] /
505  std::sqrt(resCov(Acts::eBoundLoc0, Acts::eBoundLoc0)));
506 
507  if (state.calibratedSize() >= 2) {
508  m_res_y_hit.push_back(res[Acts::eBoundLoc1]);
509  m_err_y_hit.push_back(
510  std::sqrt(V(Acts::eBoundLoc1, Acts::eBoundLoc1)));
511  m_pull_y_hit.push_back(
512  res[Acts::eBoundLoc1] /
513  std::sqrt(resCov(Acts::eBoundLoc1, Acts::eBoundLoc1)));
514  } else {
515  float nan = std::numeric_limits<float>::quiet_NaN();
516  m_res_y_hit.push_back(nan);
517  m_err_y_hit.push_back(nan);
518  m_pull_y_hit.push_back(nan);
519  }
520 
521  m_dim_hit.push_back(state.calibratedSize());
522  }
523 
524  // track parameters
525  m_eLOC0[ipar].push_back(parameters[Acts::eBoundLoc0]);
526  m_eLOC1[ipar].push_back(parameters[Acts::eBoundLoc1]);
527  m_ePHI[ipar].push_back(parameters[Acts::eBoundPhi]);
528  m_eTHETA[ipar].push_back(parameters[Acts::eBoundTheta]);
529  m_eQOP[ipar].push_back(parameters[Acts::eBoundQOverP]);
530  m_eT[ipar].push_back(parameters[Acts::eBoundTime]);
531 
532  // track parameters residual
533  m_res_eLOC0[ipar].push_back(parameters[Acts::eBoundLoc0] -
534  truthLOC0);
535  m_res_eLOC1[ipar].push_back(parameters[Acts::eBoundLoc1] -
536  truthLOC1);
537  float resPhi = Acts::detail::difference_periodic<float>(
538  parameters[Acts::eBoundPhi], truthPHI,
539  static_cast<float>(2 * M_PI));
540  m_res_ePHI[ipar].push_back(resPhi);
541  m_res_eTHETA[ipar].push_back(parameters[Acts::eBoundTheta] -
542  truthTHETA);
543  m_res_eQOP[ipar].push_back(parameters[Acts::eBoundQOverP] -
544  truthQOP);
545  m_res_eT[ipar].push_back(parameters[Acts::eBoundTime] - truthTIME);
546 
547  // track parameters error
548  // MARK: fpeMaskBegin(FLTINV, 1, #2348)
549  m_err_eLOC0[ipar].push_back(
550  std::sqrt(covariance(Acts::eBoundLoc0, Acts::eBoundLoc0)));
551  m_err_eLOC1[ipar].push_back(
552  std::sqrt(covariance(Acts::eBoundLoc1, Acts::eBoundLoc1)));
553  m_err_ePHI[ipar].push_back(
554  std::sqrt(covariance(Acts::eBoundPhi, Acts::eBoundPhi)));
555  m_err_eTHETA[ipar].push_back(
556  std::sqrt(covariance(Acts::eBoundTheta, Acts::eBoundTheta)));
557  m_err_eQOP[ipar].push_back(
558  std::sqrt(covariance(Acts::eBoundQOverP, Acts::eBoundQOverP)));
559  m_err_eT[ipar].push_back(
560  std::sqrt(covariance(Acts::eBoundTime, Acts::eBoundTime)));
561  // MARK: fpeMaskEnd(FLTINV)
562 
563  // track parameters pull
564  m_pull_eLOC0[ipar].push_back(
565  (parameters[Acts::eBoundLoc0] - truthLOC0) /
566  std::sqrt(covariance(Acts::eBoundLoc0, Acts::eBoundLoc0)));
567  m_pull_eLOC1[ipar].push_back(
568  (parameters[Acts::eBoundLoc1] - truthLOC1) /
569  std::sqrt(covariance(Acts::eBoundLoc1, Acts::eBoundLoc1)));
570  m_pull_ePHI[ipar].push_back(
571  resPhi /
572  std::sqrt(covariance(Acts::eBoundPhi, Acts::eBoundPhi)));
573  m_pull_eTHETA[ipar].push_back(
574  (parameters[Acts::eBoundTheta] - truthTHETA) /
575  std::sqrt(covariance(Acts::eBoundTheta, Acts::eBoundTheta)));
576  m_pull_eQOP[ipar].push_back(
577  (parameters[Acts::eBoundQOverP] - truthQOP) /
578  std::sqrt(covariance(Acts::eBoundQOverP, Acts::eBoundQOverP)));
579  double sigmaTime =
580  std::sqrt(covariance(Acts::eBoundTime, Acts::eBoundTime));
581  m_pull_eT[ipar].push_back(
582  sigmaTime == 0.0
583  ? std::numeric_limits<double>::quiet_NaN()
584  : (parameters[Acts::eBoundTime] - truthTIME) / sigmaTime);
585 
586  // further track parameter info
587  Acts::FreeVector freeParams =
588  Acts::detail::transformBoundToFreeParameters(surface, gctx,
589  parameters);
590  m_x[ipar].push_back(freeParams[Acts::eFreePos0]);
591  m_y[ipar].push_back(freeParams[Acts::eFreePos1]);
592  m_z[ipar].push_back(freeParams[Acts::eFreePos2]);
593  auto p = std::abs(1 / freeParams[Acts::eFreeQOverP]);
594  m_px[ipar].push_back(p * freeParams[Acts::eFreeDir0]);
595  m_py[ipar].push_back(p * freeParams[Acts::eFreeDir1]);
596  m_pz[ipar].push_back(p * freeParams[Acts::eFreeDir2]);
597  m_pT[ipar].push_back(p * std::hypot(freeParams[Acts::eFreeDir0],
598  freeParams[Acts::eFreeDir1]));
599  m_eta[ipar].push_back(Acts::VectorHelpers::eta(
600  freeParams.segment<3>(Acts::eFreeDir0)));
601  } else {
602  if (ipar == 0) {
603  // push default values if no track parameters
604  m_res_x_hit.push_back(-99.);
605  m_res_y_hit.push_back(-99.);
606  m_err_x_hit.push_back(-99.);
607  m_err_y_hit.push_back(-99.);
608  m_pull_x_hit.push_back(-99.);
609  m_pull_y_hit.push_back(-99.);
610  m_dim_hit.push_back(-99.);
611  }
612  // push default values if no track parameters
613  m_eLOC0[ipar].push_back(-99.);
614  m_eLOC1[ipar].push_back(-99.);
615  m_ePHI[ipar].push_back(-99.);
616  m_eTHETA[ipar].push_back(-99.);
617  m_eQOP[ipar].push_back(-99.);
618  m_eT[ipar].push_back(-99.);
619  m_res_eLOC0[ipar].push_back(-99.);
620  m_res_eLOC1[ipar].push_back(-99.);
621  m_res_ePHI[ipar].push_back(-99.);
622  m_res_eTHETA[ipar].push_back(-99.);
623  m_res_eQOP[ipar].push_back(-99.);
624  m_res_eT[ipar].push_back(-99.);
625  m_err_eLOC0[ipar].push_back(-99);
626  m_err_eLOC1[ipar].push_back(-99);
627  m_err_ePHI[ipar].push_back(-99);
628  m_err_eTHETA[ipar].push_back(-99);
629  m_err_eQOP[ipar].push_back(-99);
630  m_err_eT[ipar].push_back(-99);
631  m_pull_eLOC0[ipar].push_back(-99.);
632  m_pull_eLOC1[ipar].push_back(-99.);
633  m_pull_ePHI[ipar].push_back(-99.);
634  m_pull_eTHETA[ipar].push_back(-99.);
635  m_pull_eQOP[ipar].push_back(-99.);
636  m_pull_eT[ipar].push_back(-99.);
637  m_x[ipar].push_back(-99.);
638  m_y[ipar].push_back(-99.);
639  m_z[ipar].push_back(-99.);
640  m_px[ipar].push_back(-99.);
641  m_py[ipar].push_back(-99.);
642  m_pz[ipar].push_back(-99.);
643  m_pT[ipar].push_back(-99.);
644  m_eta[ipar].push_back(-99.);
645  }
646  // fill the track parameters status
647  m_hasParams[ipar].push_back(hasParams[ipar]);
648  }
649 
650  // fill the chi2
651  m_chi2.push_back(state.chi2());
652 
653  return true;
654  }); // all states
655 
656  // fill the variables for one track to tree
657  m_outputTree->Fill();
658 
659  // now reset
660  m_t_x.clear();
661  m_t_y.clear();
662  m_t_z.clear();
663  m_t_r.clear();
664  m_t_dx.clear();
665  m_t_dy.clear();
666  m_t_dz.clear();
667  m_t_eLOC0.clear();
668  m_t_eLOC1.clear();
669  m_t_ePHI.clear();
670  m_t_eTHETA.clear();
671  m_t_eQOP.clear();
672  m_t_eT.clear();
673 
674  m_volumeID.clear();
675  m_layerID.clear();
676  m_moduleID.clear();
677  m_pathLength.clear();
678  m_lx_hit.clear();
679  m_ly_hit.clear();
680  m_x_hit.clear();
681  m_y_hit.clear();
682  m_z_hit.clear();
683  m_res_x_hit.clear();
684  m_res_y_hit.clear();
685  m_err_x_hit.clear();
686  m_err_y_hit.clear();
687  m_pull_x_hit.clear();
688  m_pull_y_hit.clear();
689  m_dim_hit.clear();
690 
691  for (unsigned int ipar = 0; ipar < 4; ++ipar) {
692  m_hasParams[ipar].clear();
693  m_eLOC0[ipar].clear();
694  m_eLOC1[ipar].clear();
695  m_ePHI[ipar].clear();
696  m_eTHETA[ipar].clear();
697  m_eQOP[ipar].clear();
698  m_eT[ipar].clear();
699  m_res_eLOC0[ipar].clear();
700  m_res_eLOC1[ipar].clear();
701  m_res_ePHI[ipar].clear();
702  m_res_eTHETA[ipar].clear();
703  m_res_eQOP[ipar].clear();
704  m_res_eT[ipar].clear();
705  m_err_eLOC0[ipar].clear();
706  m_err_eLOC1[ipar].clear();
707  m_err_ePHI[ipar].clear();
708  m_err_eTHETA[ipar].clear();
709  m_err_eQOP[ipar].clear();
710  m_err_eT[ipar].clear();
711  m_pull_eLOC0[ipar].clear();
712  m_pull_eLOC1[ipar].clear();
713  m_pull_ePHI[ipar].clear();
714  m_pull_eTHETA[ipar].clear();
715  m_pull_eQOP[ipar].clear();
716  m_pull_eT[ipar].clear();
717  m_x[ipar].clear();
718  m_y[ipar].clear();
719  m_z[ipar].clear();
720  m_px[ipar].clear();
721  m_py[ipar].clear();
722  m_pz[ipar].clear();
723  m_eta[ipar].clear();
724  m_pT[ipar].clear();
725  }
726 
727  m_chi2.clear();
728  } // all subtrajectories
729  } // all trajectories
730 
731  return ProcessCode::SUCCESS;
732 }