Root/src/RootTrackStatesWriter.cpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 #include "ActsExamples/Io/Root/RootTrackStatesWriter.hpp"
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Definitions/Common.hpp"
0013 #include "Acts/Definitions/TrackParametrization.hpp"
0014 #include "Acts/EventData/MultiTrajectory.hpp"
0015 #include "Acts/EventData/TransformationHelpers.hpp"
0016 #include "Acts/EventData/VectorMultiTrajectory.hpp"
0017 #include "Acts/Geometry/GeometryContext.hpp"
0018 #include "Acts/Geometry/GeometryIdentifier.hpp"
0019 #include "Acts/Utilities/Helpers.hpp"
0020 #include "Acts/Utilities/TrackHelpers.hpp"
0021 #include "Acts/Utilities/detail/periodic.hpp"
0022 #include "ActsExamples/EventData/AverageSimHits.hpp"
0023 #include "ActsExamples/EventData/IndexSourceLink.hpp"
0024 #include "ActsExamples/EventData/Track.hpp"
0025 #include "ActsExamples/Framework/AlgorithmContext.hpp"
0026 #include "ActsExamples/Utilities/Range.hpp"
0027 #include "ActsFatras/EventData/Barcode.hpp"
0028
0029 #include <cmath>
0030 #include <ios>
0031 #include <limits>
0032 #include <numbers>
0033 #include <optional>
0034 #include <ostream>
0035 #include <stdexcept>
0036 #include <utility>
0037
0038 #include <TFile.h>
0039 #include <TTree.h>
0040
0041 namespace ActsExamples {
0042
0043 using Acts::VectorHelpers::eta;
0044 using Acts::VectorHelpers::perp;
0045 using Acts::VectorHelpers::phi;
0046 using Acts::VectorHelpers::theta;
0047
0048 RootTrackStatesWriter::RootTrackStatesWriter(
0049     const RootTrackStatesWriter::Config& config, Acts::Logging::Level level)
0050     : WriterT(config.inputTracks, "RootTrackStatesWriter", level),
0051       m_cfg(config) {
0052   // trajectories collection name is already checked by base ctor
0053   if (m_cfg.inputParticles.empty()) {
0054     throw std::invalid_argument("Missing particles input collection");
0055   }
0056   if (m_cfg.inputTrackParticleMatching.empty()) {
0057     throw std::invalid_argument("Missing input track particles matching");
0058   }
0059   if (m_cfg.inputSimHits.empty()) {
0060     throw std::invalid_argument("Missing simulated hits input collection");
0061   }
0062   if (m_cfg.inputMeasurementSimHitsMap.empty()) {
0063     throw std::invalid_argument(
0064         "Missing hit-simulated-hits map input collection");
0065   }
0066   if (m_cfg.filePath.empty()) {
0067     throw std::invalid_argument("Missing output filename");
0068   }
0069   if (m_cfg.treeName.empty()) {
0070     throw std::invalid_argument("Missing tree name");
0071   }
0072
0073   m_inputParticles.initialize(m_cfg.inputParticles);
0074   m_inputTrackParticleMatching.initialize(m_cfg.inputTrackParticleMatching);
0075   m_inputSimHits.initialize(m_cfg.inputSimHits);
0076   m_inputMeasurementSimHitsMap.initialize(m_cfg.inputMeasurementSimHitsMap);
0077
0078   // Setup ROOT I/O
0079   auto path = m_cfg.filePath;
0080   m_outputFile = TFile::Open(path.c_str(), m_cfg.fileMode.c_str());
0081   if (m_outputFile == nullptr) {
0082     throw std::ios_base::failure("Could not open '" + path + "'");
0083   }
0084   m_outputFile->cd();
0085   m_outputTree = new TTree(m_cfg.treeName.c_str(), m_cfg.treeName.c_str());
0086   if (m_outputTree == nullptr) {
0087     throw std::bad_alloc();
0088   }
0089
0090   // I/O parameters
0091   m_outputTree->Branch("event_nr", &m_eventNr);
0092   m_outputTree->Branch("track_nr", &m_trackNr);
0093
0094   m_outputTree->Branch("nStates", &m_nStates);
0095   m_outputTree->Branch("nMeasurements", &m_nMeasurements);
0096
0097   m_outputTree->Branch("volume_id", &m_volumeID);
0098   m_outputTree->Branch("layer_id", &m_layerID);
0099   m_outputTree->Branch("module_id", &m_moduleID);
0100
0101   m_outputTree->Branch("stateType", &m_stateType);
0102
0103   m_outputTree->Branch("chi2", &m_chi2);
0104
0105   m_outputTree->Branch("pathLength", &m_pathLength);
0106
0107   m_outputTree->Branch("t_x", &m_t_x);
0108   m_outputTree->Branch("t_y", &m_t_y);
0109   m_outputTree->Branch("t_z", &m_t_z);
0110   m_outputTree->Branch("t_r", &m_t_r);
0111   m_outputTree->Branch("t_dx", &m_t_dx);
0112   m_outputTree->Branch("t_dy", &m_t_dy);
0113   m_outputTree->Branch("t_dz", &m_t_dz);
0114   m_outputTree->Branch("t_eLOC0", &m_t_eLOC0);
0115   m_outputTree->Branch("t_eLOC1", &m_t_eLOC1);
0116   m_outputTree->Branch("t_ePHI", &m_t_ePHI);
0117   m_outputTree->Branch("t_eTHETA", &m_t_eTHETA);
0118   m_outputTree->Branch("t_eQOP", &m_t_eQOP);
0119   m_outputTree->Branch("t_eT", &m_t_eT);
0120   m_outputTree->Branch("particle_ids_vertex_primary", &m_particleVertexPrimary);
0121   m_outputTree->Branch("particle_ids_vertex_secondary",
0122                        &m_particleVertexSecondary);
0123   m_outputTree->Branch("particle_ids_particle", &m_particleParticle);
0124   m_outputTree->Branch("particle_ids_generation", &m_particleGeneration);
0125   m_outputTree->Branch("particle_ids_sub_particle", &m_particleSubParticle);
0126
0127   m_outputTree->Branch("dim_hit", &m_dim_hit);
0128   m_outputTree->Branch("l_x_hit", &m_lx_hit);
0129   m_outputTree->Branch("l_y_hit", &m_ly_hit);
0130   m_outputTree->Branch("g_x_hit", &m_x_hit);
0131   m_outputTree->Branch("g_y_hit", &m_y_hit);
0132   m_outputTree->Branch("g_z_hit", &m_z_hit);
0133   m_outputTree->Branch("res_x_hit", &m_res_x_hit);
0134   m_outputTree->Branch("res_y_hit", &m_res_y_hit);
0135   m_outputTree->Branch("err_x_hit", &m_err_x_hit);
0136   m_outputTree->Branch("err_y_hit", &m_err_y_hit);
0137   m_outputTree->Branch("pull_x_hit", &m_pull_x_hit);
0138   m_outputTree->Branch("pull_y_hit", &m_pull_y_hit);
0139
0140   m_outputTree->Branch("nPredicted", &m_nParams[ePredicted]);
0141   m_outputTree->Branch("predicted", &m_hasParams[ePredicted]);
0142   m_outputTree->Branch("eLOC0_prt", &m_eLOC0[ePredicted]);
0143   m_outputTree->Branch("eLOC1_prt", &m_eLOC1[ePredicted]);
0144   m_outputTree->Branch("ePHI_prt", &m_ePHI[ePredicted]);
0145   m_outputTree->Branch("eTHETA_prt", &m_eTHETA[ePredicted]);
0146   m_outputTree->Branch("eQOP_prt", &m_eQOP[ePredicted]);
0147   m_outputTree->Branch("eT_prt", &m_eT[ePredicted]);
0148   m_outputTree->Branch("res_eLOC0_prt", &m_res_eLOC0[ePredicted]);
0149   m_outputTree->Branch("res_eLOC1_prt", &m_res_eLOC1[ePredicted]);
0150   m_outputTree->Branch("res_ePHI_prt", &m_res_ePHI[ePredicted]);
0151   m_outputTree->Branch("res_eTHETA_prt", &m_res_eTHETA[ePredicted]);
0152   m_outputTree->Branch("res_eQOP_prt", &m_res_eQOP[ePredicted]);
0153   m_outputTree->Branch("res_eT_prt", &m_res_eT[ePredicted]);
0154   m_outputTree->Branch("err_eLOC0_prt", &m_err_eLOC0[ePredicted]);
0155   m_outputTree->Branch("err_eLOC1_prt", &m_err_eLOC1[ePredicted]);
0156   m_outputTree->Branch("err_ePHI_prt", &m_err_ePHI[ePredicted]);
0157   m_outputTree->Branch("err_eTHETA_prt", &m_err_eTHETA[ePredicted]);
0158   m_outputTree->Branch("err_eQOP_prt", &m_err_eQOP[ePredicted]);
0159   m_outputTree->Branch("err_eT_prt", &m_err_eT[ePredicted]);
0160   m_outputTree->Branch("pull_eLOC0_prt", &m_pull_eLOC0[ePredicted]);
0161   m_outputTree->Branch("pull_eLOC1_prt", &m_pull_eLOC1[ePredicted]);
0162   m_outputTree->Branch("pull_ePHI_prt", &m_pull_ePHI[ePredicted]);
0163   m_outputTree->Branch("pull_eTHETA_prt", &m_pull_eTHETA[ePredicted]);
0164   m_outputTree->Branch("pull_eQOP_prt", &m_pull_eQOP[ePredicted]);
0165   m_outputTree->Branch("pull_eT_prt", &m_pull_eT[ePredicted]);
0166   m_outputTree->Branch("g_x_prt", &m_x[ePredicted]);
0167   m_outputTree->Branch("g_y_prt", &m_y[ePredicted]);
0168   m_outputTree->Branch("g_z_prt", &m_z[ePredicted]);
0169   m_outputTree->Branch("px_prt", &m_px[ePredicted]);
0170   m_outputTree->Branch("py_prt", &m_py[ePredicted]);
0171   m_outputTree->Branch("pz_prt", &m_pz[ePredicted]);
0172   m_outputTree->Branch("eta_prt", &m_eta[ePredicted]);
0173   m_outputTree->Branch("pT_prt", &m_pT[ePredicted]);
0174
0175   m_outputTree->Branch("nFiltered", &m_nParams[eFiltered]);
0176   m_outputTree->Branch("filtered", &m_hasParams[eFiltered]);
0177   m_outputTree->Branch("eLOC0_flt", &m_eLOC0[eFiltered]);
0178   m_outputTree->Branch("eLOC1_flt", &m_eLOC1[eFiltered]);
0179   m_outputTree->Branch("ePHI_flt", &m_ePHI[eFiltered]);
0180   m_outputTree->Branch("eTHETA_flt", &m_eTHETA[eFiltered]);
0181   m_outputTree->Branch("eQOP_flt", &m_eQOP[eFiltered]);
0182   m_outputTree->Branch("eT_flt", &m_eT[eFiltered]);
0183   m_outputTree->Branch("res_eLOC0_flt", &m_res_eLOC0[eFiltered]);
0184   m_outputTree->Branch("res_eLOC1_flt", &m_res_eLOC1[eFiltered]);
0185   m_outputTree->Branch("res_ePHI_flt", &m_res_ePHI[eFiltered]);
0186   m_outputTree->Branch("res_eTHETA_flt", &m_res_eTHETA[eFiltered]);
0187   m_outputTree->Branch("res_eQOP_flt", &m_res_eQOP[eFiltered]);
0188   m_outputTree->Branch("res_eT_flt", &m_res_eT[eFiltered]);
0189   m_outputTree->Branch("err_eLOC0_flt", &m_err_eLOC0[eFiltered]);
0190   m_outputTree->Branch("err_eLOC1_flt", &m_err_eLOC1[eFiltered]);
0191   m_outputTree->Branch("err_ePHI_flt", &m_err_ePHI[eFiltered]);
0192   m_outputTree->Branch("err_eTHETA_flt", &m_err_eTHETA[eFiltered]);
0193   m_outputTree->Branch("err_eQOP_flt", &m_err_eQOP[eFiltered]);
0194   m_outputTree->Branch("err_eT_flt", &m_err_eT[eFiltered]);
0195   m_outputTree->Branch("pull_eLOC0_flt", &m_pull_eLOC0[eFiltered]);
0196   m_outputTree->Branch("pull_eLOC1_flt", &m_pull_eLOC1[eFiltered]);
0197   m_outputTree->Branch("pull_ePHI_flt", &m_pull_ePHI[eFiltered]);
0198   m_outputTree->Branch("pull_eTHETA_flt", &m_pull_eTHETA[eFiltered]);
0199   m_outputTree->Branch("pull_eQOP_flt", &m_pull_eQOP[eFiltered]);
0200   m_outputTree->Branch("pull_eT_flt", &m_pull_eT[eFiltered]);
0201   m_outputTree->Branch("g_x_flt", &m_x[eFiltered]);
0202   m_outputTree->Branch("g_y_flt", &m_y[eFiltered]);
0203   m_outputTree->Branch("g_z_flt", &m_z[eFiltered]);
0204   m_outputTree->Branch("px_flt", &m_px[eFiltered]);
0205   m_outputTree->Branch("py_flt", &m_py[eFiltered]);
0206   m_outputTree->Branch("pz_flt", &m_pz[eFiltered]);
0207   m_outputTree->Branch("eta_flt", &m_eta[eFiltered]);
0208   m_outputTree->Branch("pT_flt", &m_pT[eFiltered]);
0209
0210   m_outputTree->Branch("nSmoothed", &m_nParams[eSmoothed]);
0211   m_outputTree->Branch("smoothed", &m_hasParams[eSmoothed]);
0212   m_outputTree->Branch("eLOC0_smt", &m_eLOC0[eSmoothed]);
0213   m_outputTree->Branch("eLOC1_smt", &m_eLOC1[eSmoothed]);
0214   m_outputTree->Branch("ePHI_smt", &m_ePHI[eSmoothed]);
0215   m_outputTree->Branch("eTHETA_smt", &m_eTHETA[eSmoothed]);
0216   m_outputTree->Branch("eQOP_smt", &m_eQOP[eSmoothed]);
0217   m_outputTree->Branch("eT_smt", &m_eT[eSmoothed]);
0218   m_outputTree->Branch("res_eLOC0_smt", &m_res_eLOC0[eSmoothed]);
0219   m_outputTree->Branch("res_eLOC1_smt", &m_res_eLOC1[eSmoothed]);
0220   m_outputTree->Branch("res_ePHI_smt", &m_res_ePHI[eSmoothed]);
0221   m_outputTree->Branch("res_eTHETA_smt", &m_res_eTHETA[eSmoothed]);
0222   m_outputTree->Branch("res_eQOP_smt", &m_res_eQOP[eSmoothed]);
0223   m_outputTree->Branch("res_eT_smt", &m_res_eT[eSmoothed]);
0224   m_outputTree->Branch("err_eLOC0_smt", &m_err_eLOC0[eSmoothed]);
0225   m_outputTree->Branch("err_eLOC1_smt", &m_err_eLOC1[eSmoothed]);
0226   m_outputTree->Branch("err_ePHI_smt", &m_err_ePHI[eSmoothed]);
0227   m_outputTree->Branch("err_eTHETA_smt", &m_err_eTHETA[eSmoothed]);
0228   m_outputTree->Branch("err_eQOP_smt", &m_err_eQOP[eSmoothed]);
0229   m_outputTree->Branch("err_eT_smt", &m_err_eT[eSmoothed]);
0230   m_outputTree->Branch("pull_eLOC0_smt", &m_pull_eLOC0[eSmoothed]);
0231   m_outputTree->Branch("pull_eLOC1_smt", &m_pull_eLOC1[eSmoothed]);
0232   m_outputTree->Branch("pull_ePHI_smt", &m_pull_ePHI[eSmoothed]);
0233   m_outputTree->Branch("pull_eTHETA_smt", &m_pull_eTHETA[eSmoothed]);
0234   m_outputTree->Branch("pull_eQOP_smt", &m_pull_eQOP[eSmoothed]);
0235   m_outputTree->Branch("pull_eT_smt", &m_pull_eT[eSmoothed]);
0236   m_outputTree->Branch("g_x_smt", &m_x[eSmoothed]);
0237   m_outputTree->Branch("g_y_smt", &m_y[eSmoothed]);
0238   m_outputTree->Branch("g_z_smt", &m_z[eSmoothed]);
0239   m_outputTree->Branch("px_smt", &m_px[eSmoothed]);
0240   m_outputTree->Branch("py_smt", &m_py[eSmoothed]);
0241   m_outputTree->Branch("pz_smt", &m_pz[eSmoothed]);
0242   m_outputTree->Branch("eta_smt", &m_eta[eSmoothed]);
0243   m_outputTree->Branch("pT_smt", &m_pT[eSmoothed]);
0244
0245   m_outputTree->Branch("nUnbiased", &m_nParams[eUnbiased]);
0246   m_outputTree->Branch("unbiased", &m_hasParams[eUnbiased]);
0247   m_outputTree->Branch("eLOC0_ubs", &m_eLOC0[eUnbiased]);
0248   m_outputTree->Branch("eLOC1_ubs", &m_eLOC1[eUnbiased]);
0249   m_outputTree->Branch("ePHI_ubs", &m_ePHI[eUnbiased]);
0250   m_outputTree->Branch("eTHETA_ubs", &m_eTHETA[eUnbiased]);
0251   m_outputTree->Branch("eQOP_ubs", &m_eQOP[eUnbiased]);
0252   m_outputTree->Branch("eT_ubs", &m_eT[eUnbiased]);
0253   m_outputTree->Branch("res_eLOC0_ubs", &m_res_eLOC0[eUnbiased]);
0254   m_outputTree->Branch("res_eLOC1_ubs", &m_res_eLOC1[eUnbiased]);
0255   m_outputTree->Branch("res_ePHI_ubs", &m_res_ePHI[eUnbiased]);
0256   m_outputTree->Branch("res_eTHETA_ubs", &m_res_eTHETA[eUnbiased]);
0257   m_outputTree->Branch("res_eQOP_ubs", &m_res_eQOP[eUnbiased]);
0258   m_outputTree->Branch("res_eT_ubs", &m_res_eT[eUnbiased]);
0259   m_outputTree->Branch("err_eLOC0_ubs", &m_err_eLOC0[eUnbiased]);
0260   m_outputTree->Branch("err_eLOC1_ubs", &m_err_eLOC1[eUnbiased]);
0261   m_outputTree->Branch("err_ePHI_ubs", &m_err_ePHI[eUnbiased]);
0262   m_outputTree->Branch("err_eTHETA_ubs", &m_err_eTHETA[eUnbiased]);
0263   m_outputTree->Branch("err_eQOP_ubs", &m_err_eQOP[eUnbiased]);
0264   m_outputTree->Branch("err_eT_ubs", &m_err_eT[eUnbiased]);
0265   m_outputTree->Branch("pull_eLOC0_ubs", &m_pull_eLOC0[eUnbiased]);
0266   m_outputTree->Branch("pull_eLOC1_ubs", &m_pull_eLOC1[eUnbiased]);
0267   m_outputTree->Branch("pull_ePHI_ubs", &m_pull_ePHI[eUnbiased]);
0268   m_outputTree->Branch("pull_eTHETA_ubs", &m_pull_eTHETA[eUnbiased]);
0269   m_outputTree->Branch("pull_eQOP_ubs", &m_pull_eQOP[eUnbiased]);
0270   m_outputTree->Branch("pull_eT_ubs", &m_pull_eT[eUnbiased]);
0271   m_outputTree->Branch("g_x_ubs", &m_x[eUnbiased]);
0272   m_outputTree->Branch("g_y_ubs", &m_y[eUnbiased]);
0273   m_outputTree->Branch("g_z_ubs", &m_z[eUnbiased]);
0274   m_outputTree->Branch("px_ubs", &m_px[eUnbiased]);
0275   m_outputTree->Branch("py_ubs", &m_py[eUnbiased]);
0276   m_outputTree->Branch("pz_ubs", &m_pz[eUnbiased]);
0277   m_outputTree->Branch("eta_ubs", &m_eta[eUnbiased]);
0278   m_outputTree->Branch("pT_ubs", &m_pT[eUnbiased]);
0279 }
0280
0281 RootTrackStatesWriter::~RootTrackStatesWriter() {
0282   m_outputFile->Close();
0283 }
0284
0285 ProcessCode RootTrackStatesWriter::finalize() {
0286   m_outputFile->cd();
0287   m_outputTree->Write();
0288   m_outputFile->Close();
0289   return ProcessCode::SUCCESS;
0290 }
0291
0292 RootTrackStatesWriter::StateType RootTrackStatesWriter::getStateType(
0293     ConstTrackStateProxy state) {
0294   if (state.typeFlags().isOutlier()) {
0295     return StateType::eOutlier;
0296   }
0297   if (state.typeFlags().isMeasurement()) {
0298     return StateType::eMeasurement;
0299   }
0300   if (state.typeFlags().isHole()) {
0301     return StateType::eHole;
0302   }
0303   if (state.typeFlags().isMaterial()) {
0304     return StateType::eMaterial;
0305   }
0306   return StateType::eUnknown;
0307 }
0308
0309 ProcessCode RootTrackStatesWriter::writeT(const AlgorithmContext& ctx,
0310                                           const ConstTrackContainer& tracks) {
0311   constexpr float nan = std::numeric_limits<float>::quiet_NaN();
0312
0313   const Acts::GeometryContext& gctx = ctx.geoContext;
0314   // Read additional input collections
0315   const auto& particles = m_inputParticles(ctx);
0316   const auto& trackParticleMatching = m_inputTrackParticleMatching(ctx);
0317   const auto& simHits = m_inputSimHits(ctx);
0318   const auto& hitSimHitsMap = m_inputMeasurementSimHitsMap(ctx);
0319
0320   // Exclusive access to the tree while writing
0321   std::lock_guard<std::mutex> lock(m_writeMutex);
0322
0323   // Get the event number
0324   m_eventNr = ctx.eventNumber;
0325
0326   for (const auto& track : tracks) {
0327     m_trackNr = track.index();
0328
0329     // Collect the track summary info
0330     m_nMeasurements = track.nMeasurements();
0331     m_nStates = track.nTrackStates();
0332
0333     // Get the majority truth particle to this track
0334     int truthQ = 1;
0335     auto match = trackParticleMatching.find(track.index());
0336     if (match != trackParticleMatching.end() &&
0337         match->second.particle.has_value()) {
0338       // Get the barcode of the majority truth particle
0339       auto barcode = match->second.particle.value();
0340       // Find the truth particle via the barcode
0341       auto ip = particles.find(barcode);
0342       if (ip != particles.end()) {
0343         const auto& particle = *ip;
0344         ACTS_VERBOSE("Find the truth particle with barcode " << barcode << "="
0345                                                              << barcode.hash());
0346         // Get the truth particle charge
0347         truthQ = static_cast<int>(particle.charge());
0348       } else {
0349         ACTS_DEBUG("Truth particle with barcode "
0350                    << barcode << "=" << barcode.hash() << " not found!");
0351       }
0352     }
0353
0354     // Get the trackStates on the trajectory
0355     m_nParams = {0, 0, 0, 0};
0356
0357     std::vector<std::uint32_t> particleVertexPrimary;
0358     std::vector<std::uint32_t> particleVertexSecondary;
0359     std::vector<std::uint32_t> particleParticle;
0360     std::vector<std::uint32_t> particleGeneration;
0361     std::vector<std::uint32_t> particleSubParticle;
0362
0363     for (const auto& state : track.trackStatesReversed()) {
0364       const Acts::Surface& surface = state.referenceSurface();
0365
0366       // get the geometry ID
0367       const Acts::GeometryIdentifier geoID = surface.geometryId();
0368       m_volumeID.push_back(geoID.volume());
0369       m_layerID.push_back(geoID.layer());
0370       m_moduleID.push_back(geoID.sensitive());
0371
0372       m_stateType.push_back(Acts::toUnderlying(getStateType(state)));
0373
0374       // get the path length
0375       m_pathLength.push_back(state.pathLength());
0376
0377       // fill the chi2
0378       m_chi2.push_back(state.chi2());
0379
0380       // the truth track parameter at this track state
0381       Acts::BoundVector truthParams;
0382
0383       particleVertexPrimary.clear();
0384       particleVertexSecondary.clear();
0385       particleParticle.clear();
0386       particleGeneration.clear();
0387       particleSubParticle.clear();
0388
0389       if (!state.hasUncalibratedSourceLink()) {
0390         m_t_x.push_back(nan);
0391         m_t_y.push_back(nan);
0392         m_t_z.push_back(nan);
0393         m_t_r.push_back(nan);
0394         m_t_dx.push_back(nan);
0395         m_t_dy.push_back(nan);
0396         m_t_dz.push_back(nan);
0397         m_t_eLOC0.push_back(nan);
0398         m_t_eLOC1.push_back(nan);
0399         m_t_ePHI.push_back(nan);
0400         m_t_eTHETA.push_back(nan);
0401         m_t_eQOP.push_back(nan);
0402         m_t_eT.push_back(nan);
0403
0404         m_lx_hit.push_back(nan);
0405         m_ly_hit.push_back(nan);
0406         m_x_hit.push_back(nan);
0407         m_y_hit.push_back(nan);
0408         m_z_hit.push_back(nan);
0409       } else {
0410         // get the truth hits corresponding to this trackState
0411         // Use average truth in the case of multiple contributing sim hits
0412         if (state.getUncalibratedSourceLink()
0413                 .template getPtr<IndexSourceLink>() == nullptr) {
0414           continue;
0415         }
0416         const auto sl =
0417             state.getUncalibratedSourceLink().template get<IndexSourceLink>();
0418
0419         const auto hitIdx = sl.index();
0420         const auto indices = makeRange(hitSimHitsMap.equal_range(hitIdx));
0421         const auto [truthLocal, truthPos4, truthUnitDir] =
0422             averageSimHits(ctx.geoContext, surface, simHits, indices, logger());
0423
0424         // momentum averaging makes even less sense than averaging position and
0425         // direction. use the first momentum or set q/p to zero
0426         if (!indices.empty()) {
0427           // we assume that the indices are within valid ranges so we do not
0428           // need to check their validity again.
0429           const auto simHitIdx0 = indices.begin()->second;
0430           const auto& simHit0 = *simHits.nth(simHitIdx0);
0431           const double p =
0432               simHit0.momentum4Before().template segment<3>(Acts::eMom0).norm();
0433           truthParams[Acts::eBoundQOverP] = truthQ / p;
0434
0435           // extract particle ids contributed to this track state
0436           for (auto const& [key, simHitIdx] : indices) {
0437             const auto& simHit = *simHits.nth(simHitIdx);
0438             const auto barcode = simHit.particleId();
0439             particleVertexPrimary.push_back(barcode.vertexPrimary());
0440             particleVertexSecondary.push_back(barcode.vertexSecondary());
0441             particleParticle.push_back(barcode.particle());
0442             particleGeneration.push_back(barcode.generation());
0443             particleSubParticle.push_back(barcode.subParticle());
0444           }
0445         }
0446
0447         // fill the truth hit info
0448         m_t_x.push_back(Acts::clampValue<float>(truthPos4[Acts::ePos0]));
0449         m_t_y.push_back(Acts::clampValue<float>(truthPos4[Acts::ePos1]));
0450         m_t_z.push_back(Acts::clampValue<float>(truthPos4[Acts::ePos2]));
0451         m_t_r.push_back(Acts::clampValue<float>(
0452             perp(truthPos4.template segment<3>(Acts::ePos0))));
0453         m_t_dx.push_back(Acts::clampValue<float>(truthUnitDir[Acts::eMom0]));
0454         m_t_dy.push_back(Acts::clampValue<float>(truthUnitDir[Acts::eMom1]));
0455         m_t_dz.push_back(Acts::clampValue<float>(truthUnitDir[Acts::eMom2]));
0456
0457         // get the truth track parameter at this track State
0458         truthParams[Acts::eBoundLoc0] = truthLocal[Acts::ePos0];
0459         truthParams[Acts::eBoundLoc1] = truthLocal[Acts::ePos1];
0460         truthParams[Acts::eBoundPhi] = phi(truthUnitDir);
0461         truthParams[Acts::eBoundTheta] = theta(truthUnitDir);
0462         truthParams[Acts::eBoundTime] = truthPos4[Acts::eTime];
0463
0464         // fill the truth track parameter at this track State
0465         m_t_eLOC0.push_back(
0466             Acts::clampValue<float>(truthParams[Acts::eBoundLoc0]));
0467         m_t_eLOC1.push_back(
0468             Acts::clampValue<float>(truthParams[Acts::eBoundLoc1]));
0469         m_t_ePHI.push_back(
0470             Acts::clampValue<float>(truthParams[Acts::eBoundPhi]));
0471         m_t_eTHETA.push_back(
0472             Acts::clampValue<float>(truthParams[Acts::eBoundTheta]));
0473         m_t_eQOP.push_back(
0474             Acts::clampValue<float>(truthParams[Acts::eBoundQOverP]));
0475         m_t_eT.push_back(
0476             Acts::clampValue<float>(truthParams[Acts::eBoundTime]));
0477
0478         // expand the local measurements into the full bound space
0479         const Acts::BoundVector meas =
0480             state.projectorSubspaceHelper().expandVector(
0481                 state.effectiveCalibrated());
0482         // extract local and global position
0483         const Acts::Vector2 local(meas[Acts::eBoundLoc0],
0484                                   meas[Acts::eBoundLoc1]);
0485         const Acts::Vector3 global =
0486             surface.localToGlobal(ctx.geoContext, local, truthUnitDir);
0487
0488         // fill the measurement info
0489         m_lx_hit.push_back(Acts::clampValue<float>(local[Acts::ePos0]));
0490         m_ly_hit.push_back(Acts::clampValue<float>(local[Acts::ePos1]));
0491         m_x_hit.push_back(Acts::clampValue<float>(global[Acts::ePos0]));
0492         m_y_hit.push_back(Acts::clampValue<float>(global[Acts::ePos1]));
0493         m_z_hit.push_back(Acts::clampValue<float>(global[Acts::ePos2]));
0494       }
0495
0496       // lambda to get the fitted track parameters
0497       auto getTrackParams = [&](unsigned int ipar)
0498           -> std::optional<std::pair<Acts::BoundVector, Acts::BoundMatrix>> {
0499         if (ipar == ePredicted && state.hasPredicted()) {
0500           return std::pair(state.predicted(), state.predictedCovariance());
0501         }
0502         if (ipar == eFiltered && state.hasFiltered()) {
0503           return std::pair(state.filtered(), state.filteredCovariance());
0504         }
0505         if (ipar == eSmoothed && state.hasSmoothed()) {
0506           return std::pair(state.smoothed(), state.smoothedCovariance());
0507         }
0508         if (ipar == eUnbiased && state.hasSmoothed() && state.hasProjector() &&
0509             state.hasCalibrated()) {
0510           return Acts::calculateUnbiasedParametersCovariance(state);
0511         }
0512         return std::nullopt;
0513       };
0514
0515       // fill the fitted track parameters
0516       for (unsigned int ipar = 0; ipar < eSize; ++ipar) {
0517         // get the fitted track parameters
0518         const auto trackParamsOpt = getTrackParams(ipar);
0519         // fill the track parameters status
0520         m_hasParams[ipar].push_back(trackParamsOpt.has_value());
0521
0522         if (!trackParamsOpt.has_value()) {
0523           if (ipar == ePredicted) {
0524             // push default values if no track parameters
0525             m_res_x_hit.push_back(nan);
0526             m_res_y_hit.push_back(nan);
0527             m_err_x_hit.push_back(nan);
0528             m_err_y_hit.push_back(nan);
0529             m_pull_x_hit.push_back(nan);
0530             m_pull_y_hit.push_back(nan);
0531             m_dim_hit.push_back(0);
0532           }
0533
0534           // push default values if no track parameters
0535           m_eLOC0[ipar].push_back(nan);
0536           m_eLOC1[ipar].push_back(nan);
0537           m_ePHI[ipar].push_back(nan);
0538           m_eTHETA[ipar].push_back(nan);
0539           m_eQOP[ipar].push_back(nan);
0540           m_eT[ipar].push_back(nan);
0541           m_res_eLOC0[ipar].push_back(nan);
0542           m_res_eLOC1[ipar].push_back(nan);
0543           m_res_ePHI[ipar].push_back(nan);
0544           m_res_eTHETA[ipar].push_back(nan);
0545           m_res_eQOP[ipar].push_back(nan);
0546           m_res_eT[ipar].push_back(nan);
0547           m_err_eLOC0[ipar].push_back(nan);
0548           m_err_eLOC1[ipar].push_back(nan);
0549           m_err_ePHI[ipar].push_back(nan);
0550           m_err_eTHETA[ipar].push_back(nan);
0551           m_err_eQOP[ipar].push_back(nan);
0552           m_err_eT[ipar].push_back(nan);
0553           m_pull_eLOC0[ipar].push_back(nan);
0554           m_pull_eLOC1[ipar].push_back(nan);
0555           m_pull_ePHI[ipar].push_back(nan);
0556           m_pull_eTHETA[ipar].push_back(nan);
0557           m_pull_eQOP[ipar].push_back(nan);
0558           m_pull_eT[ipar].push_back(nan);
0559           m_x[ipar].push_back(nan);
0560           m_y[ipar].push_back(nan);
0561           m_z[ipar].push_back(nan);
0562           m_px[ipar].push_back(nan);
0563           m_py[ipar].push_back(nan);
0564           m_pz[ipar].push_back(nan);
0565           m_pT[ipar].push_back(nan);
0566           m_eta[ipar].push_back(nan);
0567
0568           continue;
0569         }
0570
0571         ++m_nParams[ipar];
0572         const auto& [parameters, covariance] = *trackParamsOpt;
0573
0574         // track parameters
0575         m_eLOC0[ipar].push_back(
0576             Acts::clampValue<float>(parameters[Acts::eBoundLoc0]));
0577         m_eLOC1[ipar].push_back(
0578             Acts::clampValue<float>(parameters[Acts::eBoundLoc1]));
0579         m_ePHI[ipar].push_back(
0580             Acts::clampValue<float>(parameters[Acts::eBoundPhi]));
0581         m_eTHETA[ipar].push_back(
0582             Acts::clampValue<float>(parameters[Acts::eBoundTheta]));
0583         m_eQOP[ipar].push_back(
0584             Acts::clampValue<float>(parameters[Acts::eBoundQOverP]));
0585         m_eT[ipar].push_back(
0586             Acts::clampValue<float>(parameters[Acts::eBoundTime]));
0587
0588         // track parameters error
0589         Acts::BoundVector errors;
0590         for (Eigen::Index i = 0; i < parameters.size(); ++i) {
0591           const double variance = covariance(i, i);
0592           errors[i] = variance >= 0 ? std::sqrt(variance) : nan;
0593         }
0594         m_err_eLOC0[ipar].push_back(
0595             Acts::clampValue<float>(errors[Acts::eBoundLoc0]));
0596         m_err_eLOC1[ipar].push_back(
0597             Acts::clampValue<float>(errors[Acts::eBoundLoc1]));
0598         m_err_ePHI[ipar].push_back(
0599             Acts::clampValue<float>(errors[Acts::eBoundPhi]));
0600         m_err_eTHETA[ipar].push_back(
0601             Acts::clampValue<float>(errors[Acts::eBoundTheta]));
0602         m_err_eQOP[ipar].push_back(
0603             Acts::clampValue<float>(errors[Acts::eBoundQOverP]));
0604         m_err_eT[ipar].push_back(
0605             Acts::clampValue<float>(errors[Acts::eBoundTime]));
0606
0607         // further track parameter info
0608         const Acts::FreeVector freeParams =
0609             Acts::transformBoundToFreeParameters(surface, gctx, parameters);
0610         m_x[ipar].push_back(
0611             Acts::clampValue<float>(freeParams[Acts::eFreePos0]));
0612         m_y[ipar].push_back(
0613             Acts::clampValue<float>(freeParams[Acts::eFreePos1]));
0614         m_z[ipar].push_back(
0615             Acts::clampValue<float>(freeParams[Acts::eFreePos2]));
0616         // single charge assumption
0617         const double p = std::abs(1 / freeParams[Acts::eFreeQOverP]);
0618         m_px[ipar].push_back(
0619             Acts::clampValue<float>(p * freeParams[Acts::eFreeDir0]));
0620         m_py[ipar].push_back(
0621             Acts::clampValue<float>(p * freeParams[Acts::eFreeDir1]));
0622         m_pz[ipar].push_back(
0623             Acts::clampValue<float>(p * freeParams[Acts::eFreeDir2]));
0624         m_pT[ipar].push_back(Acts::clampValue<float>(
0625             p * std::hypot(freeParams[Acts::eFreeDir0],
0626                            freeParams[Acts::eFreeDir1])));
0627         m_eta[ipar].push_back(Acts::clampValue<float>(
0628             Acts::VectorHelpers::eta(freeParams.segment<3>(Acts::eFreeDir0))));
0629
0630         if (!state.hasUncalibratedSourceLink()) {
0631           continue;
0632         }
0633
0634         // track parameters residual
0635         Acts::BoundVector residuals = parameters - truthParams;
0636         residuals[Acts::eBoundPhi] = Acts::detail::difference_periodic(
0637             parameters[Acts::eBoundPhi], truthParams[Acts::eBoundPhi],
0638             2 * std::numbers::pi);
0639         m_res_eLOC0[ipar].push_back(
0640             Acts::clampValue<float>(residuals[Acts::eBoundLoc0]));
0641         m_res_eLOC1[ipar].push_back(
0642             Acts::clampValue<float>(residuals[Acts::eBoundLoc1]));
0643         m_res_ePHI[ipar].push_back(
0644             Acts::clampValue<float>(residuals[Acts::eBoundPhi]));
0645         m_res_eTHETA[ipar].push_back(
0646             Acts::clampValue<float>(residuals[Acts::eBoundTheta]));
0647         m_res_eQOP[ipar].push_back(
0648             Acts::clampValue<float>(residuals[Acts::eBoundQOverP]));
0649         m_res_eT[ipar].push_back(
0650             Acts::clampValue<float>(residuals[Acts::eBoundTime]));
0651
0652         // track parameters pull
0653         Acts::BoundVector pulls = Acts::BoundVector::Constant(nan);
0654         for (Eigen::Index i = 0; i < parameters.size(); ++i) {
0655           pulls[i] = (!std::isnan(errors[i]) && errors[i] > 0)
0656                          ? residuals[i] / errors[i]
0657                          : nan;
0658         }
0659         m_pull_eLOC0[ipar].push_back(
0660             Acts::clampValue<float>(pulls[Acts::eBoundLoc0]));
0661         m_pull_eLOC1[ipar].push_back(
0662             Acts::clampValue<float>(pulls[Acts::eBoundLoc1]));
0663         m_pull_ePHI[ipar].push_back(
0664             Acts::clampValue<float>(pulls[Acts::eBoundPhi]));
0665         m_pull_eTHETA[ipar].push_back(
0666             Acts::clampValue<float>(pulls[Acts::eBoundTheta]));
0667         m_pull_eQOP[ipar].push_back(
0668             Acts::clampValue<float>(pulls[Acts::eBoundQOverP]));
0669         m_pull_eT[ipar].push_back(
0670             Acts::clampValue<float>(pulls[Acts::eBoundTime]));
0671
0672         if (ipar == ePredicted) {
0673           // local hit residual info
0674           const Acts::DynamicMatrix H =
0675               state.projectorSubspaceHelper().fullProjector().topLeftCorner(
0676                   state.calibratedSize(), Acts::eBoundSize);
0677           const Acts::DynamicMatrix V = state.effectiveCalibratedCovariance();
0678           const Acts::DynamicMatrix resCov = V + H * covariance * H.transpose();
0679           const Acts::DynamicVector res =
0680               state.effectiveCalibrated() - H * parameters;
0681
0682           const double resX = res[Acts::eBoundLoc0];
0683           const double errX =
0684               V(Acts::eBoundLoc0, Acts::eBoundLoc0) >= 0
0685                   ? std::sqrt(V(Acts::eBoundLoc0, Acts::eBoundLoc0))
0686                   : nan;
0687           const double pullX =
0688               resCov(Acts::eBoundLoc0, Acts::eBoundLoc0) > 0
0689                   ? resX / std::sqrt(resCov(Acts::eBoundLoc0, Acts::eBoundLoc0))
0690                   : nan;
0691
0692           m_res_x_hit.push_back(Acts::clampValue<float>(resX));
0693           m_err_x_hit.push_back(Acts::clampValue<float>(errX));
0694           m_pull_x_hit.push_back(Acts::clampValue<float>(pullX));
0695
0696           if (state.calibratedSize() >= 2) {
0697             const double resY = res[Acts::eBoundLoc1];
0698             const double errY =
0699                 V(Acts::eBoundLoc1, Acts::eBoundLoc1) >= 0
0700                     ? std::sqrt(V(Acts::eBoundLoc1, Acts::eBoundLoc1))
0701                     : nan;
0702             const double pullY =
0703                 resCov(Acts::eBoundLoc1, Acts::eBoundLoc1) > 0
0704                     ? resY /
0705                           std::sqrt(resCov(Acts::eBoundLoc1, Acts::eBoundLoc1))
0706                     : nan;
0707
0708             m_res_y_hit.push_back(Acts::clampValue<float>(resY));
0709             m_err_y_hit.push_back(Acts::clampValue<float>(errY));
0710             m_pull_y_hit.push_back(Acts::clampValue<float>(pullY));
0711           } else {
0712             m_res_y_hit.push_back(nan);
0713             m_err_y_hit.push_back(nan);
0714             m_pull_y_hit.push_back(nan);
0715           }
0716
0717           m_dim_hit.push_back(state.calibratedSize());
0718         }
0719       }
0720       m_particleVertexPrimary.push_back(std::move(particleVertexPrimary));
0721       m_particleVertexSecondary.push_back(std::move(particleVertexSecondary));
0722       m_particleParticle.push_back(std::move(particleParticle));
0723       m_particleGeneration.push_back(std::move(particleGeneration));
0724       m_particleSubParticle.push_back(std::move(particleSubParticle));
0725     }
0726
0727     // fill the variables for one track to tree
0728     m_outputTree->Fill();
0729
0730     // now reset
0731     m_volumeID.clear();
0732     m_layerID.clear();
0733     m_moduleID.clear();
0734
0735     m_stateType.clear();
0736
0737     m_chi2.clear();
0738
0739     m_pathLength.clear();
0740
0741     m_t_x.clear();
0742     m_t_y.clear();
0743     m_t_z.clear();
0744     m_t_r.clear();
0745     m_t_dx.clear();
0746     m_t_dy.clear();
0747     m_t_dz.clear();
0748     m_t_eLOC0.clear();
0749     m_t_eLOC1.clear();
0750     m_t_ePHI.clear();
0751     m_t_eTHETA.clear();
0752     m_t_eQOP.clear();
0753     m_t_eT.clear();
0754     m_particleVertexPrimary.clear();
0755     m_particleVertexSecondary.clear();
0756     m_particleParticle.clear();
0757     m_particleGeneration.clear();
0758     m_particleSubParticle.clear();
0759
0760     m_dim_hit.clear();
0761     m_lx_hit.clear();
0762     m_ly_hit.clear();
0763     m_x_hit.clear();
0764     m_y_hit.clear();
0765     m_z_hit.clear();
0766     m_res_x_hit.clear();
0767     m_res_y_hit.clear();
0768     m_err_x_hit.clear();
0769     m_err_y_hit.clear();
0770     m_pull_x_hit.clear();
0771     m_pull_y_hit.clear();
0772
0773     for (unsigned int ipar = 0; ipar < eSize; ++ipar) {
0774       m_hasParams[ipar].clear();
0775       m_eLOC0[ipar].clear();
0776       m_eLOC1[ipar].clear();
0777       m_ePHI[ipar].clear();
0778       m_eTHETA[ipar].clear();
0779       m_eQOP[ipar].clear();
0780       m_eT[ipar].clear();
0781       m_res_eLOC0[ipar].clear();
0782       m_res_eLOC1[ipar].clear();
0783       m_res_ePHI[ipar].clear();
0784       m_res_eTHETA[ipar].clear();
0785       m_res_eQOP[ipar].clear();
0786       m_res_eT[ipar].clear();
0787       m_err_eLOC0[ipar].clear();
0788       m_err_eLOC1[ipar].clear();
0789       m_err_ePHI[ipar].clear();
0790       m_err_eTHETA[ipar].clear();
0791       m_err_eQOP[ipar].clear();
0792       m_err_eT[ipar].clear();
0793       m_pull_eLOC0[ipar].clear();
0794       m_pull_eLOC1[ipar].clear();
0795       m_pull_ePHI[ipar].clear();
0796       m_pull_eTHETA[ipar].clear();
0797       m_pull_eQOP[ipar].clear();
0798       m_pull_eT[ipar].clear();
0799       m_x[ipar].clear();
0800       m_y[ipar].clear();
0801       m_z[ipar].clear();
0802       m_px[ipar].clear();
0803       m_py[ipar].clear();
0804       m_pz[ipar].clear();
0805       m_eta[ipar].clear();
0806       m_pT[ipar].clear();
0807     }
0808
0809     m_chi2.clear();
0810   }
0811
0812   return ProcessCode::SUCCESS;
0813 }
0814
0815 }  // namespace ActsExamples