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     SimVertexBarcode truthMajorityVertexId{0};
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 (auto it = truthToRecoMatching.find(truthMajorityVertexId);
0542         it != truthToRecoMatching.end()) {
0543       // This truth vertex is already matched to a reconstructed vertex so we
0544       // are dealing with a split vertex.
0545
0546       // We have to decide which of the two reconstructed vertices is the
0547       // split vertex.
0548       if (sumPt2 <= it->second.recoSumPt2) {
0549         // Since the sumPt2 is smaller we can simply call this a split vertex
0550
0551         recoToTruth.classification = RecoVertexClassification::Split;
0552
0553         // Keep the existing truth to reco matching
0554       } else {
0555         // The sumPt2 is larger, so we call the other vertex a split vertex.
0556
0557         auto& otherRecoToTruth = recoToTruthMatching.at(it->second.recoIndex);
0558         // Swap the classification
0559         recoToTruth.classification = otherRecoToTruth.classification;
0560         otherRecoToTruth.classification = RecoVertexClassification::Split;
0561
0562         // Overwrite the truth to reco matching
0563         it->second = {vtxIndex, sumPt2};
0564       }
0565     } else {
0566       truthToRecoMatching[truthMajorityVertexId] = {vtxIndex, sumPt2};
0567     }
0568   }
0569
0570   // Loop over reconstructed vertices and see if they can be matched to a true
0571   // vertex.
0572   for (const auto& [vtxIndex, vtx] : Acts::enumerate(vertices)) {
0573     // Reconstructed tracks that contribute to the reconstructed vertex
0574     const auto& tracksAtVtx = vtx.tracks();
0575     // Input tracks matched to `tracksAtVtx`
0576     std::vector<std::uint32_t> trackIndices;
0577
0578     const auto& toTruthMatching = recoToTruthMatching[vtxIndex];
0579
0580     m_recoX.push_back(vtx.fullPosition()[Acts::CoordinateIndices::eX]);
0581     m_recoY.push_back(vtx.fullPosition()[Acts::CoordinateIndices::eY]);
0582     m_recoZ.push_back(vtx.fullPosition()[Acts::CoordinateIndices::eZ]);
0583     m_recoT.push_back(vtx.fullPosition()[Acts::CoordinateIndices::eTime]);
0584
0585     double varX = vtx.fullCovariance()(Acts::CoordinateIndices::eX,
0586                                        Acts::CoordinateIndices::eX);
0587     double varY = vtx.fullCovariance()(Acts::CoordinateIndices::eY,
0588                                        Acts::CoordinateIndices::eY);
0589     double varZ = vtx.fullCovariance()(Acts::CoordinateIndices::eZ,
0590                                        Acts::CoordinateIndices::eZ);
0591     double varTime = vtx.fullCovariance()(Acts::CoordinateIndices::eTime,
0592                                           Acts::CoordinateIndices::eTime);
0593
0594     m_covXX.push_back(varX);
0595     m_covYY.push_back(varY);
0596     m_covZZ.push_back(varZ);
0597     m_covTT.push_back(varTime);
0598     m_covXY.push_back(vtx.fullCovariance()(Acts::CoordinateIndices::eX,
0599                                            Acts::CoordinateIndices::eY));
0600     m_covXZ.push_back(vtx.fullCovariance()(Acts::CoordinateIndices::eX,
0601                                            Acts::CoordinateIndices::eZ));
0602     m_covXT.push_back(vtx.fullCovariance()(Acts::CoordinateIndices::eX,
0603                                            Acts::CoordinateIndices::eTime));
0604     m_covYZ.push_back(vtx.fullCovariance()(Acts::CoordinateIndices::eY,
0605                                            Acts::CoordinateIndices::eZ));
0606     m_covYT.push_back(vtx.fullCovariance()(Acts::CoordinateIndices::eY,
0607                                            Acts::CoordinateIndices::eTime));
0608     m_covZT.push_back(vtx.fullCovariance()(Acts::CoordinateIndices::eZ,
0609                                            Acts::CoordinateIndices::eTime));
0610
0611     double sumPt2 = calcSumPt2(m_cfg, vtx);
0612     m_sumPt2.push_back(sumPt2);
0613
0614     double recoVertexTrackWeights = 0;
0615     for (const Acts::TrackAtVertex& trk : tracksAtVtx) {
0616       if (trk.trackWeight < m_cfg.minTrkWeight) {
0617         continue;
0618       }
0619       recoVertexTrackWeights += trk.trackWeight;
0620     }
0621     m_recoVertexTrackWeights.push_back(recoVertexTrackWeights);
0622
0623     unsigned int nTracksOnRecoVertex = std::count_if(
0624         tracksAtVtx.begin(), tracksAtVtx.end(), [this](const auto& trkAtVtx) {
0625           return trkAtVtx.trackWeight > m_cfg.minTrkWeight;
0626         });
0627     m_nTracksOnRecoVertex.push_back(nTracksOnRecoVertex);
0628
0629     // Saving truth information for the reconstructed vertex
0630     bool truthInfoWritten = false;
0631     std::optional<Acts::Vector4> truthPos;
0632     if (toTruthMatching.vertexId.has_value()) {
0633       auto iTruthVertex = truthVertices.find(toTruthMatching.vertexId.value());
0634       if (iTruthVertex == truthVertices.end()) {
0635         ACTS_ERROR("Truth vertex not found.");
0636         continue;
0637       }
0638       const SimVertex& truthVertex = *iTruthVertex;
0639
0640       m_vertexPrimary.push_back(truthVertex.vertexId().vertexPrimary());
0641       m_vertexSecondary.push_back(truthVertex.vertexId().vertexSecondary());
0642
0643       // Count number of reconstructible tracks on truth vertex
0644       int nTracksOnTruthVertex = 0;
0645       for (const auto& particle : selectedParticles) {
0646         if (static_cast<SimVertexBarcode>(particle.particleId().vertexId()) ==
0647             truthVertex.vertexId()) {
0648           ++nTracksOnTruthVertex;
0649         }
0650       }
0651       m_nTracksOnTruthVertex.push_back(nTracksOnTruthVertex);
0652
0653       double truthPrimaryVertexDensity =
0654           calculateTruthPrimaryVertexDensity(m_cfg, truthVertices, vtx);
0655       m_truthPrimaryVertexDensity.push_back(truthPrimaryVertexDensity);
0656
0657       double truthVertexTrackWeights =
0658           toTruthMatching.truthMajorityTrackWeights;
0659       m_truthVertexTrackWeights.push_back(truthVertexTrackWeights);
0660
0661       double truthVertexMatchRatio = toTruthMatching.matchFraction;
0662       m_truthVertexMatchRatio.push_back(truthVertexMatchRatio);
0663
0664       double recoVertexContamination = 1 - truthVertexMatchRatio;
0665       m_recoVertexContamination.push_back(recoVertexContamination);
0666
0667       RecoVertexClassification recoVertexClassification =
0668           toTruthMatching.classification;
0669       m_recoVertexClassification.push_back(
0670           static_cast<int>(recoVertexClassification));
0671
0672       if (recoVertexClassification == RecoVertexClassification::Clean) {
0673         ++m_nCleanVtx;
0674       } else if (recoVertexClassification == RecoVertexClassification::Merged) {
0675         ++m_nMergedVtx;
0676       } else if (recoVertexClassification == RecoVertexClassification::Split) {
0677         ++m_nSplitVtx;
0678       }
0679
0680       const Acts::Vector4& truePos = truthVertex.position4;
0681       truthPos = truePos;
0682       m_truthX.push_back(truePos[Acts::CoordinateIndices::eX]);
0683       m_truthY.push_back(truePos[Acts::CoordinateIndices::eY]);
0684       m_truthZ.push_back(truePos[Acts::CoordinateIndices::eZ]);
0685       m_truthT.push_back(truePos[Acts::CoordinateIndices::eTime]);
0686
0687       const Acts::Vector4 diffPos = vtx.fullPosition() - truePos;
0688       m_resX.push_back(diffPos[Acts::CoordinateIndices::eX]);
0689       m_resY.push_back(diffPos[Acts::CoordinateIndices::eY]);
0690       m_resZ.push_back(diffPos[Acts::CoordinateIndices::eZ]);
0691       m_resT.push_back(diffPos[Acts::CoordinateIndices::eTime]);
0692
0693       m_pullX.push_back(pull(diffPos[Acts::CoordinateIndices::eX], varX, "X",
0694                              true, logger()));
0695       m_pullY.push_back(pull(diffPos[Acts::CoordinateIndices::eY], varY, "Y",
0696                              true, logger()));
0697       m_pullZ.push_back(pull(diffPos[Acts::CoordinateIndices::eZ], varZ, "Z",
0698                              true, logger()));
0699       m_pullT.push_back(pull(diffPos[Acts::CoordinateIndices::eTime], varTime,
0700                              "T", true, logger()));
0701
0702       truthInfoWritten = true;
0703     }
0704     if (!truthInfoWritten) {
0705       m_vertexPrimary.push_back(-1);
0706       m_vertexSecondary.push_back(-1);
0707
0708       m_nTracksOnTruthVertex.push_back(-1);
0709
0710       m_truthPrimaryVertexDensity.push_back(nan);
0711
0712       m_truthVertexTrackWeights.push_back(nan);
0713       m_truthVertexMatchRatio.push_back(nan);
0714       m_recoVertexContamination.push_back(nan);
0715
0716       m_recoVertexClassification.push_back(
0717           static_cast<int>(RecoVertexClassification::Unknown));
0718
0719       m_truthX.push_back(nan);
0720       m_truthY.push_back(nan);
0721       m_truthZ.push_back(nan);
0722       m_truthT.push_back(nan);
0723
0724       m_resX.push_back(nan);
0725       m_resY.push_back(nan);
0726       m_resZ.push_back(nan);
0727       m_resT.push_back(nan);
0728
0729       m_pullX.push_back(nan);
0730       m_pullY.push_back(nan);
0731       m_pullZ.push_back(nan);
0732       m_pullT.push_back(nan);
0733     }
0734
0735     if (m_cfg.writeTrackInfo) {
0736       writeTrackInfo(ctx, particles, tracks, trackParticleMatching, truthPos,
0737                      tracksAtVtx);
0738     }
0739   }
0740
0741   // fill the variables
0742   m_outputTree->Fill();
0743
0744   m_nTracksOnRecoVertex.clear();
0745   m_recoVertexTrackWeights.clear();
0746   m_recoX.clear();
0747   m_recoY.clear();
0748   m_recoZ.clear();
0749   m_recoT.clear();
0750   m_covXX.clear();
0751   m_covYY.clear();
0752   m_covZZ.clear();
0753   m_covTT.clear();
0754   m_covXY.clear();
0755   m_covXZ.clear();
0756   m_covXT.clear();
0757   m_covYZ.clear();
0758   m_covYT.clear();
0759   m_covZT.clear();
0760   m_seedX.clear();
0761   m_seedY.clear();
0762   m_seedZ.clear();
0763   m_seedT.clear();
0764   m_vertexPrimary.clear();
0765   m_vertexSecondary.clear();
0766   m_nTracksOnTruthVertex.clear();
0767   m_truthPrimaryVertexDensity.clear();
0768   m_truthVertexTrackWeights.clear();
0769   m_truthVertexMatchRatio.clear();
0770   m_recoVertexContamination.clear();
0771   m_recoVertexClassification.clear();
0772   m_truthX.clear();
0773   m_truthY.clear();
0774   m_truthZ.clear();
0775   m_truthT.clear();
0776   m_resX.clear();
0777   m_resY.clear();
0778   m_resZ.clear();
0779   m_resT.clear();
0780   m_resSeedZ.clear();
0781   m_resSeedT.clear();
0782   m_pullX.clear();
0783   m_pullY.clear();
0784   m_pullZ.clear();
0785   m_pullT.clear();
0786   m_sumPt2.clear();
0787   m_trkWeight.clear();
0788   m_recoPhi.clear();
0789   m_recoTheta.clear();
0790   m_recoQOverP.clear();
0791   m_recoPhiFitted.clear();
0792   m_recoThetaFitted.clear();
0793   m_recoQOverPFitted.clear();
0794   m_trkParticleId.clear();
0795   m_truthPhi.clear();
0796   m_truthTheta.clear();
0797   m_truthQOverP.clear();
0798   m_resPhi.clear();
0799   m_resTheta.clear();
0800   m_resQOverP.clear();
0801   m_momOverlap.clear();
0802   m_resPhiFitted.clear();
0803   m_resThetaFitted.clear();
0804   m_resQOverPFitted.clear();
0805   m_momOverlapFitted.clear();
0806   m_pullPhi.clear();
0807   m_pullTheta.clear();
0808   m_pullQOverP.clear();
0809   m_pullPhiFitted.clear();
0810   m_pullThetaFitted.clear();
0811   m_pullQOverPFitted.clear();
0812
0813   return ProcessCode::SUCCESS;
0814 }
0815
0816 void VertexNTupleWriter::writeTrackInfo(
0817     const AlgorithmContext& ctx, const SimParticleContainer& particles,
0818     const ConstTrackContainer& tracks,
0819     const TrackParticleMatching& trackParticleMatching,
0820     const std::optional<Acts::Vector4>& truthPos,
0821     const std::vector<Acts::TrackAtVertex>& tracksAtVtx) {
0822   // We compare the reconstructed momenta to the true momenta at the vertex. For
0823   // this, we propagate the reconstructed tracks to the PCA of the true vertex
0824   // position. Setting up propagator:
0825   Acts::SympyStepper stepper(m_cfg.bField);
0826   using Propagator = Acts::Propagator<Acts::SympyStepper>;
0827   auto propagator = std::make_shared<Propagator>(stepper);
0828
0829   // Saving the reconstructed/truth momenta. The reconstructed momenta
0830   // are taken at the PCA to the truth vertex position -> we need to
0831   // perform a propagation.
0832
0833   // Get references to inner vectors where all track variables corresponding
0834   // to the current vertex will be saved
0835   auto& innerTrkWeight = m_trkWeight.emplace_back();
0836
0837   auto& innerRecoPhi = m_recoPhi.emplace_back();
0838   auto& innerRecoTheta = m_recoTheta.emplace_back();
0839   auto& innerRecoQOverP = m_recoQOverP.emplace_back();
0840
0841   auto& innerRecoPhiFitted = m_recoPhiFitted.emplace_back();
0842   auto& innerRecoThetaFitted = m_recoThetaFitted.emplace_back();
0843   auto& innerRecoQOverPFitted = m_recoQOverPFitted.emplace_back();
0844
0845   auto& innerTrkParticleId = m_trkParticleId.emplace_back();
0846
0847   auto& innerTruthPhi = m_truthPhi.emplace_back();
0848   auto& innerTruthTheta = m_truthTheta.emplace_back();
0849   auto& innerTruthQOverP = m_truthQOverP.emplace_back();
0850
0851   auto& innerResPhi = m_resPhi.emplace_back();
0852   auto& innerResTheta = m_resTheta.emplace_back();
0853   auto& innerResQOverP = m_resQOverP.emplace_back();
0854
0855   auto& innerResPhiFitted = m_resPhiFitted.emplace_back();
0856   auto& innerResThetaFitted = m_resThetaFitted.emplace_back();
0857   auto& innerResQOverPFitted = m_resQOverPFitted.emplace_back();
0858
0859   auto& innerMomOverlap = m_momOverlap.emplace_back();
0860   auto& innerMomOverlapFitted = m_momOverlapFitted.emplace_back();
0861
0862   auto& innerPullPhi = m_pullPhi.emplace_back();
0863   auto& innerPullTheta = m_pullTheta.emplace_back();
0864   auto& innerPullQOverP = m_pullQOverP.emplace_back();
0865
0866   auto& innerPullPhiFitted = m_pullPhiFitted.emplace_back();
0867   auto& innerPullThetaFitted = m_pullThetaFitted.emplace_back();
0868   auto& innerPullQOverPFitted = m_pullQOverPFitted.emplace_back();
0869
0870   // Perigee at the true vertex position
0871   std::shared_ptr<Acts::PerigeeSurface> perigeeSurface;
0872   if (truthPos.has_value()) {
0873     perigeeSurface =
0874         Acts::Surface::makeShared<Acts::PerigeeSurface>(truthPos->head<3>());
0875   }
0876   // Lambda for propagating the tracks to the PCA
0877   auto propagateToVtx =
0878       [&](const auto& params) -> std::optional<Acts::BoundTrackParameters> {
0879     if (!perigeeSurface) {
0880       return std::nullopt;
0881     }
0882
0883     auto intersection =
0884         perigeeSurface
0885             ->intersect(ctx.geoContext, params.position(ctx.geoContext),
0886                         params.direction(), Acts::BoundaryTolerance::Infinite())
0887             .closest();
0888
0889     // Setting the geometry/magnetic field context for the event
0890     using PropagatorOptions = Propagator::Options<>;
0891     PropagatorOptions pOptions(ctx.geoContext, ctx.magFieldContext);
0892     pOptions.direction =
0893         Acts::Direction::fromScalarZeroAsPositive(intersection.pathLength());
0894
0895     auto result = propagator->propagate(params, *perigeeSurface, pOptions);
0896     if (!result.ok()) {
0897       ACTS_ERROR("Propagation to true vertex position failed.");
0898       return std::nullopt;
0899     }
0900     auto& paramsAtVtx = *result->endParameters;
0901     return std::make_optional(paramsAtVtx);
0902   };
0903
0904   for (const Acts::TrackAtVertex& trk : tracksAtVtx) {
0905     if (trk.trackWeight < m_cfg.minTrkWeight) {
0906       continue;
0907     }
0908
0909     innerTrkWeight.push_back(trk.trackWeight);
0910
0911     Acts::Vector3 trueUnitDir = Acts::Vector3::Zero();
0912     Acts::Vector3 trueMom = Acts::Vector3::Zero();
0913
0914     const SimParticle* particle = nullptr;
0915     std::optional<ConstTrackProxy> trackOpt = findTrack(tracks, trk);
0916     if (trackOpt.has_value()) {
0917       particle =
0918           findParticle(particles, trackParticleMatching, *trackOpt, logger());
0919     }
0920
0921     if (particle != nullptr) {
0922       innerTrkParticleId.push_back(particle->particleId().value());
0923
0924       trueUnitDir = particle->direction();
0925       trueMom.head<2>() = Acts::makePhiThetaFromDirection(trueUnitDir);
0926       trueMom[2] = particle->qOverP();
0927
0928       innerTruthPhi.push_back(trueMom[0]);
0929       innerTruthTheta.push_back(trueMom[1]);
0930       innerTruthQOverP.push_back(trueMom[2]);
0931     } else {
0932       ACTS_VERBOSE("Track has no matching truth particle.");
0933
0934       innerTrkParticleId.push_back(-1);
0935
0936       innerTruthPhi.push_back(nan);
0937       innerTruthTheta.push_back(nan);
0938       innerTruthQOverP.push_back(nan);
0939     }
0940
0941     // Save track parameters before the vertex fit
0942     const auto paramsAtVtx =
0943         propagateToVtx(*(trk.originalParams.as<Acts::BoundTrackParameters>()));
0944     if (paramsAtVtx.has_value()) {
0945       Acts::Vector3 recoMom =
0946           paramsAtVtx->parameters().segment(Acts::eBoundPhi, 3);
0947       const Acts::SquareMatrix3& momCov =
0948           paramsAtVtx->covariance()->template block<3, 3>(Acts::eBoundPhi,
0949                                                           Acts::eBoundPhi);
0950       innerRecoPhi.push_back(recoMom[0]);
0951       innerRecoTheta.push_back(recoMom[1]);
0952       innerRecoQOverP.push_back(recoMom[2]);
0953
0954       if (particle != nullptr) {
0955         Acts::Vector3 diffMom = recoMom - trueMom;
0956         // Accounting for the periodicity of phi. We overwrite the
0957         // previously computed value for better readability.
0958         diffMom[0] = Acts::detail::difference_periodic(recoMom(0), trueMom(0),
0959                                                        2 * std::numbers::pi);
0960         innerResPhi.push_back(diffMom[0]);
0961         innerResTheta.push_back(diffMom[1]);
0962         innerResQOverP.push_back(diffMom[2]);
0963
0964         innerPullPhi.push_back(
0965             pull(diffMom[0], momCov(0, 0), "phi", false, logger()));
0966         innerPullTheta.push_back(
0967             pull(diffMom[1], momCov(1, 1), "theta", false, logger()));
0968         innerPullQOverP.push_back(
0969             pull(diffMom[2], momCov(2, 2), "q/p", false, logger()));
0970
0971         const auto& recoUnitDir = paramsAtVtx->direction();
0972         double overlap = trueUnitDir.dot(recoUnitDir);
0973         innerMomOverlap.push_back(overlap);
0974       } else {
0975         innerResPhi.push_back(nan);
0976         innerResTheta.push_back(nan);
0977         innerResQOverP.push_back(nan);
0978         innerPullPhi.push_back(nan);
0979         innerPullTheta.push_back(nan);
0980         innerPullQOverP.push_back(nan);
0981         innerMomOverlap.push_back(nan);
0982       }
0983     } else {
0984       innerRecoPhi.push_back(nan);
0985       innerRecoTheta.push_back(nan);
0986       innerRecoQOverP.push_back(nan);
0987       innerResPhi.push_back(nan);
0988       innerResTheta.push_back(nan);
0989       innerResQOverP.push_back(nan);
0990       innerPullPhi.push_back(nan);
0991       innerPullTheta.push_back(nan);
0992       innerPullQOverP.push_back(nan);
0993       innerMomOverlap.push_back(nan);
0994     }
0995
0996     // Save track parameters after the vertex fit
0997     const auto paramsAtVtxFitted = propagateToVtx(trk.fittedParams);
0998     if (paramsAtVtxFitted.has_value()) {
0999       Acts::Vector3 recoMomFitted =
1000           paramsAtVtxFitted->parameters().segment(Acts::eBoundPhi, 3);
1001       const Acts::SquareMatrix3& momCovFitted =
1002           paramsAtVtxFitted->covariance()->block<3, 3>(Acts::eBoundPhi,
1003                                                        Acts::eBoundPhi);
1004       innerRecoPhiFitted.push_back(recoMomFitted[0]);
1005       innerRecoThetaFitted.push_back(recoMomFitted[1]);
1006       innerRecoQOverPFitted.push_back(recoMomFitted[2]);
1007
1008       if (particle != nullptr) {
1009         Acts::Vector3 diffMomFitted = recoMomFitted - trueMom;
1010         // Accounting for the periodicity of phi. We overwrite the
1011         // previously computed value for better readability.
1012         diffMomFitted[0] = Acts::detail::difference_periodic(
1013             recoMomFitted(0), trueMom(0), 2 * std::numbers::pi);
1014         innerResPhiFitted.push_back(diffMomFitted[0]);
1015         innerResThetaFitted.push_back(diffMomFitted[1]);
1016         innerResQOverPFitted.push_back(diffMomFitted[2]);
1017
1018         innerPullPhiFitted.push_back(
1019             pull(diffMomFitted[0], momCovFitted(0, 0), "phi", true, logger()));
1020         innerPullThetaFitted.push_back(pull(
1021             diffMomFitted[1], momCovFitted(1, 1), "theta", true, logger()));
1022         innerPullQOverPFitted.push_back(
1023             pull(diffMomFitted[2], momCovFitted(2, 2), "q/p", true, logger()));
1024
1025         const auto& recoUnitDirFitted = paramsAtVtxFitted->direction();
1026         double overlapFitted = trueUnitDir.dot(recoUnitDirFitted);
1027         innerMomOverlapFitted.push_back(overlapFitted);
1028       } else {
1029         innerResPhiFitted.push_back(nan);
1030         innerResThetaFitted.push_back(nan);
1031         innerResQOverPFitted.push_back(nan);
1032         innerPullPhiFitted.push_back(nan);
1033         innerPullThetaFitted.push_back(nan);
1034         innerPullQOverPFitted.push_back(nan);
1035         innerMomOverlapFitted.push_back(nan);
1036       }
1037     } else {
1038       innerRecoPhiFitted.push_back(nan);
1039       innerRecoThetaFitted.push_back(nan);
1040       innerRecoQOverPFitted.push_back(nan);
1041       innerResPhiFitted.push_back(nan);
1042       innerResThetaFitted.push_back(nan);
1043       innerResQOverPFitted.push_back(nan);
1044       innerPullPhiFitted.push_back(nan);
1045       innerPullThetaFitted.push_back(nan);
1046       innerPullQOverPFitted.push_back(nan);
1047       innerMomOverlapFitted.push_back(nan);
1048     }
1049   }
1050 }
1051
1052 }  // namespace ActsExamples