Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-12-16 09:24:01 

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

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!