src/Validation/TrackParameterPerformanceCollector.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/Validation/TrackParameterPerformanceCollector.hpp"
0010
0011 #include "Acts/Surfaces/Surface.hpp"
0012 #include "Acts/Utilities/Logger.hpp"
0013 #include "Acts/Utilities/VectorHelpers.hpp"
0014 #include "ActsExamples/EventData/IndexSourceLink.hpp"
0015
0016 #include <stdexcept>
0017 #include <utility>
0018
0019 namespace ActsExamples {
0020
0021 TrackParameterPerformanceCollector::TrackParameterPerformanceCollector(
0022     Config cfg, std::unique_ptr<const Acts::Logger> logger)
0023     : m_cfg(std::move(cfg)),
0024       m_logger(std::move(logger)),
0025       m_resPlotTool(m_cfg.resPlotToolConfig, m_logger->level()),
0026       m_effPlotTool(m_cfg.effPlotToolConfig, m_logger->level()),
0027       m_trackSummaryPlotTool(m_cfg.trackSummaryPlotToolConfig,
0028                              m_logger->level()) {
0029   if (m_cfg.parameterSource == TrackParameterSource::Track &&
0030       (m_cfg.parameterType.has_value() || !m_cfg.geometrySelection.empty())) {
0031     throw std::invalid_argument(
0032         "Parameter type and geometry selection only apply to the TrackState "
0033         "parameter source");
0034   }
0035
0036   std::vector<Acts::GeometryHierarchyMap<unsigned int>::InputElement> elements;
0037   elements.reserve(m_cfg.geometrySelection.size());
0038   for (const Acts::GeometryIdentifier& geoId : m_cfg.geometrySelection) {
0039     elements.emplace_back(geoId, 0u);
0040   }
0041   m_geometrySelection =
0042       Acts::GeometryHierarchyMap<unsigned int>(std::move(elements));
0043 }
0044
0045 void TrackParameterPerformanceCollector::fill(
0046     const Acts::GeometryContext& geoContext, const ConstTrackContainer& tracks,
0047     const SimParticleContainer& particles,
0048     const TrackParticleMatching& trackParticleMatching,
0049     const SimHitContainer* simHits,
0050     const MeasurementSimHitsMap* measurementSimHitsMap) {
0051   // with the track-state source the comparison happens per measurement state
0052   // on the surface that state sits on, so the track itself needs no reference
0053   // surface, but the simulated hits behind the measurements are required
0054   const bool fromTrackStates =
0055       m_cfg.parameterSource == TrackParameterSource::TrackState;
0056   if (fromTrackStates &&
0057       (simHits == nullptr || measurementSimHitsMap == nullptr)) {
0058     throw std::invalid_argument(
0059         "Missing simulated hits for the TrackState parameter source");
0060   }
0061
0062   // Truth particles with corresponding reconstructed tracks
0063   std::vector<SimBarcode> reconParticleIds;
0064   reconParticleIds.reserve(tracks.size());
0065
0066   // Loop over all tracks
0067   for (const auto& track : tracks) {
0068     ++m_stats.nTotalTracks;
0069
0070     // Select reco track with fitted parameters
0071     if (!fromTrackStates && !track.hasReferenceSurface()) {
0072       ACTS_DEBUG("No fitted track parameters for track " << track.index());
0073       continue;
0074     }
0075
0076     // Get the truth-matched particle
0077     auto imatched = trackParticleMatching.find(track.index());
0078     if (imatched == trackParticleMatching.end()) {
0079       ACTS_DEBUG("No truth particle associated with track " << track.index());
0080       continue;
0081     }
0082     const auto& particleMatch = imatched->second;
0083
0084     if (!particleMatch.particle.has_value()) {
0085       ACTS_DEBUG("No truth particle associated with track " << track.index());
0086       continue;
0087     }
0088
0089     // Get the barcode of the majority truth particle
0090     SimBarcode majorityParticleId = particleMatch.particle.value();
0091
0092     // Find the truth particle via the barcode
0093     auto ip = particles.find(majorityParticleId);
0094     if (ip == particles.end()) {
0095       ACTS_DEBUG("Majority particle not found for track " << track.index());
0096       continue;
0097     }
0098
0099     // Record this majority particle ID
0100     reconParticleIds.push_back(ip->particleId());
0101
0102     if (fromTrackStates) {
0103       fillTrackStates(geoContext, track, *ip, *simHits, *measurementSimHitsMap);
0104       continue;
0105     }
0106
0107     Acts::BoundTrackParameters fittedParameters =
0108         track.createParametersAtReference();
0109
0110     // Fill residual plots
0111     m_resPlotTool.fill(geoContext, ip->initialState(), fittedParameters);
0112
0113     // Fill track summary info
0114     m_trackSummaryPlotTool.fill(fittedParameters, track.nTrackStates(),
0115                                 track.nMeasurements(), track.nOutliers(),
0116                                 track.nHoles(), track.nSharedHits());
0117   }
0118
0119   // Fill the efficiency
0120   for (const auto& particle : particles) {
0121     ++m_stats.nTotalParticles;
0122
0123     bool isReconstructed = false;
0124     if (Acts::rangeContainsValue(reconParticleIds, particle.particleId())) {
0125       isReconstructed = true;
0126       ++m_stats.nTotalMatchedTracks;
0127       ++m_stats.nTotalMatchedParticles;
0128     }
0129
0130     double minDeltaR = -1;
0131     for (const auto& closeParticle : particles) {
0132       if (closeParticle.particleId() == particle.particleId()) {
0133         continue;
0134       }
0135       double distance = Acts::VectorHelpers::deltaR(particle.direction(),
0136                                                     closeParticle.direction());
0137       if (minDeltaR == -1 || distance < minDeltaR) {
0138         minDeltaR = distance;
0139       }
0140     }
0141
0142     m_effPlotTool.fill(geoContext, particle.initialState(), minDeltaR,
0143                        isReconstructed);
0144   }
0145 }
0146
0147 void TrackParameterPerformanceCollector::fillTrackStates(
0148     const Acts::GeometryContext& geoContext, const ConstTrackProxy& track,
0149     const SimParticle& particle, const SimHitContainer& simHits,
0150     const MeasurementSimHitsMap& measurementSimHitsMap) {
0151   for (const auto& state : track.trackStatesReversed()) {
0152     if (!state.typeFlags().isMeasurement() || state.typeFlags().isOutlier()) {
0153       continue;
0154     }
0155     if (!state.hasReferenceSurface()) {
0156       continue;
0157     }
0158     const Acts::Surface& surface = state.referenceSurface();
0159
0160     if (!m_geometrySelection.empty() &&
0161         m_geometrySelection.find(surface.geometryId()) ==
0162             m_geometrySelection.end()) {
0163       continue;
0164     }
0165
0166     const std::optional<Acts::BoundTrackParameters> reco =
0167         recoParametersOnSurface(state, m_cfg.parameterType,
0168                                 track.particleHypothesis());
0169     if (!reco.has_value()) {
0170       ++m_stats.nMissingStateParameters;
0171       continue;
0172     }
0173
0174     // the source link must outlive the pointer into it
0175     const Acts::SourceLink sourceLink = state.getUncalibratedSourceLink();
0176     const auto* indexSourceLink = sourceLink.getPtr<IndexSourceLink>();
0177     if (indexSourceLink == nullptr) {
0178       ++m_stats.nMissingStateTruth;
0179       continue;
0180     }
0181
0182     const std::optional<Acts::BoundTrackParameters> truth =
0183         truthParametersOnSurface(geoContext, surface, indexSourceLink->index(),
0184                                  particle, simHits, measurementSimHitsMap,
0185                                  logger());
0186     if (!truth.has_value()) {
0187       ++m_stats.nMissingStateTruth;
0188       continue;
0189     }
0190
0191     m_resPlotTool.fill(truth.value(), reco.value());
0192   }
0193 }
0194
0195 void TrackParameterPerformanceCollector::logSummary() const {
0196   ACTS_INFO("=== Track Parameter Performance Summary ===");
0197   ACTS_INFO("Total tracks: " << m_stats.nTotalTracks);
0198   ACTS_INFO("Total matched tracks: " << m_stats.nTotalMatchedTracks);
0199   ACTS_INFO("Total particles: " << m_stats.nTotalParticles);
0200   ACTS_INFO("Total matched particles: " << m_stats.nTotalMatchedParticles);
0201
0202   if (m_cfg.parameterSource == TrackParameterSource::TrackState) {
0203     // a state counts here when it does not carry the requested parameters at
0204     // all, which is not a failure per se: an input that stores its estimate on
0205     // a single state, e.g. seeding output, skips every other state of a track
0206     ACTS_INFO("Skipped states without the requested parameters: "
0207               << m_stats.nMissingStateParameters);
0208     ACTS_INFO(
0209         "Skipped states without truth hits: " << m_stats.nMissingStateTruth);
0210   }
0211
0212   if (m_stats.nTotalTracks > 0) {
0213     double efficiency =
0214         static_cast<double>(m_stats.nTotalMatchedTracks) / m_stats.nTotalTracks;
0215     ACTS_INFO("Track efficiency: " << efficiency * 100 << "%");
0216   }
0217 }
0218
0219 template <std::size_t Dim>
0220 void TrackParameterPerformanceCollector::addFittedProfiles(
0221     const std::map<std::string, Acts::Experimental::Histogram<Dim>>& histMap,
0222     const std::string& meanPrefix, const std::string& widthPrefix,
0223     std::vector<Acts::Experimental::Histogram<Dim - 1>>& out) const {
0224   for (const auto& [name, hist] : histMap) {
0225     // Extract the suffix from the histogram name (e.g., "_d0_vs_eta")
0226     const std::string& baseName = hist.name();
0227     const std::string suffix = baseName.substr(baseName.find('_'));
0228
0229     auto profiles = extractMeanWidthProfiles(
0230         m_cfg.fitFunction, hist, meanPrefix + suffix, widthPrefix + suffix,
0231         m_cfg.fitMinEntries, m_cfg.fitSigmaRange, m_cfg.fitIterations,
0232         logger());
0233     if (profiles.fitFailureFraction >=
0234         m_cfg.warningThresholdFitFailureFraction) {
0235       ACTS_WARNING("Fit failures for " << baseName << ": "
0236                                        << profiles.fitFailureFraction * 100
0237                                        << "%");
0238     }
0239
0240     out.push_back(std::move(profiles.mean));
0241     out.push_back(std::move(profiles.width));
0242   }
0243 }
0244
0245 TrackParameterPerformanceCollector::FittedProfiles
0246 TrackParameterPerformanceCollector::fitProfiles() const {
0247   FittedProfiles profiles;
0248
0249   if (!m_cfg.fitFunction) {
0250     ACTS_WARNING(
0251         "No fit function configured; skipping mean/width profile "
0252         "extraction");
0253     return profiles;
0254   }
0255
0256   addFittedProfiles<2>(m_resPlotTool.resVsEta(), "resmean", "reswidth",
0257                        profiles.profiles1);
0258   addFittedProfiles<2>(m_resPlotTool.resVsPt(), "resmean", "reswidth",
0259                        profiles.profiles1);
0260   addFittedProfiles<3>(m_resPlotTool.resVsEtaPhi(), "resmean", "reswidth",
0261                        profiles.profiles2);
0262   addFittedProfiles<3>(m_resPlotTool.resVsEtaPt(), "resmean", "reswidth",
0263                        profiles.profiles2);
0264
0265   addFittedProfiles<2>(m_resPlotTool.pullVsEta(), "pullmean", "pullwidth",
0266                        profiles.profiles1);
0267   addFittedProfiles<2>(m_resPlotTool.pullVsPt(), "pullmean", "pullwidth",
0268                        profiles.profiles1);
0269   addFittedProfiles<3>(m_resPlotTool.pullVsEtaPhi(), "pullmean", "pullwidth",
0270                        profiles.profiles2);
0271   addFittedProfiles<3>(m_resPlotTool.pullVsEtaPt(), "pullmean", "pullwidth",
0272                        profiles.profiles2);
0273
0274   return profiles;
0275 }
0276
0277 }  // namespace ActsExamples