Root/src/RootTrackStatesWriter.cpp

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