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