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