Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-10-23 08:23:30

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/VertexNTupleWriter.hpp"
0010 
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Definitions/TrackParametrization.hpp"
0013 #include "Acts/EventData/GenericBoundTrackParameters.hpp"
0014 #include "Acts/EventData/TrackParameters.hpp"
0015 #include "Acts/Propagator/Propagator.hpp"
0016 #include "Acts/Propagator/SympyStepper.hpp"
0017 #include "Acts/Surfaces/PerigeeSurface.hpp"
0018 #include "Acts/Surfaces/Surface.hpp"
0019 #include "Acts/Utilities/Intersection.hpp"
0020 #include "Acts/Utilities/Logger.hpp"
0021 #include "Acts/Utilities/UnitVectors.hpp"
0022 #include "Acts/Vertexing/TrackAtVertex.hpp"
0023 #include "ActsExamples/EventData/SimParticle.hpp"
0024 #include "ActsExamples/EventData/SimVertex.hpp"
0025 #include "ActsExamples/EventData/Track.hpp"
0026 #include "ActsExamples/EventData/TruthMatching.hpp"
0027 #include "ActsFatras/EventData/Barcode.hpp"
0028 
0029 #include <algorithm>
0030 #include <cmath>
0031 #include <cstddef>
0032 #include <cstdint>
0033 #include <ios>
0034 #include <limits>
0035 #include <map>
0036 #include <memory>
0037 #include <numbers>
0038 #include <optional>
0039 #include <ostream>
0040 #include <set>
0041 #include <stdexcept>
0042 #include <utility>
0043 
0044 #include <TFile.h>
0045 #include <TTree.h>
0046 
0047 using Acts::VectorHelpers::eta;
0048 using Acts::VectorHelpers::perp;
0049 using Acts::VectorHelpers::phi;
0050 using Acts::VectorHelpers::theta;
0051 
0052 namespace ActsExamples {
0053 
0054 namespace {
0055 
0056 constexpr double nan = std::numeric_limits<double>::quiet_NaN();
0057 
0058 std::uint32_t getNumberOfReconstructableVertices(
0059     const SimParticleContainer& collection) {
0060   // map for finding frequency
0061   std::map<std::uint32_t, std::uint32_t> fmap;
0062 
0063   std::vector<std::uint32_t> reconstructableTruthVertices;
0064 
0065   // traverse the array for frequency
0066   for (const SimParticle& p : collection) {
0067     std::uint32_t generation = p.particleId().generation();
0068     if (generation > 0) {
0069       // truthparticle from secondary vtx
0070       continue;
0071     }
0072     std::uint32_t priVtxId = p.particleId().vertexPrimary();
0073     fmap[priVtxId]++;
0074   }
0075 
0076   // iterate over the map
0077   for (const auto& [priVtxId, occurrence] : fmap) {
0078     // Require at least 2 tracks
0079     if (occurrence > 1) {
0080       reconstructableTruthVertices.push_back(priVtxId);
0081     }
0082   }
0083 
0084   return reconstructableTruthVertices.size();
0085 }
0086 
0087 std::uint32_t getNumberOfTruePriVertices(
0088     const SimParticleContainer& collection) {
0089   // Vector to store indices of all primary vertices
0090   std::set<std::uint32_t> allPriVtxIds;
0091   for (const SimParticle& p : collection) {
0092     std::uint32_t priVtxId = p.particleId().vertexPrimary();
0093     std::uint32_t generation = p.particleId().generation();
0094     if (generation > 0) {
0095       // truthparticle from secondary vtx
0096       continue;
0097     }
0098     // Insert to set, removing duplicates
0099     allPriVtxIds.insert(priVtxId);
0100   }
0101   // Size of set corresponds to total number of primary vertices
0102   return allPriVtxIds.size();
0103 }
0104 
0105 double calcSumPt2(const VertexNTupleWriter::Config& config,
0106                   const Acts::Vertex& vtx) {
0107   double sumPt2 = 0;
0108   for (const auto& trk : vtx.tracks()) {
0109     if (trk.trackWeight > config.minTrkWeight) {
0110       double pt = trk.originalParams.as<Acts::BoundTrackParameters>()
0111                       ->transverseMomentum();
0112       sumPt2 += pt * pt;
0113     }
0114   }
0115   return sumPt2;
0116 }
0117 
0118 /// Helper function for computing the pull
0119 double pull(double diff, double variance, const std::string& variableStr,
0120             bool afterFit, const Acts::Logger& logger) {
0121   if (variance <= 0) {
0122     std::string tempStr;
0123     if (afterFit) {
0124       tempStr = "after";
0125     } else {
0126       tempStr = "before";
0127     }
0128     ACTS_WARNING("Nonpositive variance " << tempStr << " vertex fit: Var("
0129                                          << variableStr << ") = " << variance
0130                                          << " <= 0.");
0131     return nan;
0132   }
0133   double std = std::sqrt(variance);
0134   return diff / std;
0135 }
0136 
0137 double calculateTruthPrimaryVertexDensity(
0138     const VertexNTupleWriter::Config& config,
0139     const SimVertexContainer& truthVertices, const Acts::Vertex& vtx) {
0140   double z = vtx.fullPosition()[Acts::CoordinateIndices::eZ];
0141   int count = 0;
0142   for (const SimVertex& truthVertex : truthVertices) {
0143     if (truthVertex.vertexId().vertexSecondary() != 0) {
0144       continue;
0145     }
0146     double zTruth = truthVertex.position4[Acts::CoordinateIndices::eZ];
0147     if (std::abs(z - zTruth) <= config.vertexDensityWindow) {
0148       ++count;
0149     }
0150   }
0151   return count / (2 * config.vertexDensityWindow);
0152 }
0153 
0154 const SimParticle* findParticle(
0155     const SimParticleContainer& particles,
0156     const TrackParticleMatching& trackParticleMatching, ConstTrackProxy track,
0157     const Acts::Logger& logger) {
0158   // Get the truth-matched particle
0159   auto imatched = trackParticleMatching.find(track.index());
0160   if (imatched == trackParticleMatching.end() ||
0161       !imatched->second.particle.has_value()) {
0162     ACTS_DEBUG("No truth particle associated with this track, index = "
0163                << track.index() << " tip index = " << track.tipIndex());
0164     return nullptr;
0165   }
0166 
0167   const TrackMatchEntry& particleMatch = imatched->second;
0168 
0169   auto iparticle = particles.find(SimBarcode{particleMatch.particle.value()});
0170   if (iparticle == particles.end()) {
0171     ACTS_DEBUG(
0172         "Truth particle found but not monitored with this track, index = "
0173         << track.index() << " tip index = " << track.tipIndex()
0174         << " and this barcode = " << particleMatch.particle.value());
0175     return {};
0176   }
0177 
0178   return &(*iparticle);
0179 }
0180 
0181 std::optional<ConstTrackProxy> findTrack(const ConstTrackContainer& tracks,
0182                                          const Acts::TrackAtVertex& trkAtVtx) {
0183   // Track parameters before the vertex fit
0184   const Acts::BoundTrackParameters& origTrack =
0185       *trkAtVtx.originalParams.as<Acts::BoundTrackParameters>();
0186 
0187   // Finding the matching parameters in the container of all track
0188   // parameters. This allows us to identify the corresponding particle.
0189   // TODO this should not be necessary if the tracks at vertex would keep
0190   // this information
0191   for (ConstTrackProxy inputTrk : tracks) {
0192     const auto& params = inputTrk.parameters();
0193 
0194     if (origTrack.parameters() == params) {
0195       return inputTrk;
0196     }
0197   }
0198 
0199   return std::nullopt;
0200 }
0201 
0202 }  // namespace
0203 
0204 VertexNTupleWriter::VertexNTupleWriter(const VertexNTupleWriter::Config& config,
0205                                        Acts::Logging::Level level)
0206     : WriterT(config.inputVertices, "VertexNTupleWriter", level),
0207       m_cfg(config) {
0208   if (m_cfg.filePath.empty()) {
0209     throw std::invalid_argument("Missing output filename");
0210   }
0211   if (m_cfg.treeName.empty()) {
0212     throw std::invalid_argument("Missing tree name");
0213   }
0214   if (m_cfg.inputTruthVertices.empty()) {
0215     throw std::invalid_argument("Collection with truth vertices missing");
0216   }
0217   if (m_cfg.inputParticles.empty()) {
0218     throw std::invalid_argument("Collection with particles missing");
0219   }
0220   if (m_cfg.inputSelectedParticles.empty()) {
0221     throw std::invalid_argument("Collection with selected particles missing");
0222   }
0223 
0224   m_inputTracks.maybeInitialize(m_cfg.inputTracks);
0225   m_inputTruthVertices.initialize(m_cfg.inputTruthVertices);
0226   m_inputParticles.initialize(m_cfg.inputParticles);
0227   m_inputSelectedParticles.initialize(m_cfg.inputSelectedParticles);
0228   m_inputTrackParticleMatching.maybeInitialize(
0229       m_cfg.inputTrackParticleMatching);
0230 
0231   if (m_cfg.writeTrackInfo && !m_inputTracks.isInitialized()) {
0232     throw std::invalid_argument("Missing input tracks collection");
0233   }
0234   if (m_cfg.writeTrackInfo && !m_inputTrackParticleMatching.isInitialized()) {
0235     throw std::invalid_argument("Missing input track particles matching");
0236   }
0237 
0238   // Setup ROOT I/O
0239   m_outputFile = TFile::Open(m_cfg.filePath.c_str(), m_cfg.fileMode.c_str());
0240   if (m_outputFile == nullptr) {
0241     throw std::ios_base::failure("Could not open '" + m_cfg.filePath + "'");
0242   }
0243   m_outputFile->cd();
0244   m_outputTree = new TTree(m_cfg.treeName.c_str(), m_cfg.treeName.c_str());
0245   if (m_outputTree == nullptr) {
0246     throw std::bad_alloc();
0247   }
0248 
0249   m_outputTree->Branch("event_nr", &m_eventNr);
0250 
0251   m_outputTree->Branch("nRecoVtx", &m_nRecoVtx);
0252   m_outputTree->Branch("nTrueVtx", &m_nTrueVtx);
0253   m_outputTree->Branch("nCleanVtx", &m_nCleanVtx);
0254   m_outputTree->Branch("nMergedVtx", &m_nMergedVtx);
0255   m_outputTree->Branch("nSplitVtx", &m_nSplitVtx);
0256   m_outputTree->Branch("nVtxDetectorAcceptance", &m_nVtxDetAcceptance);
0257   m_outputTree->Branch("nVtxReconstructable", &m_nVtxReconstructable);
0258 
0259   m_outputTree->Branch("nTracksRecoVtx", &m_nTracksOnRecoVertex);
0260 
0261   m_outputTree->Branch("recoVertexTrackWeights", &m_recoVertexTrackWeights);
0262 
0263   m_outputTree->Branch("sumPt2", &m_sumPt2);
0264 
0265   m_outputTree->Branch("recoX", &m_recoX);
0266   m_outputTree->Branch("recoY", &m_recoY);
0267   m_outputTree->Branch("recoZ", &m_recoZ);
0268   m_outputTree->Branch("recoT", &m_recoT);
0269 
0270   m_outputTree->Branch("covXX", &m_covXX);
0271   m_outputTree->Branch("covYY", &m_covYY);
0272   m_outputTree->Branch("covZZ", &m_covZZ);
0273   m_outputTree->Branch("covTT", &m_covTT);
0274   m_outputTree->Branch("covXY", &m_covXY);
0275   m_outputTree->Branch("covXZ", &m_covXZ);
0276   m_outputTree->Branch("covXT", &m_covXT);
0277   m_outputTree->Branch("covYZ", &m_covYZ);
0278   m_outputTree->Branch("covYT", &m_covYT);
0279   m_outputTree->Branch("covZT", &m_covZT);
0280 
0281   m_outputTree->Branch("seedX", &m_seedX);
0282   m_outputTree->Branch("seedY", &m_seedY);
0283   m_outputTree->Branch("seedZ", &m_seedZ);
0284   m_outputTree->Branch("seedT", &m_seedT);
0285 
0286   m_outputTree->Branch("vertex_primary", &m_vertexPrimary);
0287   m_outputTree->Branch("vertex_secondary", &m_vertexSecondary);
0288 
0289   m_outputTree->Branch("nTracksTruthVtx", &m_nTracksOnTruthVertex);
0290 
0291   m_outputTree->Branch("truthPrimaryVertexDensity",
0292                        &m_truthPrimaryVertexDensity);
0293 
0294   m_outputTree->Branch("truthVertexTrackWeights", &m_truthVertexTrackWeights);
0295   m_outputTree->Branch("truthVertexMatchRatio", &m_truthVertexMatchRatio);
0296   m_outputTree->Branch("recoVertexContamination", &m_recoVertexContamination);
0297 
0298   m_outputTree->Branch("recoVertexClassification", &m_recoVertexClassification);
0299 
0300   m_outputTree->Branch("truthX", &m_truthX);
0301   m_outputTree->Branch("truthY", &m_truthY);
0302   m_outputTree->Branch("truthZ", &m_truthZ);
0303   m_outputTree->Branch("truthT", &m_truthT);
0304 
0305   m_outputTree->Branch("resX", &m_resX);
0306   m_outputTree->Branch("resY", &m_resY);
0307   m_outputTree->Branch("resZ", &m_resZ);
0308   m_outputTree->Branch("resT", &m_resT);
0309 
0310   m_outputTree->Branch("resSeedZ", &m_resSeedZ);
0311   m_outputTree->Branch("resSeedT", &m_resSeedT);
0312 
0313   m_outputTree->Branch("pullX", &m_pullX);
0314   m_outputTree->Branch("pullY", &m_pullY);
0315   m_outputTree->Branch("pullZ", &m_pullZ);
0316   m_outputTree->Branch("pullT", &m_pullT);
0317 
0318   if (m_cfg.writeTrackInfo) {
0319     m_outputTree->Branch("trk_weight", &m_trkWeight);
0320 
0321     m_outputTree->Branch("trk_recoPhi", &m_recoPhi);
0322     m_outputTree->Branch("trk_recoTheta", &m_recoTheta);
0323     m_outputTree->Branch("trk_recoQOverP", &m_recoQOverP);
0324 
0325     m_outputTree->Branch("trk_recoPhiFitted", &m_recoPhiFitted);
0326     m_outputTree->Branch("trk_recoThetaFitted", &m_recoThetaFitted);
0327     m_outputTree->Branch("trk_recoQOverPFitted", &m_recoQOverPFitted);
0328 
0329     m_outputTree->Branch("trk_particleId", &m_trkParticleId);
0330 
0331     m_outputTree->Branch("trk_truthPhi", &m_truthPhi);
0332     m_outputTree->Branch("trk_truthTheta", &m_truthTheta);
0333     m_outputTree->Branch("trk_truthQOverP", &m_truthQOverP);
0334 
0335     m_outputTree->Branch("trk_resPhi", &m_resPhi);
0336     m_outputTree->Branch("trk_resTheta", &m_resTheta);
0337     m_outputTree->Branch("trk_resQOverP", &m_resQOverP);
0338     m_outputTree->Branch("trk_momOverlap", &m_momOverlap);
0339 
0340     m_outputTree->Branch("trk_resPhiFitted", &m_resPhiFitted);
0341     m_outputTree->Branch("trk_resThetaFitted", &m_resThetaFitted);
0342     m_outputTree->Branch("trk_resQOverPFitted", &m_resQOverPFitted);
0343     m_outputTree->Branch("trk_momOverlapFitted", &m_momOverlapFitted);
0344 
0345     m_outputTree->Branch("trk_pullPhi", &m_pullPhi);
0346     m_outputTree->Branch("trk_pullTheta", &m_pullTheta);
0347     m_outputTree->Branch("trk_pullQOverP", &m_pullQOverP);
0348 
0349     m_outputTree->Branch("trk_pullPhiFitted", &m_pullPhiFitted);
0350     m_outputTree->Branch("trk_pullThetaFitted", &m_pullThetaFitted);
0351     m_outputTree->Branch("trk_pullQOverPFitted", &m_pullQOverPFitted);
0352   }
0353 }
0354 
0355 VertexNTupleWriter::~VertexNTupleWriter() {
0356   if (m_outputFile != nullptr) {
0357     m_outputFile->Close();
0358   }
0359 }
0360 
0361 ProcessCode VertexNTupleWriter::finalize() {
0362   m_outputFile->cd();
0363   m_outputTree->Write();
0364   m_outputFile->Close();
0365 
0366   return ProcessCode::SUCCESS;
0367 }
0368 
0369 ProcessCode VertexNTupleWriter::writeT(
0370     const AlgorithmContext& ctx, const std::vector<Acts::Vertex>& vertices) {
0371   // In case we do not have any tracks in the vertex, we create empty
0372   // collections
0373   const static TrackParticleMatching emptyTrackParticleMatching;
0374   const static Acts::ConstVectorTrackContainer emptyConstTrackContainer;
0375   const static Acts::ConstVectorMultiTrajectory emptyConstTrackStateContainer;
0376   const static ConstTrackContainer emptyTracks(
0377       std::make_shared<Acts::ConstVectorTrackContainer>(
0378           emptyConstTrackContainer),
0379       std::make_shared<Acts::ConstVectorMultiTrajectory>(
0380           emptyConstTrackStateContainer));
0381 
0382   // Read truth vertex input collection
0383   const SimVertexContainer& truthVertices = m_inputTruthVertices(ctx);
0384   // Read truth particle input collection
0385   const SimParticleContainer& particles = m_inputParticles(ctx);
0386   const SimParticleContainer& selectedParticles = m_inputSelectedParticles(ctx);
0387   const TrackParticleMatching& trackParticleMatching =
0388       (m_inputTrackParticleMatching.isInitialized()
0389            ? m_inputTrackParticleMatching(ctx)
0390            : emptyTrackParticleMatching);
0391   const ConstTrackContainer& tracks =
0392       (m_inputTracks.isInitialized() ? m_inputTracks(ctx) : emptyTracks);
0393   SimParticleContainer recoParticles;
0394 
0395   for (ConstTrackProxy track : tracks) {
0396     if (!track.hasReferenceSurface()) {
0397       ACTS_DEBUG("No reference surface on this track, index = "
0398                  << track.index() << " tip index = " << track.tipIndex());
0399       continue;
0400     }
0401 
0402     if (const SimParticle* particle =
0403             findParticle(particles, trackParticleMatching, track, logger());
0404         particle != nullptr) {
0405       recoParticles.insert(*particle);
0406     }
0407   }
0408   if (tracks.size() == 0) {
0409     // if not using tracks, then all truth particles are associated with the
0410     // vertex
0411     recoParticles = particles;
0412   }
0413 
0414   // Exclusive access to the tree while writing
0415   std::lock_guard<std::mutex> lock(m_writeMutex);
0416 
0417   m_nRecoVtx = vertices.size();
0418   m_nCleanVtx = 0;
0419   m_nMergedVtx = 0;
0420   m_nSplitVtx = 0;
0421 
0422   ACTS_DEBUG("Number of reco vertices in event: " << m_nRecoVtx);
0423 
0424   // Get number of generated true primary vertices
0425   m_nTrueVtx = getNumberOfTruePriVertices(particles);
0426   // Get number of detector-accepted true primary vertices
0427   m_nVtxDetAcceptance = getNumberOfTruePriVertices(selectedParticles);
0428 
0429   ACTS_DEBUG("Number of truth particles in event : " << particles.size());
0430   ACTS_DEBUG("Number of truth primary vertices : " << m_nTrueVtx);
0431   ACTS_DEBUG("Number of detector-accepted truth primary vertices : "
0432              << m_nVtxDetAcceptance);
0433 
0434   // Get the event number
0435   m_eventNr = ctx.eventNumber;
0436 
0437   // Get number of track-associated true primary vertices
0438   m_nVtxReconstructable = getNumberOfReconstructableVertices(recoParticles);
0439 
0440   ACTS_DEBUG("Number of reconstructed tracks : " << tracks.size());
0441   ACTS_DEBUG("Number of reco track-associated truth particles in event : "
0442              << recoParticles.size());
0443   ACTS_DEBUG("Maximum number of reconstructible primary vertices : "
0444              << m_nVtxReconstructable);
0445 
0446   struct ToTruthMatching {
0447     std::optional<SimVertexBarcode> vertexId;
0448     double totalTrackWeight{};
0449     double truthMajorityTrackWeights{};
0450     double matchFraction{};
0451 
0452     RecoVertexClassification classification{RecoVertexClassification::Unknown};
0453   };
0454   struct ToRecoMatching {
0455     std::size_t recoIndex{};
0456 
0457     double recoSumPt2{};
0458   };
0459 
0460   std::vector<ToTruthMatching> recoToTruthMatching;
0461   std::map<SimVertexBarcode, ToRecoMatching> truthToRecoMatching;
0462 
0463   // Do truth matching for each reconstructed vertex
0464   for (const auto& [vtxIndex, vtx] : Acts::enumerate(vertices)) {
0465     // Reconstructed tracks that contribute to the reconstructed vertex
0466     const std::vector<Acts::TrackAtVertex>& tracksAtVtx = vtx.tracks();
0467 
0468     // Containers for storing truth particles and truth vertices that
0469     // contribute to the reconstructed vertex
0470     std::vector<std::pair<SimVertexBarcode, double>> contributingTruthVertices;
0471 
0472     double totalTrackWeight = 0;
0473     if (!tracksAtVtx.empty()) {
0474       for (const Acts::TrackAtVertex& trk : tracksAtVtx) {
0475         if (trk.trackWeight < m_cfg.minTrkWeight) {
0476           continue;
0477         }
0478 
0479         totalTrackWeight += trk.trackWeight;
0480 
0481         std::optional<ConstTrackProxy> trackOpt = findTrack(tracks, trk);
0482         if (!trackOpt.has_value()) {
0483           ACTS_DEBUG("Track has no matching input track.");
0484           continue;
0485         }
0486         const ConstTrackProxy& inputTrk = *trackOpt;
0487         const SimParticle* particle =
0488             findParticle(particles, trackParticleMatching, inputTrk, logger());
0489         if (particle == nullptr) {
0490           ACTS_VERBOSE("Track has no matching truth particle.");
0491         } else {
0492           contributingTruthVertices.emplace_back(
0493               SimBarcode{particle->particleId()}.vertexId(), trk.trackWeight);
0494         }
0495       }
0496     } else {
0497       // If no tracks at the vertex, then use all truth particles
0498       for (auto& particle : recoParticles) {
0499         contributingTruthVertices.emplace_back(
0500             SimBarcode{particle.particleId()}.vertexId(), 1.);
0501       }
0502     }
0503 
0504     // Find true vertex that contributes most to the reconstructed vertex
0505     std::map<SimVertexBarcode, std::pair<int, double>> fmap;
0506     for (const auto& [vtxId, weight] : contributingTruthVertices) {
0507       ++fmap[vtxId].first;
0508       fmap[vtxId].second += weight;
0509     }
0510     double truthMajorityVertexTrackWeights = 0;
0511     std::optional<SimVertexBarcode> truthMajorityVertexId = std::nullopt;
0512     for (const auto& [vtxId, counter] : fmap) {
0513       if (counter.second > truthMajorityVertexTrackWeights) {
0514         truthMajorityVertexId = vtxId;
0515         truthMajorityVertexTrackWeights = counter.second;
0516       }
0517     }
0518 
0519     double sumPt2 = calcSumPt2(m_cfg, vtx);
0520 
0521     double vertexMatchFraction =
0522         (totalTrackWeight > 0
0523              ? truthMajorityVertexTrackWeights / totalTrackWeight
0524              : 0.);
0525     RecoVertexClassification recoVertexClassification =
0526         RecoVertexClassification::Unknown;
0527 
0528     if (vertexMatchFraction >= m_cfg.vertexMatchThreshold) {
0529       recoVertexClassification = RecoVertexClassification::Clean;
0530     } else {
0531       recoVertexClassification = RecoVertexClassification::Merged;
0532     }
0533 
0534     recoToTruthMatching.push_back({truthMajorityVertexId, totalTrackWeight,
0535                                    truthMajorityVertexTrackWeights,
0536                                    vertexMatchFraction,
0537                                    recoVertexClassification});
0538 
0539     auto& recoToTruth = recoToTruthMatching.back();
0540 
0541     // We have to decide if this reco vertex is a split vertex.
0542     if (!truthMajorityVertexId.has_value()) {
0543       // No truth vertex matched to this reconstructed vertex
0544       ACTS_DEBUG("No truth vertex matched to this reconstructed vertex.");
0545       continue;
0546     }
0547 
0548     if (auto it = truthToRecoMatching.find(truthMajorityVertexId.value());
0549         it != truthToRecoMatching.end()) {
0550       // This truth vertex is already matched to a reconstructed vertex so we
0551       // are dealing with a split vertex.
0552 
0553       // We have to decide which of the two reconstructed vertices is the
0554       // split vertex.
0555       if (sumPt2 <= it->second.recoSumPt2) {
0556         // Since the sumPt2 is smaller we can simply call this a split vertex
0557 
0558         recoToTruth.classification = RecoVertexClassification::Split;
0559 
0560         // Keep the existing truth to reco matching
0561       } else {
0562         // The sumPt2 is larger, so we call the other vertex a split vertex.
0563 
0564         auto& otherRecoToTruth = recoToTruthMatching.at(it->second.recoIndex);
0565         // Swap the classification
0566         recoToTruth.classification = otherRecoToTruth.classification;
0567         otherRecoToTruth.classification = RecoVertexClassification::Split;
0568 
0569         // Overwrite the truth to reco matching
0570         it->second = {vtxIndex, sumPt2};
0571       }
0572     } else {
0573       truthToRecoMatching[truthMajorityVertexId.value()] = {vtxIndex, sumPt2};
0574     }
0575   }
0576 
0577   // Loop over reconstructed vertices and see if they can be matched to a true
0578   // vertex.
0579   for (const auto& [vtxIndex, vtx] : Acts::enumerate(vertices)) {
0580     // Reconstructed tracks that contribute to the reconstructed vertex
0581     const auto& tracksAtVtx = vtx.tracks();
0582     // Input tracks matched to `tracksAtVtx`
0583     std::vector<std::uint32_t> trackIndices;
0584 
0585     const auto& toTruthMatching = recoToTruthMatching[vtxIndex];
0586 
0587     m_recoX.push_back(vtx.fullPosition()[Acts::CoordinateIndices::eX]);
0588     m_recoY.push_back(vtx.fullPosition()[Acts::CoordinateIndices::eY]);
0589     m_recoZ.push_back(vtx.fullPosition()[Acts::CoordinateIndices::eZ]);
0590     m_recoT.push_back(vtx.fullPosition()[Acts::CoordinateIndices::eTime]);
0591 
0592     double varX = vtx.fullCovariance()(Acts::CoordinateIndices::eX,
0593                                        Acts::CoordinateIndices::eX);
0594     double varY = vtx.fullCovariance()(Acts::CoordinateIndices::eY,
0595                                        Acts::CoordinateIndices::eY);
0596     double varZ = vtx.fullCovariance()(Acts::CoordinateIndices::eZ,
0597                                        Acts::CoordinateIndices::eZ);
0598     double varTime = vtx.fullCovariance()(Acts::CoordinateIndices::eTime,
0599                                           Acts::CoordinateIndices::eTime);
0600 
0601     m_covXX.push_back(varX);
0602     m_covYY.push_back(varY);
0603     m_covZZ.push_back(varZ);
0604     m_covTT.push_back(varTime);
0605     m_covXY.push_back(vtx.fullCovariance()(Acts::CoordinateIndices::eX,
0606                                            Acts::CoordinateIndices::eY));
0607     m_covXZ.push_back(vtx.fullCovariance()(Acts::CoordinateIndices::eX,
0608                                            Acts::CoordinateIndices::eZ));
0609     m_covXT.push_back(vtx.fullCovariance()(Acts::CoordinateIndices::eX,
0610                                            Acts::CoordinateIndices::eTime));
0611     m_covYZ.push_back(vtx.fullCovariance()(Acts::CoordinateIndices::eY,
0612                                            Acts::CoordinateIndices::eZ));
0613     m_covYT.push_back(vtx.fullCovariance()(Acts::CoordinateIndices::eY,
0614                                            Acts::CoordinateIndices::eTime));
0615     m_covZT.push_back(vtx.fullCovariance()(Acts::CoordinateIndices::eZ,
0616                                            Acts::CoordinateIndices::eTime));
0617 
0618     double sumPt2 = calcSumPt2(m_cfg, vtx);
0619     m_sumPt2.push_back(sumPt2);
0620 
0621     double recoVertexTrackWeights = 0;
0622     for (const Acts::TrackAtVertex& trk : tracksAtVtx) {
0623       if (trk.trackWeight < m_cfg.minTrkWeight) {
0624         continue;
0625       }
0626       recoVertexTrackWeights += trk.trackWeight;
0627     }
0628     m_recoVertexTrackWeights.push_back(recoVertexTrackWeights);
0629 
0630     unsigned int nTracksOnRecoVertex = std::count_if(
0631         tracksAtVtx.begin(), tracksAtVtx.end(), [this](const auto& trkAtVtx) {
0632           return trkAtVtx.trackWeight > m_cfg.minTrkWeight;
0633         });
0634     m_nTracksOnRecoVertex.push_back(nTracksOnRecoVertex);
0635 
0636     // Saving truth information for the reconstructed vertex
0637     bool truthInfoWritten = false;
0638     std::optional<Acts::Vector4> truthPos;
0639     if (toTruthMatching.vertexId.has_value()) {
0640       auto iTruthVertex = truthVertices.find(toTruthMatching.vertexId.value());
0641       if (iTruthVertex == truthVertices.end()) {
0642         ACTS_ERROR("Truth vertex not found for id: "
0643                    << toTruthMatching.vertexId.value());
0644         continue;
0645       }
0646       const SimVertex& truthVertex = *iTruthVertex;
0647 
0648       m_vertexPrimary.push_back(truthVertex.vertexId().vertexPrimary());
0649       m_vertexSecondary.push_back(truthVertex.vertexId().vertexSecondary());
0650 
0651       // Count number of reconstructible tracks on truth vertex
0652       int nTracksOnTruthVertex = 0;
0653       for (const auto& particle : selectedParticles) {
0654         if (static_cast<SimVertexBarcode>(particle.particleId().vertexId()) ==
0655             truthVertex.vertexId()) {
0656           ++nTracksOnTruthVertex;
0657         }
0658       }
0659       m_nTracksOnTruthVertex.push_back(nTracksOnTruthVertex);
0660 
0661       double truthPrimaryVertexDensity =
0662           calculateTruthPrimaryVertexDensity(m_cfg, truthVertices, vtx);
0663       m_truthPrimaryVertexDensity.push_back(truthPrimaryVertexDensity);
0664 
0665       double truthVertexTrackWeights =
0666           toTruthMatching.truthMajorityTrackWeights;
0667       m_truthVertexTrackWeights.push_back(truthVertexTrackWeights);
0668 
0669       double truthVertexMatchRatio = toTruthMatching.matchFraction;
0670       m_truthVertexMatchRatio.push_back(truthVertexMatchRatio);
0671 
0672       double recoVertexContamination = 1 - truthVertexMatchRatio;
0673       m_recoVertexContamination.push_back(recoVertexContamination);
0674 
0675       RecoVertexClassification recoVertexClassification =
0676           toTruthMatching.classification;
0677       m_recoVertexClassification.push_back(
0678           static_cast<int>(recoVertexClassification));
0679 
0680       if (recoVertexClassification == RecoVertexClassification::Clean) {
0681         ++m_nCleanVtx;
0682       } else if (recoVertexClassification == RecoVertexClassification::Merged) {
0683         ++m_nMergedVtx;
0684       } else if (recoVertexClassification == RecoVertexClassification::Split) {
0685         ++m_nSplitVtx;
0686       }
0687 
0688       const Acts::Vector4& truePos = truthVertex.position4;
0689       truthPos = truePos;
0690       m_truthX.push_back(truePos[Acts::CoordinateIndices::eX]);
0691       m_truthY.push_back(truePos[Acts::CoordinateIndices::eY]);
0692       m_truthZ.push_back(truePos[Acts::CoordinateIndices::eZ]);
0693       m_truthT.push_back(truePos[Acts::CoordinateIndices::eTime]);
0694 
0695       const Acts::Vector4 diffPos = vtx.fullPosition() - truePos;
0696       m_resX.push_back(diffPos[Acts::CoordinateIndices::eX]);
0697       m_resY.push_back(diffPos[Acts::CoordinateIndices::eY]);
0698       m_resZ.push_back(diffPos[Acts::CoordinateIndices::eZ]);
0699       m_resT.push_back(diffPos[Acts::CoordinateIndices::eTime]);
0700 
0701       m_pullX.push_back(pull(diffPos[Acts::CoordinateIndices::eX], varX, "X",
0702                              true, logger()));
0703       m_pullY.push_back(pull(diffPos[Acts::CoordinateIndices::eY], varY, "Y",
0704                              true, logger()));
0705       m_pullZ.push_back(pull(diffPos[Acts::CoordinateIndices::eZ], varZ, "Z",
0706                              true, logger()));
0707       m_pullT.push_back(pull(diffPos[Acts::CoordinateIndices::eTime], varTime,
0708                              "T", true, logger()));
0709 
0710       truthInfoWritten = true;
0711     }
0712     if (!truthInfoWritten) {
0713       m_vertexPrimary.push_back(-1);
0714       m_vertexSecondary.push_back(-1);
0715 
0716       m_nTracksOnTruthVertex.push_back(-1);
0717 
0718       m_truthPrimaryVertexDensity.push_back(nan);
0719 
0720       m_truthVertexTrackWeights.push_back(nan);
0721       m_truthVertexMatchRatio.push_back(nan);
0722       m_recoVertexContamination.push_back(nan);
0723 
0724       m_recoVertexClassification.push_back(
0725           static_cast<int>(RecoVertexClassification::Unknown));
0726 
0727       m_truthX.push_back(nan);
0728       m_truthY.push_back(nan);
0729       m_truthZ.push_back(nan);
0730       m_truthT.push_back(nan);
0731 
0732       m_resX.push_back(nan);
0733       m_resY.push_back(nan);
0734       m_resZ.push_back(nan);
0735       m_resT.push_back(nan);
0736 
0737       m_pullX.push_back(nan);
0738       m_pullY.push_back(nan);
0739       m_pullZ.push_back(nan);
0740       m_pullT.push_back(nan);
0741     }
0742 
0743     if (m_cfg.writeTrackInfo) {
0744       writeTrackInfo(ctx, particles, tracks, trackParticleMatching, truthPos,
0745                      tracksAtVtx);
0746     }
0747   }
0748 
0749   // fill the variables
0750   m_outputTree->Fill();
0751 
0752   m_nTracksOnRecoVertex.clear();
0753   m_recoVertexTrackWeights.clear();
0754   m_recoX.clear();
0755   m_recoY.clear();
0756   m_recoZ.clear();
0757   m_recoT.clear();
0758   m_covXX.clear();
0759   m_covYY.clear();
0760   m_covZZ.clear();
0761   m_covTT.clear();
0762   m_covXY.clear();
0763   m_covXZ.clear();
0764   m_covXT.clear();
0765   m_covYZ.clear();
0766   m_covYT.clear();
0767   m_covZT.clear();
0768   m_seedX.clear();
0769   m_seedY.clear();
0770   m_seedZ.clear();
0771   m_seedT.clear();
0772   m_vertexPrimary.clear();
0773   m_vertexSecondary.clear();
0774   m_nTracksOnTruthVertex.clear();
0775   m_truthPrimaryVertexDensity.clear();
0776   m_truthVertexTrackWeights.clear();
0777   m_truthVertexMatchRatio.clear();
0778   m_recoVertexContamination.clear();
0779   m_recoVertexClassification.clear();
0780   m_truthX.clear();
0781   m_truthY.clear();
0782   m_truthZ.clear();
0783   m_truthT.clear();
0784   m_resX.clear();
0785   m_resY.clear();
0786   m_resZ.clear();
0787   m_resT.clear();
0788   m_resSeedZ.clear();
0789   m_resSeedT.clear();
0790   m_pullX.clear();
0791   m_pullY.clear();
0792   m_pullZ.clear();
0793   m_pullT.clear();
0794   m_sumPt2.clear();
0795   m_trkWeight.clear();
0796   m_recoPhi.clear();
0797   m_recoTheta.clear();
0798   m_recoQOverP.clear();
0799   m_recoPhiFitted.clear();
0800   m_recoThetaFitted.clear();
0801   m_recoQOverPFitted.clear();
0802   m_trkParticleId.clear();
0803   m_truthPhi.clear();
0804   m_truthTheta.clear();
0805   m_truthQOverP.clear();
0806   m_resPhi.clear();
0807   m_resTheta.clear();
0808   m_resQOverP.clear();
0809   m_momOverlap.clear();
0810   m_resPhiFitted.clear();
0811   m_resThetaFitted.clear();
0812   m_resQOverPFitted.clear();
0813   m_momOverlapFitted.clear();
0814   m_pullPhi.clear();
0815   m_pullTheta.clear();
0816   m_pullQOverP.clear();
0817   m_pullPhiFitted.clear();
0818   m_pullThetaFitted.clear();
0819   m_pullQOverPFitted.clear();
0820 
0821   return ProcessCode::SUCCESS;
0822 }
0823 
0824 void VertexNTupleWriter::writeTrackInfo(
0825     const AlgorithmContext& ctx, const SimParticleContainer& particles,
0826     const ConstTrackContainer& tracks,
0827     const TrackParticleMatching& trackParticleMatching,
0828     const std::optional<Acts::Vector4>& truthPos,
0829     const std::vector<Acts::TrackAtVertex>& tracksAtVtx) {
0830   // We compare the reconstructed momenta to the true momenta at the vertex. For
0831   // this, we propagate the reconstructed tracks to the PCA of the true vertex
0832   // position. Setting up propagator:
0833   Acts::SympyStepper stepper(m_cfg.bField);
0834   using Propagator = Acts::Propagator<Acts::SympyStepper>;
0835   auto propagator = std::make_shared<Propagator>(stepper);
0836 
0837   // Saving the reconstructed/truth momenta. The reconstructed momenta
0838   // are taken at the PCA to the truth vertex position -> we need to
0839   // perform a propagation.
0840 
0841   // Get references to inner vectors where all track variables corresponding
0842   // to the current vertex will be saved
0843   auto& innerTrkWeight = m_trkWeight.emplace_back();
0844 
0845   auto& innerRecoPhi = m_recoPhi.emplace_back();
0846   auto& innerRecoTheta = m_recoTheta.emplace_back();
0847   auto& innerRecoQOverP = m_recoQOverP.emplace_back();
0848 
0849   auto& innerRecoPhiFitted = m_recoPhiFitted.emplace_back();
0850   auto& innerRecoThetaFitted = m_recoThetaFitted.emplace_back();
0851   auto& innerRecoQOverPFitted = m_recoQOverPFitted.emplace_back();
0852 
0853   auto& innerTrkParticleId = m_trkParticleId.emplace_back();
0854 
0855   auto& innerTruthPhi = m_truthPhi.emplace_back();
0856   auto& innerTruthTheta = m_truthTheta.emplace_back();
0857   auto& innerTruthQOverP = m_truthQOverP.emplace_back();
0858 
0859   auto& innerResPhi = m_resPhi.emplace_back();
0860   auto& innerResTheta = m_resTheta.emplace_back();
0861   auto& innerResQOverP = m_resQOverP.emplace_back();
0862 
0863   auto& innerResPhiFitted = m_resPhiFitted.emplace_back();
0864   auto& innerResThetaFitted = m_resThetaFitted.emplace_back();
0865   auto& innerResQOverPFitted = m_resQOverPFitted.emplace_back();
0866 
0867   auto& innerMomOverlap = m_momOverlap.emplace_back();
0868   auto& innerMomOverlapFitted = m_momOverlapFitted.emplace_back();
0869 
0870   auto& innerPullPhi = m_pullPhi.emplace_back();
0871   auto& innerPullTheta = m_pullTheta.emplace_back();
0872   auto& innerPullQOverP = m_pullQOverP.emplace_back();
0873 
0874   auto& innerPullPhiFitted = m_pullPhiFitted.emplace_back();
0875   auto& innerPullThetaFitted = m_pullThetaFitted.emplace_back();
0876   auto& innerPullQOverPFitted = m_pullQOverPFitted.emplace_back();
0877 
0878   // Perigee at the true vertex position
0879   std::shared_ptr<Acts::PerigeeSurface> perigeeSurface;
0880   if (truthPos.has_value()) {
0881     perigeeSurface =
0882         Acts::Surface::makeShared<Acts::PerigeeSurface>(truthPos->head<3>());
0883   }
0884   // Lambda for propagating the tracks to the PCA
0885   auto propagateToVtx =
0886       [&](const auto& params) -> std::optional<Acts::BoundTrackParameters> {
0887     if (!perigeeSurface) {
0888       return std::nullopt;
0889     }
0890 
0891     Acts::Intersection3D intersection =
0892         perigeeSurface
0893             ->intersect(ctx.geoContext, params.position(ctx.geoContext),
0894                         params.direction(), Acts::BoundaryTolerance::Infinite())
0895             .closest();
0896 
0897     // Setting the geometry/magnetic field context for the event
0898     using PropagatorOptions = Propagator::Options<>;
0899     PropagatorOptions pOptions(ctx.geoContext, ctx.magFieldContext);
0900     pOptions.direction =
0901         Acts::Direction::fromScalarZeroAsPositive(intersection.pathLength());
0902 
0903     auto result = propagator->propagate(params, *perigeeSurface, pOptions);
0904     if (!result.ok()) {
0905       ACTS_ERROR("Propagation to true vertex position failed.");
0906       return std::nullopt;
0907     }
0908     auto& paramsAtVtx = *result->endParameters;
0909     return std::make_optional(paramsAtVtx);
0910   };
0911 
0912   for (const Acts::TrackAtVertex& trk : tracksAtVtx) {
0913     if (trk.trackWeight < m_cfg.minTrkWeight) {
0914       continue;
0915     }
0916 
0917     innerTrkWeight.push_back(trk.trackWeight);
0918 
0919     Acts::Vector3 trueUnitDir = Acts::Vector3::Zero();
0920     Acts::Vector3 trueMom = Acts::Vector3::Zero();
0921 
0922     const SimParticle* particle = nullptr;
0923     std::optional<ConstTrackProxy> trackOpt = findTrack(tracks, trk);
0924     if (trackOpt.has_value()) {
0925       particle =
0926           findParticle(particles, trackParticleMatching, *trackOpt, logger());
0927     }
0928 
0929     if (particle != nullptr) {
0930       innerTrkParticleId.push_back(particle->particleId().asVector());
0931 
0932       trueUnitDir = particle->direction();
0933       trueMom.head<2>() = Acts::makePhiThetaFromDirection(trueUnitDir);
0934       trueMom[2] = particle->qOverP();
0935 
0936       innerTruthPhi.push_back(trueMom[0]);
0937       innerTruthTheta.push_back(trueMom[1]);
0938       innerTruthQOverP.push_back(trueMom[2]);
0939     } else {
0940       ACTS_VERBOSE("Track has no matching truth particle.");
0941 
0942       innerTrkParticleId.push_back(ActsFatras::Barcode::Invalid().asVector());
0943 
0944       innerTruthPhi.push_back(nan);
0945       innerTruthTheta.push_back(nan);
0946       innerTruthQOverP.push_back(nan);
0947     }
0948 
0949     // Save track parameters before the vertex fit
0950     const auto paramsAtVtx =
0951         propagateToVtx(*(trk.originalParams.as<Acts::BoundTrackParameters>()));
0952     if (paramsAtVtx.has_value()) {
0953       Acts::Vector3 recoMom =
0954           paramsAtVtx->parameters().segment(Acts::eBoundPhi, 3);
0955       const Acts::SquareMatrix3& momCov =
0956           paramsAtVtx->covariance()->template block<3, 3>(Acts::eBoundPhi,
0957                                                           Acts::eBoundPhi);
0958       innerRecoPhi.push_back(recoMom[0]);
0959       innerRecoTheta.push_back(recoMom[1]);
0960       innerRecoQOverP.push_back(recoMom[2]);
0961 
0962       if (particle != nullptr) {
0963         Acts::Vector3 diffMom = recoMom - trueMom;
0964         // Accounting for the periodicity of phi. We overwrite the
0965         // previously computed value for better readability.
0966         diffMom[0] = Acts::detail::difference_periodic(recoMom(0), trueMom(0),
0967                                                        2 * std::numbers::pi);
0968         innerResPhi.push_back(diffMom[0]);
0969         innerResTheta.push_back(diffMom[1]);
0970         innerResQOverP.push_back(diffMom[2]);
0971 
0972         innerPullPhi.push_back(
0973             pull(diffMom[0], momCov(0, 0), "phi", false, logger()));
0974         innerPullTheta.push_back(
0975             pull(diffMom[1], momCov(1, 1), "theta", false, logger()));
0976         innerPullQOverP.push_back(
0977             pull(diffMom[2], momCov(2, 2), "q/p", false, logger()));
0978 
0979         const auto& recoUnitDir = paramsAtVtx->direction();
0980         double overlap = trueUnitDir.dot(recoUnitDir);
0981         innerMomOverlap.push_back(overlap);
0982       } else {
0983         innerResPhi.push_back(nan);
0984         innerResTheta.push_back(nan);
0985         innerResQOverP.push_back(nan);
0986         innerPullPhi.push_back(nan);
0987         innerPullTheta.push_back(nan);
0988         innerPullQOverP.push_back(nan);
0989         innerMomOverlap.push_back(nan);
0990       }
0991     } else {
0992       innerRecoPhi.push_back(nan);
0993       innerRecoTheta.push_back(nan);
0994       innerRecoQOverP.push_back(nan);
0995       innerResPhi.push_back(nan);
0996       innerResTheta.push_back(nan);
0997       innerResQOverP.push_back(nan);
0998       innerPullPhi.push_back(nan);
0999       innerPullTheta.push_back(nan);
1000       innerPullQOverP.push_back(nan);
1001       innerMomOverlap.push_back(nan);
1002     }
1003 
1004     // Save track parameters after the vertex fit
1005     const auto paramsAtVtxFitted = propagateToVtx(trk.fittedParams);
1006     if (paramsAtVtxFitted.has_value()) {
1007       Acts::Vector3 recoMomFitted =
1008           paramsAtVtxFitted->parameters().segment(Acts::eBoundPhi, 3);
1009       const Acts::SquareMatrix3& momCovFitted =
1010           paramsAtVtxFitted->covariance()->block<3, 3>(Acts::eBoundPhi,
1011                                                        Acts::eBoundPhi);
1012       innerRecoPhiFitted.push_back(recoMomFitted[0]);
1013       innerRecoThetaFitted.push_back(recoMomFitted[1]);
1014       innerRecoQOverPFitted.push_back(recoMomFitted[2]);
1015 
1016       if (particle != nullptr) {
1017         Acts::Vector3 diffMomFitted = recoMomFitted - trueMom;
1018         // Accounting for the periodicity of phi. We overwrite the
1019         // previously computed value for better readability.
1020         diffMomFitted[0] = Acts::detail::difference_periodic(
1021             recoMomFitted(0), trueMom(0), 2 * std::numbers::pi);
1022         innerResPhiFitted.push_back(diffMomFitted[0]);
1023         innerResThetaFitted.push_back(diffMomFitted[1]);
1024         innerResQOverPFitted.push_back(diffMomFitted[2]);
1025 
1026         innerPullPhiFitted.push_back(
1027             pull(diffMomFitted[0], momCovFitted(0, 0), "phi", true, logger()));
1028         innerPullThetaFitted.push_back(pull(
1029             diffMomFitted[1], momCovFitted(1, 1), "theta", true, logger()));
1030         innerPullQOverPFitted.push_back(
1031             pull(diffMomFitted[2], momCovFitted(2, 2), "q/p", true, logger()));
1032 
1033         const auto& recoUnitDirFitted = paramsAtVtxFitted->direction();
1034         double overlapFitted = trueUnitDir.dot(recoUnitDirFitted);
1035         innerMomOverlapFitted.push_back(overlapFitted);
1036       } else {
1037         innerResPhiFitted.push_back(nan);
1038         innerResThetaFitted.push_back(nan);
1039         innerResQOverPFitted.push_back(nan);
1040         innerPullPhiFitted.push_back(nan);
1041         innerPullThetaFitted.push_back(nan);
1042         innerPullQOverPFitted.push_back(nan);
1043         innerMomOverlapFitted.push_back(nan);
1044       }
1045     } else {
1046       innerRecoPhiFitted.push_back(nan);
1047       innerRecoThetaFitted.push_back(nan);
1048       innerRecoQOverPFitted.push_back(nan);
1049       innerResPhiFitted.push_back(nan);
1050       innerResThetaFitted.push_back(nan);
1051       innerResQOverPFitted.push_back(nan);
1052       innerPullPhiFitted.push_back(nan);
1053       innerPullThetaFitted.push_back(nan);
1054       innerPullQOverPFitted.push_back(nan);
1055       innerMomOverlapFitted.push_back(nan);
1056     }
1057   }
1058 }
1059 
1060 }  // namespace ActsExamples