Root/src/RootVertexNTupleWriter.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/RootVertexNTupleWriter.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 RootVertexNTupleWriter::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 RootVertexNTupleWriter::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 RootVertexNTupleWriter::RootVertexNTupleWriter(
0205     const RootVertexNTupleWriter::Config& config, Acts::Logging::Level level)
0206     : WriterT(config.inputVertices, "RootVertexNTupleWriter", 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_truthPhi", &m_truthPhi);
0330     m_outputTree->Branch("trk_truthTheta", &m_truthTheta);
0331     m_outputTree->Branch("trk_truthQOverP", &m_truthQOverP);
0332
0333     m_outputTree->Branch("trk_resPhi", &m_resPhi);
0334     m_outputTree->Branch("trk_resTheta", &m_resTheta);
0335     m_outputTree->Branch("trk_resQOverP", &m_resQOverP);
0336     m_outputTree->Branch("trk_momOverlap", &m_momOverlap);
0337
0338     m_outputTree->Branch("trk_resPhiFitted", &m_resPhiFitted);
0339     m_outputTree->Branch("trk_resThetaFitted", &m_resThetaFitted);
0340     m_outputTree->Branch("trk_resQOverPFitted", &m_resQOverPFitted);
0341     m_outputTree->Branch("trk_momOverlapFitted", &m_momOverlapFitted);
0342
0343     m_outputTree->Branch("trk_pullPhi", &m_pullPhi);
0344     m_outputTree->Branch("trk_pullTheta", &m_pullTheta);
0345     m_outputTree->Branch("trk_pullQOverP", &m_pullQOverP);
0346
0347     m_outputTree->Branch("trk_pullPhiFitted", &m_pullPhiFitted);
0348     m_outputTree->Branch("trk_pullThetaFitted", &m_pullThetaFitted);
0349     m_outputTree->Branch("trk_pullQOverPFitted", &m_pullQOverPFitted);
0350   }
0351 }
0352
0353 RootVertexNTupleWriter::~RootVertexNTupleWriter() {
0354   if (m_outputFile != nullptr) {
0355     m_outputFile->Close();
0356   }
0357 }
0358
0359 ProcessCode RootVertexNTupleWriter::finalize() {
0360   m_outputFile->cd();
0361   m_outputTree->Write();
0362   m_outputFile->Close();
0363
0364   return ProcessCode::SUCCESS;
0365 }
0366
0367 ProcessCode RootVertexNTupleWriter::writeT(
0368     const AlgorithmContext& ctx, const std::vector<Acts::Vertex>& vertices) {
0369   // In case we do not have any tracks in the vertex, we create empty
0370   // collections
0371   const static TrackParticleMatching emptyTrackParticleMatching;
0372   const static Acts::ConstVectorTrackContainer emptyConstTrackContainer;
0373   const static Acts::ConstVectorMultiTrajectory emptyConstTrackStateContainer;
0374   const static ConstTrackContainer emptyTracks(
0375       std::make_shared<Acts::ConstVectorTrackContainer>(
0376           emptyConstTrackContainer),
0377       std::make_shared<Acts::ConstVectorMultiTrajectory>(
0378           emptyConstTrackStateContainer));
0379
0380   // Read truth vertex input collection
0381   const SimVertexContainer& truthVertices = m_inputTruthVertices(ctx);
0382   // Read truth particle input collection
0383   const SimParticleContainer& particles = m_inputParticles(ctx);
0384   const SimParticleContainer& selectedParticles = m_inputSelectedParticles(ctx);
0385   const TrackParticleMatching& trackParticleMatching =
0386       (m_inputTrackParticleMatching.isInitialized()
0387            ? m_inputTrackParticleMatching(ctx)
0388            : emptyTrackParticleMatching);
0389   const ConstTrackContainer& tracks =
0390       (m_inputTracks.isInitialized() ? m_inputTracks(ctx) : emptyTracks);
0391   SimParticleContainer recoParticles;
0392
0393   for (ConstTrackProxy track : tracks) {
0394     if (!track.hasReferenceSurface()) {
0395       ACTS_DEBUG("No reference surface on this track, index = "
0396                  << track.index() << " tip index = " << track.tipIndex());
0397       continue;
0398     }
0399
0400     if (const SimParticle* particle =
0401             findParticle(particles, trackParticleMatching, track, logger());
0402         particle != nullptr) {
0403       recoParticles.insert(*particle);
0404     }
0405   }
0406   if (tracks.size() == 0) {
0407     // if not using tracks, then all truth particles are associated with the
0408     // vertex
0409     recoParticles = particles;
0410   }
0411
0412   // Exclusive access to the tree while writing
0413   std::lock_guard<std::mutex> lock(m_writeMutex);
0414
0415   m_nRecoVtx = vertices.size();
0416   m_nCleanVtx = 0;
0417   m_nMergedVtx = 0;
0418   m_nSplitVtx = 0;
0419
0420   ACTS_DEBUG("Number of reco vertices in event: " << m_nRecoVtx);
0421
0422   // Get number of generated true primary vertices
0423   m_nTrueVtx = getNumberOfTruePriVertices(particles);
0424   // Get number of detector-accepted true primary vertices
0425   m_nVtxDetAcceptance = getNumberOfTruePriVertices(selectedParticles);
0426
0427   ACTS_DEBUG("Number of truth particles in event : " << particles.size());
0428   ACTS_DEBUG("Number of truth primary vertices : " << m_nTrueVtx);
0429   ACTS_DEBUG("Number of detector-accepted truth primary vertices : "
0430              << m_nVtxDetAcceptance);
0431
0432   // Get the event number
0433   m_eventNr = ctx.eventNumber;
0434
0435   // Get number of track-associated true primary vertices
0436   m_nVtxReconstructable = getNumberOfReconstructableVertices(recoParticles);
0437
0438   ACTS_DEBUG("Number of reconstructed tracks : " << tracks.size());
0439   ACTS_DEBUG("Number of reco track-associated truth particles in event : "
0440              << recoParticles.size());
0441   ACTS_DEBUG("Maximum number of reconstructible primary vertices : "
0442              << m_nVtxReconstructable);
0443
0444   struct ToTruthMatching {
0445     std::optional<SimVertexBarcode> vertexId;
0446     double totalTrackWeight{};
0447     double truthMajorityTrackWeights{};
0448     double matchFraction{};
0449
0450     RecoVertexClassification classification{RecoVertexClassification::Unknown};
0451   };
0452   struct ToRecoMatching {
0453     std::size_t recoIndex{};
0454
0455     double recoSumPt2{};
0456   };
0457
0458   std::vector<ToTruthMatching> recoToTruthMatching;
0459   std::map<SimVertexBarcode, ToRecoMatching> truthToRecoMatching;
0460
0461   // Do truth matching for each reconstructed vertex
0462   for (const auto& [vtxIndex, vtx] : Acts::enumerate(vertices)) {
0463     // Reconstructed tracks that contribute to the reconstructed vertex
0464     const std::vector<Acts::TrackAtVertex>& tracksAtVtx = vtx.tracks();
0465
0466     // Containers for storing truth particles and truth vertices that
0467     // contribute to the reconstructed vertex
0468     std::vector<std::pair<SimVertexBarcode, double>> contributingTruthVertices;
0469
0470     double totalTrackWeight = 0;
0471     if (!tracksAtVtx.empty()) {
0472       for (const Acts::TrackAtVertex& trk : tracksAtVtx) {
0473         if (trk.trackWeight < m_cfg.minTrkWeight) {
0474           continue;
0475         }
0476
0477         totalTrackWeight += trk.trackWeight;
0478
0479         std::optional<ConstTrackProxy> trackOpt = findTrack(tracks, trk);
0480         if (!trackOpt.has_value()) {
0481           ACTS_DEBUG("Track has no matching input track.");
0482           continue;
0483         }
0484         const ConstTrackProxy& inputTrk = *trackOpt;
0485         const SimParticle* particle =
0486             findParticle(particles, trackParticleMatching, inputTrk, logger());
0487         if (particle == nullptr) {
0488           ACTS_VERBOSE("Track has no matching truth particle.");
0489         } else {
0490           contributingTruthVertices.emplace_back(
0491               SimBarcode{particle->particleId()}.vertexId(), trk.trackWeight);
0492         }
0493       }
0494     } else {
0495       // If no tracks at the vertex, then use all truth particles
0496       for (auto& particle : recoParticles) {
0497         contributingTruthVertices.emplace_back(
0498             SimBarcode{particle.particleId()}.vertexId(), 1.);
0499       }
0500     }
0501
0502     // Find true vertex that contributes most to the reconstructed vertex
0503     std::map<SimVertexBarcode, std::pair<int, double>> fmap;
0504     for (const auto& [vtxId, weight] : contributingTruthVertices) {
0505       ++fmap[vtxId].first;
0506       fmap[vtxId].second += weight;
0507     }
0508     double truthMajorityVertexTrackWeights = 0;
0509     std::optional<SimVertexBarcode> truthMajorityVertexId = std::nullopt;
0510     for (const auto& [vtxId, counter] : fmap) {
0511       if (counter.second > truthMajorityVertexTrackWeights) {
0512         truthMajorityVertexId = vtxId;
0513         truthMajorityVertexTrackWeights = counter.second;
0514       }
0515     }
0516
0517     double sumPt2 = calcSumPt2(m_cfg, vtx);
0518
0519     double vertexMatchFraction =
0520         (totalTrackWeight > 0
0521              ? truthMajorityVertexTrackWeights / totalTrackWeight
0522              : 0.);
0523     RecoVertexClassification recoVertexClassification =
0524         RecoVertexClassification::Unknown;
0525
0526     if (vertexMatchFraction >= m_cfg.vertexMatchThreshold) {
0527       recoVertexClassification = RecoVertexClassification::Clean;
0528     } else {
0529       recoVertexClassification = RecoVertexClassification::Merged;
0530     }
0531
0532     recoToTruthMatching.push_back({truthMajorityVertexId, totalTrackWeight,
0533                                    truthMajorityVertexTrackWeights,
0534                                    vertexMatchFraction,
0535                                    recoVertexClassification});
0536
0537     auto& recoToTruth = recoToTruthMatching.back();
0538
0539     // We have to decide if this reco vertex is a split vertex.
0540     if (!truthMajorityVertexId.has_value()) {
0541       // No truth vertex matched to this reconstructed vertex
0542       ACTS_DEBUG("No truth vertex matched to this reconstructed vertex.");
0543       continue;
0544     }
0545
0546     if (auto it = truthToRecoMatching.find(truthMajorityVertexId.value());
0547         it != truthToRecoMatching.end()) {
0548       // This truth vertex is already matched to a reconstructed vertex so we
0549       // are dealing with a split vertex.
0550
0551       // We have to decide which of the two reconstructed vertices is the
0552       // split vertex.
0553       if (sumPt2 <= it->second.recoSumPt2) {
0554         // Since the sumPt2 is smaller we can simply call this a split vertex
0555
0556         recoToTruth.classification = RecoVertexClassification::Split;
0557
0558         // Keep the existing truth to reco matching
0559       } else {
0560         // The sumPt2 is larger, so we call the other vertex a split vertex.
0561
0562         auto& otherRecoToTruth = recoToTruthMatching.at(it->second.recoIndex);
0563         // Swap the classification
0564         recoToTruth.classification = otherRecoToTruth.classification;
0565         otherRecoToTruth.classification = RecoVertexClassification::Split;
0566
0567         // Overwrite the truth to reco matching
0568         it->second = {vtxIndex, sumPt2};
0569       }
0570     } else {
0571       truthToRecoMatching[truthMajorityVertexId.value()] = {vtxIndex, sumPt2};
0572     }
0573   }
0574
0575   // Loop over reconstructed vertices and see if they can be matched to a true
0576   // vertex.
0577   for (const auto& [vtxIndex, vtx] : Acts::enumerate(vertices)) {
0578     // Reconstructed tracks that contribute to the reconstructed vertex
0579     const auto& tracksAtVtx = vtx.tracks();
0580     // Input tracks matched to `tracksAtVtx`
0581     std::vector<std::uint32_t> trackIndices;
0582
0583     const auto& toTruthMatching = recoToTruthMatching[vtxIndex];
0584
0585     m_recoX.push_back(vtx.fullPosition()[Acts::CoordinateIndices::eX]);
0586     m_recoY.push_back(vtx.fullPosition()[Acts::CoordinateIndices::eY]);
0587     m_recoZ.push_back(vtx.fullPosition()[Acts::CoordinateIndices::eZ]);
0588     m_recoT.push_back(vtx.fullPosition()[Acts::CoordinateIndices::eTime]);
0589
0590     double varX = vtx.fullCovariance()(Acts::CoordinateIndices::eX,
0591                                        Acts::CoordinateIndices::eX);
0592     double varY = vtx.fullCovariance()(Acts::CoordinateIndices::eY,
0593                                        Acts::CoordinateIndices::eY);
0594     double varZ = vtx.fullCovariance()(Acts::CoordinateIndices::eZ,
0595                                        Acts::CoordinateIndices::eZ);
0596     double varTime = vtx.fullCovariance()(Acts::CoordinateIndices::eTime,
0597                                           Acts::CoordinateIndices::eTime);
0598
0599     m_covXX.push_back(varX);
0600     m_covYY.push_back(varY);
0601     m_covZZ.push_back(varZ);
0602     m_covTT.push_back(varTime);
0603     m_covXY.push_back(vtx.fullCovariance()(Acts::CoordinateIndices::eX,
0604                                            Acts::CoordinateIndices::eY));
0605     m_covXZ.push_back(vtx.fullCovariance()(Acts::CoordinateIndices::eX,
0606                                            Acts::CoordinateIndices::eZ));
0607     m_covXT.push_back(vtx.fullCovariance()(Acts::CoordinateIndices::eX,
0608                                            Acts::CoordinateIndices::eTime));
0609     m_covYZ.push_back(vtx.fullCovariance()(Acts::CoordinateIndices::eY,
0610                                            Acts::CoordinateIndices::eZ));
0611     m_covYT.push_back(vtx.fullCovariance()(Acts::CoordinateIndices::eY,
0612                                            Acts::CoordinateIndices::eTime));
0613     m_covZT.push_back(vtx.fullCovariance()(Acts::CoordinateIndices::eZ,
0614                                            Acts::CoordinateIndices::eTime));
0615
0616     double sumPt2 = calcSumPt2(m_cfg, vtx);
0617     m_sumPt2.push_back(sumPt2);
0618
0619     double recoVertexTrackWeights = 0;
0620     for (const Acts::TrackAtVertex& trk : tracksAtVtx) {
0621       if (trk.trackWeight < m_cfg.minTrkWeight) {
0622         continue;
0623       }
0624       recoVertexTrackWeights += trk.trackWeight;
0625     }
0626     m_recoVertexTrackWeights.push_back(recoVertexTrackWeights);
0627
0628     unsigned int nTracksOnRecoVertex = std::count_if(
0629         tracksAtVtx.begin(), tracksAtVtx.end(), [this](const auto& trkAtVtx) {
0630           return trkAtVtx.trackWeight > m_cfg.minTrkWeight;
0631         });
0632     m_nTracksOnRecoVertex.push_back(nTracksOnRecoVertex);
0633
0634     // Saving truth information for the reconstructed vertex
0635     bool truthInfoWritten = false;
0636     std::optional<Acts::Vector4> truthPos;
0637     if (toTruthMatching.vertexId.has_value()) {
0638       auto iTruthVertex = truthVertices.find(toTruthMatching.vertexId.value());
0639       if (iTruthVertex == truthVertices.end()) {
0640         ACTS_ERROR("Truth vertex not found for id: "
0641                    << toTruthMatching.vertexId.value());
0642         continue;
0643       }
0644       const SimVertex& truthVertex = *iTruthVertex;
0645
0646       m_vertexPrimary.push_back(truthVertex.vertexId().vertexPrimary());
0647       m_vertexSecondary.push_back(truthVertex.vertexId().vertexSecondary());
0648
0649       // Count number of reconstructible tracks on truth vertex
0650       int nTracksOnTruthVertex = 0;
0651       for (const auto& particle : selectedParticles) {
0652         if (static_cast<SimVertexBarcode>(particle.particleId().vertexId()) ==
0653             truthVertex.vertexId()) {
0654           ++nTracksOnTruthVertex;
0655         }
0656       }
0657       m_nTracksOnTruthVertex.push_back(nTracksOnTruthVertex);
0658
0659       double truthPrimaryVertexDensity =
0660           calculateTruthPrimaryVertexDensity(m_cfg, truthVertices, vtx);
0661       m_truthPrimaryVertexDensity.push_back(truthPrimaryVertexDensity);
0662
0663       double truthVertexTrackWeights =
0664           toTruthMatching.truthMajorityTrackWeights;
0665       m_truthVertexTrackWeights.push_back(truthVertexTrackWeights);
0666
0667       double truthVertexMatchRatio = toTruthMatching.matchFraction;
0668       m_truthVertexMatchRatio.push_back(truthVertexMatchRatio);
0669
0670       double recoVertexContamination = 1 - truthVertexMatchRatio;
0671       m_recoVertexContamination.push_back(recoVertexContamination);
0672
0673       RecoVertexClassification recoVertexClassification =
0674           toTruthMatching.classification;
0675       m_recoVertexClassification.push_back(
0676           static_cast<int>(recoVertexClassification));
0677
0678       if (recoVertexClassification == RecoVertexClassification::Clean) {
0679         ++m_nCleanVtx;
0680       } else if (recoVertexClassification == RecoVertexClassification::Merged) {
0681         ++m_nMergedVtx;
0682       } else if (recoVertexClassification == RecoVertexClassification::Split) {
0683         ++m_nSplitVtx;
0684       }
0685
0686       const Acts::Vector4& truePos = truthVertex.position4;
0687       truthPos = truePos;
0688       m_truthX.push_back(truePos[Acts::CoordinateIndices::eX]);
0689       m_truthY.push_back(truePos[Acts::CoordinateIndices::eY]);
0690       m_truthZ.push_back(truePos[Acts::CoordinateIndices::eZ]);
0691       m_truthT.push_back(truePos[Acts::CoordinateIndices::eTime]);
0692
0693       const Acts::Vector4 diffPos = vtx.fullPosition() - truePos;
0694       m_resX.push_back(diffPos[Acts::CoordinateIndices::eX]);
0695       m_resY.push_back(diffPos[Acts::CoordinateIndices::eY]);
0696       m_resZ.push_back(diffPos[Acts::CoordinateIndices::eZ]);
0697       m_resT.push_back(diffPos[Acts::CoordinateIndices::eTime]);
0698
0699       m_pullX.push_back(pull(diffPos[Acts::CoordinateIndices::eX], varX, "X",
0700                              true, logger()));
0701       m_pullY.push_back(pull(diffPos[Acts::CoordinateIndices::eY], varY, "Y",
0702                              true, logger()));
0703       m_pullZ.push_back(pull(diffPos[Acts::CoordinateIndices::eZ], varZ, "Z",
0704                              true, logger()));
0705       m_pullT.push_back(pull(diffPos[Acts::CoordinateIndices::eTime], varTime,
0706                              "T", true, logger()));
0707
0708       truthInfoWritten = true;
0709     }
0710     if (!truthInfoWritten) {
0711       m_vertexPrimary.push_back(-1);
0712       m_vertexSecondary.push_back(-1);
0713
0714       m_nTracksOnTruthVertex.push_back(-1);
0715
0716       m_truthPrimaryVertexDensity.push_back(nan);
0717
0718       m_truthVertexTrackWeights.push_back(nan);
0719       m_truthVertexMatchRatio.push_back(nan);
0720       m_recoVertexContamination.push_back(nan);
0721
0722       m_recoVertexClassification.push_back(
0723           static_cast<int>(RecoVertexClassification::Unknown));
0724
0725       m_truthX.push_back(nan);
0726       m_truthY.push_back(nan);
0727       m_truthZ.push_back(nan);
0728       m_truthT.push_back(nan);
0729
0730       m_resX.push_back(nan);
0731       m_resY.push_back(nan);
0732       m_resZ.push_back(nan);
0733       m_resT.push_back(nan);
0734
0735       m_pullX.push_back(nan);
0736       m_pullY.push_back(nan);
0737       m_pullZ.push_back(nan);
0738       m_pullT.push_back(nan);
0739     }
0740
0741     if (m_cfg.writeTrackInfo) {
0742       writeTrackInfo(ctx, particles, tracks, trackParticleMatching, truthPos,
0743                      tracksAtVtx);
0744     }
0745   }
0746
0747   // fill the variables
0748   m_outputTree->Fill();
0749
0750   m_nTracksOnRecoVertex.clear();
0751   m_recoVertexTrackWeights.clear();
0752   m_recoX.clear();
0753   m_recoY.clear();
0754   m_recoZ.clear();
0755   m_recoT.clear();
0756   m_covXX.clear();
0757   m_covYY.clear();
0758   m_covZZ.clear();
0759   m_covTT.clear();
0760   m_covXY.clear();
0761   m_covXZ.clear();
0762   m_covXT.clear();
0763   m_covYZ.clear();
0764   m_covYT.clear();
0765   m_covZT.clear();
0766   m_seedX.clear();
0767   m_seedY.clear();
0768   m_seedZ.clear();
0769   m_seedT.clear();
0770   m_vertexPrimary.clear();
0771   m_vertexSecondary.clear();
0772   m_nTracksOnTruthVertex.clear();
0773   m_truthPrimaryVertexDensity.clear();
0774   m_truthVertexTrackWeights.clear();
0775   m_truthVertexMatchRatio.clear();
0776   m_recoVertexContamination.clear();
0777   m_recoVertexClassification.clear();
0778   m_truthX.clear();
0779   m_truthY.clear();
0780   m_truthZ.clear();
0781   m_truthT.clear();
0782   m_resX.clear();
0783   m_resY.clear();
0784   m_resZ.clear();
0785   m_resT.clear();
0786   m_resSeedZ.clear();
0787   m_resSeedT.clear();
0788   m_pullX.clear();
0789   m_pullY.clear();
0790   m_pullZ.clear();
0791   m_pullT.clear();
0792   m_sumPt2.clear();
0793   m_trkWeight.clear();
0794   m_recoPhi.clear();
0795   m_recoTheta.clear();
0796   m_recoQOverP.clear();
0797   m_recoPhiFitted.clear();
0798   m_recoThetaFitted.clear();
0799   m_recoQOverPFitted.clear();
0800   m_trkParticleIdVertexPrimary.clear();
0801   m_trkParticleIdVertexSecondary.clear();
0802   m_trkParticleIdParticle.clear();
0803   m_trkParticleIdGeneration.clear();
0804   m_trkParticleIdSubParticle.clear();
0805   m_truthPhi.clear();
0806   m_truthTheta.clear();
0807   m_truthQOverP.clear();
0808   m_resPhi.clear();
0809   m_resTheta.clear();
0810   m_resQOverP.clear();
0811   m_momOverlap.clear();
0812   m_resPhiFitted.clear();
0813   m_resThetaFitted.clear();
0814   m_resQOverPFitted.clear();
0815   m_momOverlapFitted.clear();
0816   m_pullPhi.clear();
0817   m_pullTheta.clear();
0818   m_pullQOverP.clear();
0819   m_pullPhiFitted.clear();
0820   m_pullThetaFitted.clear();
0821   m_pullQOverPFitted.clear();
0822
0823   return ProcessCode::SUCCESS;
0824 }
0825
0826 void RootVertexNTupleWriter::writeTrackInfo(
0827     const AlgorithmContext& ctx, const SimParticleContainer& particles,
0828     const ConstTrackContainer& tracks,
0829     const TrackParticleMatching& trackParticleMatching,
0830     const std::optional<Acts::Vector4>& truthPos,
0831     const std::vector<Acts::TrackAtVertex>& tracksAtVtx) {
0832   // We compare the reconstructed momenta to the true momenta at the vertex. For
0833   // this, we propagate the reconstructed tracks to the PCA of the true vertex
0834   // position. Setting up propagator:
0835   Acts::SympyStepper stepper(m_cfg.bField);
0836   using Propagator = Acts::Propagator<Acts::SympyStepper>;
0837   auto propagator = std::make_shared<Propagator>(stepper);
0838
0839   // Saving the reconstructed/truth momenta. The reconstructed momenta
0840   // are taken at the PCA to the truth vertex position -> we need to
0841   // perform a propagation.
0842
0843   // Get references to inner vectors where all track variables corresponding
0844   // to the current vertex will be saved
0845   auto& innerTrkWeight = m_trkWeight.emplace_back();
0846
0847   auto& innerRecoPhi = m_recoPhi.emplace_back();
0848   auto& innerRecoTheta = m_recoTheta.emplace_back();
0849   auto& innerRecoQOverP = m_recoQOverP.emplace_back();
0850
0851   auto& innerRecoPhiFitted = m_recoPhiFitted.emplace_back();
0852   auto& innerRecoThetaFitted = m_recoThetaFitted.emplace_back();
0853   auto& innerRecoQOverPFitted = m_recoQOverPFitted.emplace_back();
0854
0855   auto& innerTrkParticleIdVertexPrimary =
0856       m_trkParticleIdVertexPrimary.emplace_back();
0857   auto& innerTrkParticleIdVertexSecondary =
0858       m_trkParticleIdVertexSecondary.emplace_back();
0859   auto& innerTrkParticleIdParticle = m_trkParticleIdParticle.emplace_back();
0860   auto& innerTrkParticleIdGeneration = m_trkParticleIdGeneration.emplace_back();
0861   auto& innerTrkParticleIdSubParticle =
0862       m_trkParticleIdSubParticle.emplace_back();
0863
0864   auto& innerTruthPhi = m_truthPhi.emplace_back();
0865   auto& innerTruthTheta = m_truthTheta.emplace_back();
0866   auto& innerTruthQOverP = m_truthQOverP.emplace_back();
0867
0868   auto& innerResPhi = m_resPhi.emplace_back();
0869   auto& innerResTheta = m_resTheta.emplace_back();
0870   auto& innerResQOverP = m_resQOverP.emplace_back();
0871
0872   auto& innerResPhiFitted = m_resPhiFitted.emplace_back();
0873   auto& innerResThetaFitted = m_resThetaFitted.emplace_back();
0874   auto& innerResQOverPFitted = m_resQOverPFitted.emplace_back();
0875
0876   auto& innerMomOverlap = m_momOverlap.emplace_back();
0877   auto& innerMomOverlapFitted = m_momOverlapFitted.emplace_back();
0878
0879   auto& innerPullPhi = m_pullPhi.emplace_back();
0880   auto& innerPullTheta = m_pullTheta.emplace_back();
0881   auto& innerPullQOverP = m_pullQOverP.emplace_back();
0882
0883   auto& innerPullPhiFitted = m_pullPhiFitted.emplace_back();
0884   auto& innerPullThetaFitted = m_pullThetaFitted.emplace_back();
0885   auto& innerPullQOverPFitted = m_pullQOverPFitted.emplace_back();
0886
0887   // Perigee at the true vertex position
0888   std::shared_ptr<Acts::PerigeeSurface> perigeeSurface;
0889   if (truthPos.has_value()) {
0890     perigeeSurface =
0891         Acts::Surface::makeShared<Acts::PerigeeSurface>(truthPos->head<3>());
0892   }
0893   // Lambda for propagating the tracks to the PCA
0894   auto propagateToVtx =
0895       [&](const auto& params) -> std::optional<Acts::BoundTrackParameters> {
0896     if (!perigeeSurface) {
0897       return std::nullopt;
0898     }
0899
0900     Acts::Intersection3D intersection =
0901         perigeeSurface
0902             ->intersect(ctx.geoContext, params.position(ctx.geoContext),
0903                         params.direction(), Acts::BoundaryTolerance::Infinite())
0904             .closest();
0905
0906     // Setting the geometry/magnetic field context for the event
0907     using PropagatorOptions = Propagator::Options<>;
0908     PropagatorOptions pOptions(ctx.geoContext, ctx.magFieldContext);
0909     pOptions.direction =
0910         Acts::Direction::fromScalarZeroAsPositive(intersection.pathLength());
0911
0912     auto result = propagator->propagate(params, *perigeeSurface, pOptions);
0913     if (!result.ok()) {
0914       ACTS_ERROR("Propagation to true vertex position failed.");
0915       return std::nullopt;
0916     }
0917     auto& paramsAtVtx = *result->endParameters;
0918     return std::make_optional(paramsAtVtx);
0919   };
0920
0921   for (const Acts::TrackAtVertex& trk : tracksAtVtx) {
0922     if (trk.trackWeight < m_cfg.minTrkWeight) {
0923       continue;
0924     }
0925
0926     innerTrkWeight.push_back(trk.trackWeight);
0927
0928     Acts::Vector3 trueUnitDir = Acts::Vector3::Zero();
0929     Acts::Vector3 trueMom = Acts::Vector3::Zero();
0930
0931     const SimParticle* particle = nullptr;
0932     std::optional<ConstTrackProxy> trackOpt = findTrack(tracks, trk);
0933     if (trackOpt.has_value()) {
0934       particle =
0935           findParticle(particles, trackParticleMatching, *trackOpt, logger());
0936     }
0937
0938     if (particle != nullptr) {
0939       const auto barcode = particle->particleId();
0940       innerTrkParticleIdVertexPrimary.push_back(barcode.vertexPrimary());
0941       innerTrkParticleIdVertexSecondary.push_back(barcode.vertexSecondary());
0942       innerTrkParticleIdParticle.push_back(barcode.particle());
0943       innerTrkParticleIdGeneration.push_back(barcode.generation());
0944       innerTrkParticleIdSubParticle.push_back(barcode.subParticle());
0945
0946       trueUnitDir = particle->direction();
0947       trueMom.head<2>() = Acts::makePhiThetaFromDirection(trueUnitDir);
0948       trueMom[2] = particle->qOverP();
0949
0950       innerTruthPhi.push_back(trueMom[0]);
0951       innerTruthTheta.push_back(trueMom[1]);
0952       innerTruthQOverP.push_back(trueMom[2]);
0953     } else {
0954       ACTS_VERBOSE("Track has no matching truth particle.");
0955
0956       innerTrkParticleIdVertexPrimary.push_back(0);
0957       innerTrkParticleIdVertexSecondary.push_back(0);
0958       innerTrkParticleIdParticle.push_back(0);
0959       innerTrkParticleIdGeneration.push_back(0);
0960       innerTrkParticleIdSubParticle.push_back(0);
0961
0962       innerTruthPhi.push_back(nan);
0963       innerTruthTheta.push_back(nan);
0964       innerTruthQOverP.push_back(nan);
0965     }
0966
0967     // Save track parameters before the vertex fit
0968     const auto paramsAtVtx =
0969         propagateToVtx(*(trk.originalParams.as<Acts::BoundTrackParameters>()));
0970     if (paramsAtVtx.has_value()) {
0971       Acts::Vector3 recoMom =
0972           paramsAtVtx->parameters().segment(Acts::eBoundPhi, 3);
0973       const Acts::SquareMatrix3& momCov =
0974           paramsAtVtx->covariance()->template block<3, 3>(Acts::eBoundPhi,
0975                                                           Acts::eBoundPhi);
0976       innerRecoPhi.push_back(recoMom[0]);
0977       innerRecoTheta.push_back(recoMom[1]);
0978       innerRecoQOverP.push_back(recoMom[2]);
0979
0980       if (particle != nullptr) {
0981         Acts::Vector3 diffMom = recoMom - trueMom;
0982         // Accounting for the periodicity of phi. We overwrite the
0983         // previously computed value for better readability.
0984         diffMom[0] = Acts::detail::difference_periodic(recoMom(0), trueMom(0),
0985                                                        2 * std::numbers::pi);
0986         innerResPhi.push_back(diffMom[0]);
0987         innerResTheta.push_back(diffMom[1]);
0988         innerResQOverP.push_back(diffMom[2]);
0989
0990         innerPullPhi.push_back(
0991             pull(diffMom[0], momCov(0, 0), "phi", false, logger()));
0992         innerPullTheta.push_back(
0993             pull(diffMom[1], momCov(1, 1), "theta", false, logger()));
0994         innerPullQOverP.push_back(
0995             pull(diffMom[2], momCov(2, 2), "q/p", false, logger()));
0996
0997         const auto& recoUnitDir = paramsAtVtx->direction();
0998         double overlap = trueUnitDir.dot(recoUnitDir);
0999         innerMomOverlap.push_back(overlap);
1000       } else {
1001         innerResPhi.push_back(nan);
1002         innerResTheta.push_back(nan);
1003         innerResQOverP.push_back(nan);
1004         innerPullPhi.push_back(nan);
1005         innerPullTheta.push_back(nan);
1006         innerPullQOverP.push_back(nan);
1007         innerMomOverlap.push_back(nan);
1008       }
1009     } else {
1010       innerRecoPhi.push_back(nan);
1011       innerRecoTheta.push_back(nan);
1012       innerRecoQOverP.push_back(nan);
1013       innerResPhi.push_back(nan);
1014       innerResTheta.push_back(nan);
1015       innerResQOverP.push_back(nan);
1016       innerPullPhi.push_back(nan);
1017       innerPullTheta.push_back(nan);
1018       innerPullQOverP.push_back(nan);
1019       innerMomOverlap.push_back(nan);
1020     }
1021
1022     // Save track parameters after the vertex fit
1023     const auto paramsAtVtxFitted = propagateToVtx(trk.fittedParams);
1024     if (paramsAtVtxFitted.has_value()) {
1025       Acts::Vector3 recoMomFitted =
1026           paramsAtVtxFitted->parameters().segment(Acts::eBoundPhi, 3);
1027       const Acts::SquareMatrix3& momCovFitted =
1028           paramsAtVtxFitted->covariance()->block<3, 3>(Acts::eBoundPhi,
1029                                                        Acts::eBoundPhi);
1030       innerRecoPhiFitted.push_back(recoMomFitted[0]);
1031       innerRecoThetaFitted.push_back(recoMomFitted[1]);
1032       innerRecoQOverPFitted.push_back(recoMomFitted[2]);
1033
1034       if (particle != nullptr) {
1035         Acts::Vector3 diffMomFitted = recoMomFitted - trueMom;
1036         // Accounting for the periodicity of phi. We overwrite the
1037         // previously computed value for better readability.
1038         diffMomFitted[0] = Acts::detail::difference_periodic(
1039             recoMomFitted(0), trueMom(0), 2 * std::numbers::pi);
1040         innerResPhiFitted.push_back(diffMomFitted[0]);
1041         innerResThetaFitted.push_back(diffMomFitted[1]);
1042         innerResQOverPFitted.push_back(diffMomFitted[2]);
1043
1044         innerPullPhiFitted.push_back(
1045             pull(diffMomFitted[0], momCovFitted(0, 0), "phi", true, logger()));
1046         innerPullThetaFitted.push_back(pull(
1047             diffMomFitted[1], momCovFitted(1, 1), "theta", true, logger()));
1048         innerPullQOverPFitted.push_back(
1049             pull(diffMomFitted[2], momCovFitted(2, 2), "q/p", true, logger()));
1050
1051         const auto& recoUnitDirFitted = paramsAtVtxFitted->direction();
1052         double overlapFitted = trueUnitDir.dot(recoUnitDirFitted);
1053         innerMomOverlapFitted.push_back(overlapFitted);
1054       } else {
1055         innerResPhiFitted.push_back(nan);
1056         innerResThetaFitted.push_back(nan);
1057         innerResQOverPFitted.push_back(nan);
1058         innerPullPhiFitted.push_back(nan);
1059         innerPullThetaFitted.push_back(nan);
1060         innerPullQOverPFitted.push_back(nan);
1061         innerMomOverlapFitted.push_back(nan);
1062       }
1063     } else {
1064       innerRecoPhiFitted.push_back(nan);
1065       innerRecoThetaFitted.push_back(nan);
1066       innerRecoQOverPFitted.push_back(nan);
1067       innerResPhiFitted.push_back(nan);
1068       innerResThetaFitted.push_back(nan);
1069       innerResQOverPFitted.push_back(nan);
1070       innerPullPhiFitted.push_back(nan);
1071       innerPullThetaFitted.push_back(nan);
1072       innerPullQOverPFitted.push_back(nan);
1073       innerMomOverlapFitted.push_back(nan);
1074     }
1075   }
1076 }
1077
1078 }  // namespace ActsExamples