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