Root/src/RootTrackSummaryWriter.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/RootTrackSummaryWriter.hpp"
0010
0011 #include "Acts/Definitions/TrackParametrization.hpp"
0012 #include "Acts/EventData/AnyTrackProxy.hpp"
0013 #include "Acts/EventData/VectorMultiTrajectory.hpp"
0014 #include "Acts/Surfaces/Surface.hpp"
0015 #include "Acts/TrackFitting/GsfOptions.hpp"
0016 #include "Acts/Utilities/Intersection.hpp"
0017 #include "Acts/Utilities/Result.hpp"
0018 #include "Acts/Utilities/detail/periodic.hpp"
0019 #include "ActsExamples/EventData/TruthMatching.hpp"
0020 #include "ActsExamples/Framework/AlgorithmContext.hpp"
0021 #include "ActsExamples/Framework/WriterT.hpp"
0022 #include "ActsExamples/Validation/TrackClassification.hpp"
0023 #include "ActsFatras/EventData/Barcode.hpp"
0024
0025 #include <array>
0026 #include <cmath>
0027 #include <cstddef>
0028 #include <cstdint>
0029 #include <ios>
0030 #include <limits>
0031 #include <numbers>
0032 #include <optional>
0033 #include <ostream>
0034 #include <stdexcept>
0035
0036 #include <TFile.h>
0037 #include <TTree.h>
0038
0039 using Acts::VectorHelpers::eta;
0040 using Acts::VectorHelpers::perp;
0041 using Acts::VectorHelpers::phi;
0042 using Acts::VectorHelpers::theta;
0043
0044 namespace ActsExamples {
0045
0046 RootTrackSummaryWriter::RootTrackSummaryWriter(
0047     const RootTrackSummaryWriter::Config& config, Acts::Logging::Level level)
0048     : WriterT(config.inputTracks, "RootTrackSummaryWriter", level),
0049       m_cfg(config) {
0050   // tracks collection name is already checked by base ctor
0051   if (m_cfg.filePath.empty()) {
0052     throw std::invalid_argument("Missing output filename");
0053   }
0054   if (m_cfg.treeName.empty()) {
0055     throw std::invalid_argument("Missing tree name");
0056   }
0057
0058   m_inputParticles.maybeInitialize(m_cfg.inputParticles);
0059   m_inputTrackParticleMatching.maybeInitialize(
0060       m_cfg.inputTrackParticleMatching);
0061   if (m_cfg.writeJets) {
0062     m_inputJets.maybeInitialize(m_cfg.inputJets);
0063   }
0064
0065   // Setup ROOT I/O
0066   auto path = m_cfg.filePath;
0067   m_outputFile = TFile::Open(path.c_str(), m_cfg.fileMode.c_str());
0068   if (m_outputFile == nullptr) {
0069     throw std::ios_base::failure("Could not open '" + path + "'");
0070   }
0071   m_outputFile->cd();
0072   m_outputTree = new TTree(m_cfg.treeName.c_str(), m_cfg.treeName.c_str());
0073   if (m_outputTree == nullptr) {
0074     throw std::bad_alloc();
0075   }
0076
0077   // I/O parameters
0078   m_outputTree->Branch("event_nr", &m_eventNr);
0079   m_outputTree->Branch("track_nr", &m_trackNr);
0080
0081   m_outputTree->Branch("nStates", &m_nStates);
0082   m_outputTree->Branch("nMeasurements", &m_nMeasurements);
0083   m_outputTree->Branch("nOutliers", &m_nOutliers);
0084   m_outputTree->Branch("nHoles", &m_nHoles);
0085   m_outputTree->Branch("nSharedHits", &m_nSharedHits);
0086   m_outputTree->Branch("chi2Sum", &m_chi2Sum);
0087   m_outputTree->Branch("NDF", &m_NDF);
0088   m_outputTree->Branch("measurementChi2", &m_measurementChi2);
0089   m_outputTree->Branch("outlierChi2", &m_outlierChi2);
0090   m_outputTree->Branch("measurementVolume", &m_measurementVolume);
0091   m_outputTree->Branch("measurementLayer", &m_measurementLayer);
0092   m_outputTree->Branch("outlierVolume", &m_outlierVolume);
0093   m_outputTree->Branch("outlierLayer", &m_outlierLayer);
0094
0095   m_outputTree->Branch("nMajorityHits", &m_nMajorityHits);
0096   m_outputTree->Branch("majorityParticleId_vertex_primary",
0097                        &m_majorityParticleVertexPrimary);
0098   m_outputTree->Branch("majorityParticleId_vertex_secondary",
0099                        &m_majorityParticleVertexSecondary);
0100   m_outputTree->Branch("majorityParticleId_particle",
0101                        &m_majorityParticleParticle);
0102   m_outputTree->Branch("majorityParticleId_generation",
0103                        &m_majorityParticleGeneration);
0104   m_outputTree->Branch("majorityParticleId_sub_particle",
0105                        &m_majorityParticleSubParticle);
0106   m_outputTree->Branch("trackClassification", &m_trackClassification);
0107   m_outputTree->Branch("t_charge", &m_t_charge);
0108   m_outputTree->Branch("t_time", &m_t_time);
0109   m_outputTree->Branch("t_vx", &m_t_vx);
0110   m_outputTree->Branch("t_vy", &m_t_vy);
0111   m_outputTree->Branch("t_vz", &m_t_vz);
0112   m_outputTree->Branch("t_px", &m_t_px);
0113   m_outputTree->Branch("t_py", &m_t_py);
0114   m_outputTree->Branch("t_pz", &m_t_pz);
0115   m_outputTree->Branch("t_theta", &m_t_theta);
0116   m_outputTree->Branch("t_phi", &m_t_phi);
0117   m_outputTree->Branch("t_eta", &m_t_eta);
0118   m_outputTree->Branch("t_p", &m_t_p);
0119   m_outputTree->Branch("t_pT", &m_t_pT);
0120   m_outputTree->Branch("t_d0", &m_t_d0);
0121   m_outputTree->Branch("t_z0", &m_t_z0);
0122   m_outputTree->Branch("t_prodR", &m_t_prodR);
0123
0124   m_outputTree->Branch("hasFittedParams", &m_hasFittedParams);
0125   m_outputTree->Branch("eLOC0_fit", &m_eLOC0_fit);
0126   m_outputTree->Branch("eLOC1_fit", &m_eLOC1_fit);
0127   m_outputTree->Branch("ePHI_fit", &m_ePHI_fit);
0128   m_outputTree->Branch("eTHETA_fit", &m_eTHETA_fit);
0129   m_outputTree->Branch("eQOP_fit", &m_eQOP_fit);
0130   m_outputTree->Branch("eT_fit", &m_eT_fit);
0131   m_outputTree->Branch("err_eLOC0_fit", &m_err_eLOC0_fit);
0132   m_outputTree->Branch("err_eLOC1_fit", &m_err_eLOC1_fit);
0133   m_outputTree->Branch("err_ePHI_fit", &m_err_ePHI_fit);
0134   m_outputTree->Branch("err_eTHETA_fit", &m_err_eTHETA_fit);
0135   m_outputTree->Branch("err_eQOP_fit", &m_err_eQOP_fit);
0136   m_outputTree->Branch("err_eT_fit", &m_err_eT_fit);
0137   m_outputTree->Branch("res_eLOC0_fit", &m_res_eLOC0_fit);
0138   m_outputTree->Branch("res_eLOC1_fit", &m_res_eLOC1_fit);
0139   m_outputTree->Branch("res_ePHI_fit", &m_res_ePHI_fit);
0140   m_outputTree->Branch("res_eTHETA_fit", &m_res_eTHETA_fit);
0141   m_outputTree->Branch("res_eQOP_fit", &m_res_eQOP_fit);
0142   m_outputTree->Branch("res_eT_fit", &m_res_eT_fit);
0143   m_outputTree->Branch("pull_eLOC0_fit", &m_pull_eLOC0_fit);
0144   m_outputTree->Branch("pull_eLOC1_fit", &m_pull_eLOC1_fit);
0145   m_outputTree->Branch("pull_ePHI_fit", &m_pull_ePHI_fit);
0146   m_outputTree->Branch("pull_eTHETA_fit", &m_pull_eTHETA_fit);
0147   m_outputTree->Branch("pull_eQOP_fit", &m_pull_eQOP_fit);
0148   m_outputTree->Branch("pull_eT_fit", &m_pull_eT_fit);
0149
0150   if (m_cfg.writeGsfSpecific) {
0151     m_outputTree->Branch("max_material_fwd", &m_gsf_max_material_fwd);
0152     m_outputTree->Branch("sum_material_fwd", &m_gsf_sum_material_fwd);
0153   }
0154
0155   if (m_cfg.writeCovMat) {
0156     // create one branch for every entry of covariance matrix
0157     // one block for every row of the matrix, every entry gets own branch
0158     m_outputTree->Branch("cov_eLOC0_eLOC0", &m_cov_eLOC0_eLOC0);
0159     m_outputTree->Branch("cov_eLOC0_eLOC1", &m_cov_eLOC0_eLOC1);
0160     m_outputTree->Branch("cov_eLOC0_ePHI", &m_cov_eLOC0_ePHI);
0161     m_outputTree->Branch("cov_eLOC0_eTHETA", &m_cov_eLOC0_eTHETA);
0162     m_outputTree->Branch("cov_eLOC0_eQOP", &m_cov_eLOC0_eQOP);
0163     m_outputTree->Branch("cov_eLOC0_eT", &m_cov_eLOC0_eT);
0164
0165     m_outputTree->Branch("cov_eLOC1_eLOC0", &m_cov_eLOC1_eLOC0);
0166     m_outputTree->Branch("cov_eLOC1_eLOC1", &m_cov_eLOC1_eLOC1);
0167     m_outputTree->Branch("cov_eLOC1_ePHI", &m_cov_eLOC1_ePHI);
0168     m_outputTree->Branch("cov_eLOC1_eTHETA", &m_cov_eLOC1_eTHETA);
0169     m_outputTree->Branch("cov_eLOC1_eQOP", &m_cov_eLOC1_eQOP);
0170     m_outputTree->Branch("cov_eLOC1_eT", &m_cov_eLOC1_eT);
0171
0172     m_outputTree->Branch("cov_ePHI_eLOC0", &m_cov_ePHI_eLOC0);
0173     m_outputTree->Branch("cov_ePHI_eLOC1", &m_cov_ePHI_eLOC1);
0174     m_outputTree->Branch("cov_ePHI_ePHI", &m_cov_ePHI_ePHI);
0175     m_outputTree->Branch("cov_ePHI_eTHETA", &m_cov_ePHI_eTHETA);
0176     m_outputTree->Branch("cov_ePHI_eQOP", &m_cov_ePHI_eQOP);
0177     m_outputTree->Branch("cov_ePHI_eT", &m_cov_ePHI_eT);
0178
0179     m_outputTree->Branch("cov_eTHETA_eLOC0", &m_cov_eTHETA_eLOC0);
0180     m_outputTree->Branch("cov_eTHETA_eLOC1", &m_cov_eTHETA_eLOC1);
0181     m_outputTree->Branch("cov_eTHETA_ePHI", &m_cov_eTHETA_ePHI);
0182     m_outputTree->Branch("cov_eTHETA_eTHETA", &m_cov_eTHETA_eTHETA);
0183     m_outputTree->Branch("cov_eTHETA_eQOP", &m_cov_eTHETA_eQOP);
0184     m_outputTree->Branch("cov_eTHETA_eT", &m_cov_eTHETA_eT);
0185
0186     m_outputTree->Branch("cov_eQOP_eLOC0", &m_cov_eQOP_eLOC0);
0187     m_outputTree->Branch("cov_eQOP_eLOC1", &m_cov_eQOP_eLOC1);
0188     m_outputTree->Branch("cov_eQOP_ePHI", &m_cov_eQOP_ePHI);
0189     m_outputTree->Branch("cov_eQOP_eTHETA", &m_cov_eQOP_eTHETA);
0190     m_outputTree->Branch("cov_eQOP_eQOP", &m_cov_eQOP_eQOP);
0191     m_outputTree->Branch("cov_eQOP_eT", &m_cov_eQOP_eT);
0192
0193     m_outputTree->Branch("cov_eT_eLOC0", &m_cov_eT_eLOC0);
0194     m_outputTree->Branch("cov_eT_eLOC1", &m_cov_eT_eLOC1);
0195     m_outputTree->Branch("cov_eT_ePHI", &m_cov_eT_ePHI);
0196     m_outputTree->Branch("cov_eT_eTHETA", &m_cov_eT_eTHETA);
0197     m_outputTree->Branch("cov_eT_eQOP", &m_cov_eT_eQOP);
0198     m_outputTree->Branch("cov_eT_eT", &m_cov_eT_eT);
0199   }
0200
0201   if (m_cfg.writeGx2fSpecific) {
0202     m_outputTree->Branch("nUpdatesGx2f", &m_nUpdatesGx2f);
0203   }
0204
0205   if (m_cfg.writeJets) {
0206     m_outputTree->Branch("nJets", &m_nJets);
0207     m_outputTree->Branch("jet_pt", &m_jet_pt);
0208     m_outputTree->Branch("jet_eta", &m_jet_eta);
0209     m_outputTree->Branch("jet_phi", &m_jet_phi);
0210     m_outputTree->Branch("jet_label", &m_jet_label);
0211     m_outputTree->Branch("ntracks_per_jets", &m_ntracks_per_jets);
0212   }
0213 }
0214
0215 RootTrackSummaryWriter::~RootTrackSummaryWriter() {
0216   m_outputFile->Close();
0217 }
0218
0219 ProcessCode RootTrackSummaryWriter::finalize() {
0220   m_outputFile->cd();
0221   m_outputTree->Write();
0222   m_outputFile->Close();
0223
0224   if (m_cfg.writeCovMat) {
0225     ACTS_INFO("Wrote full covariance matrix to tree");
0226   }
0227   ACTS_INFO("Wrote parameters of tracks to tree '" << m_cfg.treeName << "' in '"
0228                                                    << m_cfg.filePath << "'");
0229
0230   return ProcessCode::SUCCESS;
0231 }
0232
0233 ProcessCode RootTrackSummaryWriter::writeT(const AlgorithmContext& ctx,
0234                                            const ConstTrackContainer& tracks) {
0235   // In case we do not have truth info, we bind to a empty collection
0236   const static SimParticleContainer emptyParticles;
0237   const static TrackParticleMatching emptyTrackParticleMatching;
0238
0239   const auto& particles =
0240       m_inputParticles.isInitialized() ? m_inputParticles(ctx) : emptyParticles;
0241   const auto& trackParticleMatching =
0242       m_inputTrackParticleMatching.isInitialized()
0243           ? m_inputTrackParticleMatching(ctx)
0244           : emptyTrackParticleMatching;
0245
0246   // For each particle within a track, how many hits did it contribute
0247   std::vector<ParticleHitCount> particleHitCounts;
0248
0249   // Exclusive access to the tree while writing
0250   std::lock_guard<std::mutex> lock(m_writeMutex);
0251
0252   // Get the event number
0253   m_eventNr = ctx.eventNumber;
0254
0255   std::vector<ActsExamples::TruthJet> jets;
0256   std::unordered_map<std::size_t, std::vector<std::int32_t>>
0257       jetToTrackIndicesMap;
0258
0259   // Fill the jet vector if requested
0260   if (m_cfg.writeJets) {
0261     auto& inputJets = m_inputJets(ctx);
0262     jets = inputJets;
0263
0264     // Loop over jets and fill jet kinematic variables
0265     for (std::size_t ijet = 0; ijet < jets.size(); ++ijet) {
0266       m_nJets.push_back(jets.size());
0267       Acts::Vector4 jet_4mom = jets[ijet].fourMomentum();
0268       Acts::Vector3 jet_3mom{jet_4mom[0], jet_4mom[1], jet_4mom[2]};
0269
0270       float jet_theta = theta(jet_3mom);
0271
0272       m_jet_pt.push_back(perp(jet_4mom));
0273       m_jet_eta.push_back(std::atanh(std::cos(jet_theta)));
0274       m_jet_phi.push_back(phi(jet_4mom));
0275       m_jet_label.push_back(static_cast<int>(jets[ijet].jetLabel()));
0276       m_ntracks_per_jets.push_back(jets[ijet].associatedTracks().size());
0277     }
0278   }
0279
0280   for (const auto& track : tracks) {
0281     m_trackNr.push_back(track.index());
0282
0283     // Collect the trajectory summary info
0284     m_nStates.push_back(track.nTrackStates());
0285     m_nMeasurements.push_back(track.nMeasurements());
0286     m_nOutliers.push_back(track.nOutliers());
0287     m_nHoles.push_back(track.nHoles());
0288     m_nSharedHits.push_back(track.nSharedHits());
0289     m_chi2Sum.push_back(track.chi2());
0290     m_NDF.push_back(track.nDoF());
0291
0292     {
0293       std::vector<double> measurementChi2;
0294       std::vector<std::uint32_t> measurementVolume;
0295       std::vector<std::uint32_t> measurementLayer;
0296       std::vector<double> outlierChi2;
0297       std::vector<std::uint32_t> outlierVolume;
0298       std::vector<std::uint32_t> outlierLayer;
0299       for (const auto& state : track.trackStatesReversed()) {
0300         const auto& geoID = state.referenceSurface().geometryId();
0301         const auto& volume = geoID.volume();
0302         const auto& layer = geoID.layer();
0303         if (state.typeFlags().isOutlier()) {
0304           outlierChi2.push_back(state.chi2());
0305           outlierVolume.push_back(volume);
0306           outlierLayer.push_back(layer);
0307         } else if (state.typeFlags().isMeasurement()) {
0308           measurementChi2.push_back(state.chi2());
0309           measurementVolume.push_back(volume);
0310           measurementLayer.push_back(layer);
0311         }
0312       }
0313       m_measurementChi2.push_back(std::move(measurementChi2));
0314       m_measurementVolume.push_back(std::move(measurementVolume));
0315       m_measurementLayer.push_back(std::move(measurementLayer));
0316       m_outlierChi2.push_back(std::move(outlierChi2));
0317       m_outlierVolume.push_back(std::move(outlierVolume));
0318       m_outlierLayer.push_back(std::move(outlierLayer));
0319     }
0320
0321     // Initialize the truth particle info
0322     SimBarcode majorityParticleId{};
0323     TrackMatchClassification trackClassification =
0324         TrackMatchClassification::Unknown;
0325     unsigned int nMajorityHits = std::numeric_limits<unsigned int>::max();
0326     int t_charge = std::numeric_limits<int>::max();
0327     float t_time = NaNfloat;
0328     float t_vx = NaNfloat;
0329     float t_vy = NaNfloat;
0330     float t_vz = NaNfloat;
0331     float t_px = NaNfloat;
0332     float t_py = NaNfloat;
0333     float t_pz = NaNfloat;
0334     float t_theta = NaNfloat;
0335     float t_phi = NaNfloat;
0336     float t_eta = NaNfloat;
0337     float t_p = NaNfloat;
0338     float t_pT = NaNfloat;
0339     float t_d0 = NaNfloat;
0340     float t_z0 = NaNfloat;
0341     float t_qop = NaNfloat;
0342     float t_prodR = NaNfloat;
0343
0344     // Get the perigee surface
0345     const Acts::Surface* pSurface =
0346         track.hasReferenceSurface() ? &track.referenceSurface() : nullptr;
0347
0348     // Get the majority truth particle to this track
0349     auto match = trackParticleMatching.find(track.index());
0350     if (match != trackParticleMatching.end()) {
0351       trackClassification = match->second.classification;
0352     }
0353     bool foundMajorityParticle = false;
0354     // Get the truth particle info
0355     if (match != trackParticleMatching.end() &&
0356         match->second.particle.has_value()) {
0357       // Get the barcode of the majority truth particle
0358       majorityParticleId = match->second.particle.value();
0359       nMajorityHits = match->second.contributingParticles.front().hitCount;
0360
0361       // Find the truth particle via the barcode
0362       auto ip = particles.find(majorityParticleId);
0363       if (ip != particles.end()) {
0364         foundMajorityParticle = true;
0365
0366         const auto& particle = *ip;
0367         ACTS_VERBOSE("Find the truth particle with barcode "
0368                      << majorityParticleId << "=" << majorityParticleId.hash());
0369         // Get the truth particle info at vertex
0370         t_p = particle.absoluteMomentum();
0371         t_charge = static_cast<int>(particle.charge());
0372         t_time = particle.time();
0373         t_vx = particle.position().x();
0374         t_vy = particle.position().y();
0375         t_vz = particle.position().z();
0376         t_px = t_p * particle.direction().x();
0377         t_py = t_p * particle.direction().y();
0378         t_pz = t_p * particle.direction().z();
0379         t_theta = theta(particle.direction());
0380         t_phi = phi(particle.direction());
0381         t_eta = eta(particle.direction());
0382         t_pT = t_p * perp(particle.direction());
0383         t_qop = particle.qOverP();
0384         t_prodR = std::sqrt(t_vx * t_vx + t_vy * t_vy);
0385
0386         if (pSurface != nullptr) {
0387           // The perigee surface is not aligned, so either context works
0388           Acts::Intersection3D intersection =
0389               pSurface
0390                   ->intersect(ctx.recoGeoContext, particle.position(),
0391                               particle.direction(),
0392                               Acts::BoundaryTolerance::Infinite())
0393                   .closest();
0394           auto position = intersection.position();
0395
0396           // get the truth perigee parameter
0397           auto lpResult = pSurface->globalToLocal(ctx.recoGeoContext, position,
0398                                                   particle.direction());
0399           if (lpResult.ok()) {
0400             t_d0 = lpResult.value()[Acts::BoundIndices::eBoundLoc0];
0401             t_z0 = lpResult.value()[Acts::BoundIndices::eBoundLoc1];
0402           } else {
0403             ACTS_ERROR("Global to local transformation did not succeed.");
0404           }
0405         }
0406       } else {
0407         ACTS_DEBUG("Truth particle with barcode "
0408                    << majorityParticleId << "=" << majorityParticleId.hash()
0409                    << " not found in the input collection!");
0410       }
0411     }
0412     if (!foundMajorityParticle) {
0413       ACTS_DEBUG("Truth particle for track " << track.tipIndex()
0414                                              << " not found!");
0415     }
0416
0417     // Push the corresponding truth particle info for the track.
0418     // Always push back even if majority particle not found
0419     m_majorityParticleVertexPrimary.push_back(
0420         majorityParticleId.vertexPrimary());
0421     m_majorityParticleVertexSecondary.push_back(
0422         majorityParticleId.vertexSecondary());
0423     m_majorityParticleParticle.push_back(majorityParticleId.particle());
0424     m_majorityParticleGeneration.push_back(majorityParticleId.generation());
0425     m_majorityParticleSubParticle.push_back(majorityParticleId.subParticle());
0426     m_trackClassification.push_back(static_cast<int>(trackClassification));
0427     m_nMajorityHits.push_back(nMajorityHits);
0428     m_t_charge.push_back(t_charge);
0429     m_t_time.push_back(t_time);
0430     m_t_vx.push_back(t_vx);
0431     m_t_vy.push_back(t_vy);
0432     m_t_vz.push_back(t_vz);
0433     m_t_px.push_back(t_px);
0434     m_t_py.push_back(t_py);
0435     m_t_pz.push_back(t_pz);
0436     m_t_theta.push_back(t_theta);
0437     m_t_phi.push_back(t_phi);
0438     m_t_eta.push_back(t_eta);
0439     m_t_p.push_back(t_p);
0440     m_t_pT.push_back(t_pT);
0441     m_t_d0.push_back(t_d0);
0442     m_t_z0.push_back(t_z0);
0443     m_t_prodR.push_back(t_prodR);
0444
0445     // Initialize the fitted track parameters info
0446     std::array<float, Acts::eBoundSize> param = {NaNfloat, NaNfloat, NaNfloat,
0447                                                  NaNfloat, NaNfloat, NaNfloat};
0448     std::array<float, Acts::eBoundSize> error = {NaNfloat, NaNfloat, NaNfloat,
0449                                                  NaNfloat, NaNfloat, NaNfloat};
0450
0451     // get entries of covariance matrix. If no entry, return NaN
0452     auto getCov = [&](auto i, auto j) { return track.covariance()(i, j); };
0453
0454     bool hasFittedParams = track.hasReferenceSurface();
0455     if (hasFittedParams) {
0456       const auto& parameter = track.parameters();
0457       for (unsigned int i = 0; i < Acts::eBoundSize; ++i) {
0458         param[i] = parameter[i];
0459       }
0460
0461       for (unsigned int i = 0; i < Acts::eBoundSize; ++i) {
0462         double variance = getCov(i, i);
0463         error[i] = variance >= 0 ? std::sqrt(variance) : NaNfloat;
0464       }
0465     }
0466
0467     std::array<float, Acts::eBoundSize> res = {NaNfloat, NaNfloat, NaNfloat,
0468                                                NaNfloat, NaNfloat, NaNfloat};
0469     std::array<float, Acts::eBoundSize> pull = {NaNfloat, NaNfloat, NaNfloat,
0470                                                 NaNfloat, NaNfloat, NaNfloat};
0471     if (foundMajorityParticle && hasFittedParams) {
0472       res = {param[Acts::eBoundLoc0] - t_d0,
0473              param[Acts::eBoundLoc1] - t_z0,
0474              Acts::detail::difference_periodic(
0475                  param[Acts::eBoundPhi], t_phi,
0476                  static_cast<float>(2 * std::numbers::pi)),
0477              param[Acts::eBoundTheta] - t_theta,
0478              param[Acts::eBoundQOverP] - t_qop,
0479              param[Acts::eBoundTime] - t_time};
0480
0481       for (unsigned int i = 0; i < Acts::eBoundSize; ++i) {
0482         pull[i] = res[i] / error[i];
0483       }
0484     }
0485
0486     // Push the fitted track parameters.
0487     // Always push back even if no fitted track parameters
0488     m_eLOC0_fit.push_back(param[Acts::eBoundLoc0]);
0489     m_eLOC1_fit.push_back(param[Acts::eBoundLoc1]);
0490     m_ePHI_fit.push_back(param[Acts::eBoundPhi]);
0491     m_eTHETA_fit.push_back(param[Acts::eBoundTheta]);
0492     m_eQOP_fit.push_back(param[Acts::eBoundQOverP]);
0493     m_eT_fit.push_back(param[Acts::eBoundTime]);
0494
0495     m_res_eLOC0_fit.push_back(res[Acts::eBoundLoc0]);
0496     m_res_eLOC1_fit.push_back(res[Acts::eBoundLoc1]);
0497     m_res_ePHI_fit.push_back(res[Acts::eBoundPhi]);
0498     m_res_eTHETA_fit.push_back(res[Acts::eBoundTheta]);
0499     m_res_eQOP_fit.push_back(res[Acts::eBoundQOverP]);
0500     m_res_eT_fit.push_back(res[Acts::eBoundTime]);
0501
0502     m_err_eLOC0_fit.push_back(error[Acts::eBoundLoc0]);
0503     m_err_eLOC1_fit.push_back(error[Acts::eBoundLoc1]);
0504     m_err_ePHI_fit.push_back(error[Acts::eBoundPhi]);
0505     m_err_eTHETA_fit.push_back(error[Acts::eBoundTheta]);
0506     m_err_eQOP_fit.push_back(error[Acts::eBoundQOverP]);
0507     m_err_eT_fit.push_back(error[Acts::eBoundTime]);
0508
0509     m_pull_eLOC0_fit.push_back(pull[Acts::eBoundLoc0]);
0510     m_pull_eLOC1_fit.push_back(pull[Acts::eBoundLoc1]);
0511     m_pull_ePHI_fit.push_back(pull[Acts::eBoundPhi]);
0512     m_pull_eTHETA_fit.push_back(pull[Acts::eBoundTheta]);
0513     m_pull_eQOP_fit.push_back(pull[Acts::eBoundQOverP]);
0514     m_pull_eT_fit.push_back(pull[Acts::eBoundTime]);
0515
0516     m_hasFittedParams.push_back(hasFittedParams);
0517
0518     if (m_cfg.writeGsfSpecific) {
0519       using namespace Acts::GsfConstants;
0520       if (tracks.hasColumn(Acts::hashString(kFwdMaxMaterialXOverX0))) {
0521         m_gsf_max_material_fwd.push_back(
0522             track.template component<double>(kFwdMaxMaterialXOverX0));
0523       } else {
0524         m_gsf_max_material_fwd.push_back(NaNfloat);
0525       }
0526
0527       if (tracks.hasColumn(Acts::hashString(kFwdSumMaterialXOverX0))) {
0528         m_gsf_sum_material_fwd.push_back(
0529             track.template component<double>(kFwdSumMaterialXOverX0));
0530       } else {
0531         m_gsf_sum_material_fwd.push_back(NaNfloat);
0532       }
0533     }
0534
0535     if (m_cfg.writeCovMat) {
0536       // write all entries of covariance matrix to output file
0537       // one branch for every entry of the matrix.
0538       m_cov_eLOC0_eLOC0.push_back(getCov(0, 0));
0539       m_cov_eLOC0_eLOC1.push_back(getCov(0, 1));
0540       m_cov_eLOC0_ePHI.push_back(getCov(0, 2));
0541       m_cov_eLOC0_eTHETA.push_back(getCov(0, 3));
0542       m_cov_eLOC0_eQOP.push_back(getCov(0, 4));
0543       m_cov_eLOC0_eT.push_back(getCov(0, 5));
0544
0545       m_cov_eLOC1_eLOC0.push_back(getCov(1, 0));
0546       m_cov_eLOC1_eLOC1.push_back(getCov(1, 1));
0547       m_cov_eLOC1_ePHI.push_back(getCov(1, 2));
0548       m_cov_eLOC1_eTHETA.push_back(getCov(1, 3));
0549       m_cov_eLOC1_eQOP.push_back(getCov(1, 4));
0550       m_cov_eLOC1_eT.push_back(getCov(1, 5));
0551
0552       m_cov_ePHI_eLOC0.push_back(getCov(2, 0));
0553       m_cov_ePHI_eLOC1.push_back(getCov(2, 1));
0554       m_cov_ePHI_ePHI.push_back(getCov(2, 2));
0555       m_cov_ePHI_eTHETA.push_back(getCov(2, 3));
0556       m_cov_ePHI_eQOP.push_back(getCov(2, 4));
0557       m_cov_ePHI_eT.push_back(getCov(2, 5));
0558
0559       m_cov_eTHETA_eLOC0.push_back(getCov(3, 0));
0560       m_cov_eTHETA_eLOC1.push_back(getCov(3, 1));
0561       m_cov_eTHETA_ePHI.push_back(getCov(3, 2));
0562       m_cov_eTHETA_eTHETA.push_back(getCov(3, 3));
0563       m_cov_eTHETA_eQOP.push_back(getCov(3, 4));
0564       m_cov_eTHETA_eT.push_back(getCov(3, 5));
0565
0566       m_cov_eQOP_eLOC0.push_back(getCov(4, 0));
0567       m_cov_eQOP_eLOC1.push_back(getCov(4, 1));
0568       m_cov_eQOP_ePHI.push_back(getCov(4, 2));
0569       m_cov_eQOP_eTHETA.push_back(getCov(4, 3));
0570       m_cov_eQOP_eQOP.push_back(getCov(4, 4));
0571       m_cov_eQOP_eT.push_back(getCov(4, 5));
0572
0573       m_cov_eT_eLOC0.push_back(getCov(5, 0));
0574       m_cov_eT_eLOC1.push_back(getCov(5, 1));
0575       m_cov_eT_ePHI.push_back(getCov(5, 2));
0576       m_cov_eT_eTHETA.push_back(getCov(5, 3));
0577       m_cov_eT_eQOP.push_back(getCov(5, 4));
0578       m_cov_eT_eT.push_back(getCov(5, 5));
0579     }
0580
0581     if (m_cfg.writeGx2fSpecific) {
0582       if (tracks.hasColumn(Acts::hashString("Gx2fnUpdateColumn"))) {
0583         int nUpdate = static_cast<int>(
0584             track.template component<std::uint32_t,
0585                                      Acts::hashString("Gx2fnUpdateColumn")>());
0586         m_nUpdatesGx2f.push_back(nUpdate);
0587       } else {
0588         m_nUpdatesGx2f.push_back(-1);
0589       }
0590     }
0591   }
0592
0593   // fill the variables
0594   m_outputTree->Fill();
0595
0596   m_trackNr.clear();
0597   m_nStates.clear();
0598   m_nMeasurements.clear();
0599   m_nOutliers.clear();
0600   m_nHoles.clear();
0601   m_nSharedHits.clear();
0602   m_chi2Sum.clear();
0603   m_NDF.clear();
0604   m_measurementChi2.clear();
0605   m_outlierChi2.clear();
0606   m_measurementVolume.clear();
0607   m_measurementLayer.clear();
0608   m_outlierVolume.clear();
0609   m_outlierLayer.clear();
0610
0611   m_nMajorityHits.clear();
0612   m_majorityParticleVertexPrimary.clear();
0613   m_majorityParticleVertexSecondary.clear();
0614   m_majorityParticleParticle.clear();
0615   m_majorityParticleGeneration.clear();
0616   m_majorityParticleSubParticle.clear();
0617   m_trackClassification.clear();
0618   m_t_charge.clear();
0619   m_t_time.clear();
0620   m_t_vx.clear();
0621   m_t_vy.clear();
0622   m_t_vz.clear();
0623   m_t_px.clear();
0624   m_t_py.clear();
0625   m_t_pz.clear();
0626   m_t_theta.clear();
0627   m_t_phi.clear();
0628   m_t_p.clear();
0629   m_t_pT.clear();
0630   m_t_eta.clear();
0631   m_t_d0.clear();
0632   m_t_z0.clear();
0633   m_t_prodR.clear();
0634
0635   m_hasFittedParams.clear();
0636   m_eLOC0_fit.clear();
0637   m_eLOC1_fit.clear();
0638   m_ePHI_fit.clear();
0639   m_eTHETA_fit.clear();
0640   m_eQOP_fit.clear();
0641   m_eT_fit.clear();
0642   m_err_eLOC0_fit.clear();
0643   m_err_eLOC1_fit.clear();
0644   m_err_ePHI_fit.clear();
0645   m_err_eTHETA_fit.clear();
0646   m_err_eQOP_fit.clear();
0647   m_err_eT_fit.clear();
0648   m_res_eLOC0_fit.clear();
0649   m_res_eLOC1_fit.clear();
0650   m_res_ePHI_fit.clear();
0651   m_res_eTHETA_fit.clear();
0652   m_res_eQOP_fit.clear();
0653   m_res_eT_fit.clear();
0654   m_pull_eLOC0_fit.clear();
0655   m_pull_eLOC1_fit.clear();
0656   m_pull_ePHI_fit.clear();
0657   m_pull_eTHETA_fit.clear();
0658   m_pull_eQOP_fit.clear();
0659   m_pull_eT_fit.clear();
0660
0661   m_gsf_max_material_fwd.clear();
0662   m_gsf_sum_material_fwd.clear();
0663
0664   if (m_cfg.writeCovMat) {
0665     m_cov_eLOC0_eLOC0.clear();
0666     m_cov_eLOC0_eLOC1.clear();
0667     m_cov_eLOC0_ePHI.clear();
0668     m_cov_eLOC0_eTHETA.clear();
0669     m_cov_eLOC0_eQOP.clear();
0670     m_cov_eLOC0_eT.clear();
0671
0672     m_cov_eLOC1_eLOC0.clear();
0673     m_cov_eLOC1_eLOC1.clear();
0674     m_cov_eLOC1_ePHI.clear();
0675     m_cov_eLOC1_eTHETA.clear();
0676     m_cov_eLOC1_eQOP.clear();
0677     m_cov_eLOC1_eT.clear();
0678
0679     m_cov_ePHI_eLOC0.clear();
0680     m_cov_ePHI_eLOC1.clear();
0681     m_cov_ePHI_ePHI.clear();
0682     m_cov_ePHI_eTHETA.clear();
0683     m_cov_ePHI_eQOP.clear();
0684     m_cov_ePHI_eT.clear();
0685
0686     m_cov_eTHETA_eLOC0.clear();
0687     m_cov_eTHETA_eLOC1.clear();
0688     m_cov_eTHETA_ePHI.clear();
0689     m_cov_eTHETA_eTHETA.clear();
0690     m_cov_eTHETA_eQOP.clear();
0691     m_cov_eTHETA_eT.clear();
0692
0693     m_cov_eQOP_eLOC0.clear();
0694     m_cov_eQOP_eLOC1.clear();
0695     m_cov_eQOP_ePHI.clear();
0696     m_cov_eQOP_eTHETA.clear();
0697     m_cov_eQOP_eQOP.clear();
0698     m_cov_eQOP_eT.clear();
0699
0700     m_cov_eT_eLOC0.clear();
0701     m_cov_eT_eLOC1.clear();
0702     m_cov_eT_ePHI.clear();
0703     m_cov_eT_eTHETA.clear();
0704     m_cov_eT_eQOP.clear();
0705     m_cov_eT_eT.clear();
0706   }
0707
0708   m_nUpdatesGx2f.clear();
0709
0710   if (m_cfg.writeJets) {
0711     m_nJets.clear();
0712     m_jet_pt.clear();
0713     m_jet_eta.clear();
0714     m_jet_phi.clear();
0715     m_jet_label.clear();
0716     m_ntracks_per_jets.clear();
0717   }
0718
0719   return ProcessCode::SUCCESS;
0720 }
0721
0722 }  // namespace ActsExamples