diff --git a/Examples/Framework/CMakeLists.txt b/Examples/Framework/CMakeLists.txt index 70342a78b5f..59352e29393 100644 --- a/Examples/Framework/CMakeLists.txt +++ b/Examples/Framework/CMakeLists.txt @@ -26,6 +26,7 @@ acts_add_library( src/Validation/EffPlotTool.cpp src/Validation/FakePlotTool.cpp src/Validation/HistogramFit.cpp + src/Validation/ParametersOnSurface.cpp src/Validation/ResPlotTool.cpp src/Validation/TrackClassification.cpp src/Validation/PatternRecognitionPerformanceCollector.cpp diff --git a/Examples/Framework/include/ActsExamples/Validation/ParametersOnSurface.hpp b/Examples/Framework/include/ActsExamples/Validation/ParametersOnSurface.hpp new file mode 100644 index 00000000000..edd63ae2db9 --- /dev/null +++ b/Examples/Framework/include/ActsExamples/Validation/ParametersOnSurface.hpp @@ -0,0 +1,80 @@ +// This file is part of the ACTS project. +// +// Copyright (C) 2016 CERN for the benefit of the ACTS project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at https://mozilla.org/MPL/2.0/. + +#pragma once + +#include "Acts/EventData/BoundTrackParameters.hpp" +#include "Acts/Geometry/GeometryContext.hpp" +#include "Acts/Utilities/Logger.hpp" +#include "ActsExamples/EventData/Index.hpp" +#include "ActsExamples/EventData/SimHit.hpp" +#include "ActsExamples/EventData/SimParticle.hpp" +#include "ActsExamples/EventData/Track.hpp" +#include "ActsExamples/EventData/TruthMatching.hpp" + +#include + +namespace Acts { +class Surface; +} + +namespace ActsExamples { + +/// The track-state parameters to extract. +enum class TrackParameterType { + /// Parameters before the measurement update + Predicted, + /// Parameters after the measurement update + Filtered, + /// Smoothed parameters + Smoothed, + /// Smoothed parameters with the state's own measurement removed + Unbiased, +}; + +/// Compute the truth bound track parameters on a surface from the simulated +/// hits of a measurement. +/// +/// Position, direction, and time are averaged over the simulated hits, the +/// momentum is taken from the first one, and charge and hypothesis from the +/// given particle. +/// +/// @note The hits are averaged regardless of which particle produced them, so +/// a merged cluster gives a mixture. Same as `RootTrackParameterWriter`. +/// +/// @param gctx The geometry context +/// @param surface The surface to express the truth parameters on +/// @param measurementIndex The index of the measurement on the surface +/// @param particle The truth particle +/// @param simHits The simulated hits container +/// @param measurementSimHitsMap Map from measurement index to simulated hits +/// @param logger A logger for messages +/// @return The parameters without covariance, or nullopt without truth hits +std::optional truthParametersOnSurface( + const Acts::GeometryContext& gctx, const Acts::Surface& surface, + Index measurementIndex, const SimParticle& particle, + const SimHitContainer& simHits, + const MeasurementSimHitsMap& measurementSimHitsMap, + const Acts::Logger& logger); + +/// Extract the reconstructed bound track parameters on the reference surface +/// of a track state. +/// +/// @param state The track state to extract the parameters from +/// @param parameterType Which parameters; if not set, the best available ones +/// (smoothed, filtered, or predicted). `Unbiased` is never picked +/// implicitly and has to be requested +/// @param hypothesis The particle hypothesis for the parameters +/// @return The parameters with covariance, or nullopt without a reference +/// surface or matching parameters +std::optional recoParametersOnSurface( + const ConstTrackStateProxy& state, + std::optional parameterType, + const Acts::ParticleHypothesis& hypothesis); + +} // namespace ActsExamples diff --git a/Examples/Framework/include/ActsExamples/Validation/ResPlotTool.hpp b/Examples/Framework/include/ActsExamples/Validation/ResPlotTool.hpp index d46323a5f29..65c5196802f 100644 --- a/Examples/Framework/include/ActsExamples/Validation/ResPlotTool.hpp +++ b/Examples/Framework/include/ActsExamples/Validation/ResPlotTool.hpp @@ -8,6 +8,7 @@ #pragma once +#include "Acts/EventData/BoundTrackParameters.hpp" #include "Acts/Geometry/GeometryContext.hpp" #include "Acts/Utilities/Histogram.hpp" #include "Acts/Utilities/Logger.hpp" @@ -22,8 +23,7 @@ namespace ActsExamples { /// Tools to make hists to show residual, i.e. smoothed_parameter - /// truth_parameter, and pull, i.e. (smoothed_parameter - -/// truth_parameter)/smoothed_paramter_error, of track parameters at perigee -/// surface +/// truth_parameter)/smoothed_paramter_error, of bound track parameters class ResPlotTool { public: using AxisVariant = Acts::Experimental::AxisVariant; @@ -71,6 +71,13 @@ class ResPlotTool { const SimParticleState& truthParticle, const Acts::BoundTrackParameters& fittedParamters); + /// Fill from truth parameters on the same surface as the fitted ones. + /// + /// @param truthParameters the truth bound parameters + /// @param fittedParameters the fitted parameters + void fill(const Acts::BoundTrackParameters& truthParameters, + const Acts::BoundTrackParameters& fittedParameters); + const std::map& res() const { return m_res; } const std::map& resVsEta() const { return m_resVsEta; @@ -97,6 +104,16 @@ class ResPlotTool { } private: + /// Truth quantities the histograms are binned in. Each `fill` overload + /// derives them from its own truth source to avoid lossy conversions. + struct TruthBinning { + double eta; + double phi; + double pt; + double charge; + double absCharge; + }; + Config m_cfg; std::unique_ptr m_logger; @@ -123,6 +140,10 @@ class ResPlotTool { /// Pull vs eta-pT scatter plot std::map m_pullVsEtaPt; + void fill(const Acts::BoundVector& truthVector, + const TruthBinning& truthBinning, + const Acts::BoundTrackParameters& fittedParameters); + void fillResidual(const std::string& paramName, double residual, double truthEta, double truthPhi, double truthPt); void fillPull(const std::string& paramName, double pull, double truthEta, diff --git a/Examples/Framework/include/ActsExamples/Validation/TrackParameterPerformanceCollector.hpp b/Examples/Framework/include/ActsExamples/Validation/TrackParameterPerformanceCollector.hpp index a46338f506c..4175a685d37 100644 --- a/Examples/Framework/include/ActsExamples/Validation/TrackParameterPerformanceCollector.hpp +++ b/Examples/Framework/include/ActsExamples/Validation/TrackParameterPerformanceCollector.hpp @@ -9,30 +9,50 @@ #pragma once #include "Acts/Geometry/GeometryContext.hpp" +#include "Acts/Geometry/GeometryHierarchyMap.hpp" +#include "Acts/Geometry/GeometryIdentifier.hpp" #include "Acts/Utilities/Histogram.hpp" #include "Acts/Utilities/Logger.hpp" +#include "ActsExamples/EventData/SimHit.hpp" #include "ActsExamples/EventData/SimParticle.hpp" #include "ActsExamples/EventData/Track.hpp" #include "ActsExamples/EventData/TruthMatching.hpp" #include "ActsExamples/Validation/EffPlotTool.hpp" #include "ActsExamples/Validation/HistogramFit.hpp" +#include "ActsExamples/Validation/ParametersOnSurface.hpp" #include "ActsExamples/Validation/ResPlotTool.hpp" #include "ActsExamples/Validation/TrackSummaryPlotTool.hpp" #include #include +#include #include #include namespace ActsExamples { -/// Collects performance histograms of the track parameters at the track -/// reference surface, without any file I/O. +/// Where to take the reconstructed parameters from. +enum class TrackParameterSource { + /// The track parameters at the track reference surface, i.e. the fitter + /// output as delivered. + Track, + /// The parameters of the individual track states, on the surface they sit + /// on. Compared against the simulated hits of the state's measurement. + TrackState, +}; + +/// Collects performance histograms of the track parameters, without any file +/// I/O. /// /// Collects residual/pull histograms, efficiency plots, and track summary /// information for track fitting performance evaluation. The Gaussian fit /// backend is supplied by the caller via @c Config::fitFunction. /// +/// With `parameterSource = Track` the track parameters at the track reference +/// surface are compared to the truth particle. With `TrackState` every +/// selected measurement state is compared to the truth on its own surface, +/// which is what makes per-sensor estimates, e.g. from a seed, measurable. +/// /// @note The caller must ensure exclusive access (e.g. hold a mutex) when /// calling fill(). This class applies no locking of its own. class TrackParameterPerformanceCollector { @@ -42,6 +62,15 @@ class TrackParameterPerformanceCollector { EffPlotTool::Config effPlotToolConfig; TrackSummaryPlotTool::Config trackSummaryPlotToolConfig; + /// Where to take the reconstructed parameters from. + TrackParameterSource parameterSource = TrackParameterSource::Track; + /// Which track-state parameters to use. If not set, the best available + /// ones (smoothed, filtered, or predicted). `TrackState` source only. + std::optional parameterType; + /// If non-empty, only track states in these geometry regions are used. + /// `TrackState` source only. + std::vector geometrySelection; + /// The Gaussian fit backend used by @c fitProfiles. If unset, /// @c fitProfiles logs a warning and returns no profiles. HistogramFitFunction fitFunction; @@ -63,11 +92,21 @@ class TrackParameterPerformanceCollector { /// Fill histograms for one event. /// + /// @param geoContext the geometry context + /// @param tracks the input tracks + /// @param particles the truth particles + /// @param trackParticleMatching the track to particle matching + /// @param simHits the simulated hits, required for `TrackState` + /// @param measurementSimHitsMap the measurement to simulated hits map, + /// required for `TrackState` + /// /// @note The caller must ensure exclusive access (e.g. hold a mutex). void fill(const Acts::GeometryContext& geoContext, const ConstTrackContainer& tracks, const SimParticleContainer& particles, - const TrackParticleMatching& trackParticleMatching); + const TrackParticleMatching& trackParticleMatching, + const SimHitContainer* simHits = nullptr, + const MeasurementSimHitsMap* measurementSimHitsMap = nullptr); /// Summary count statistics accumulated across all filled events. struct Stats { @@ -76,6 +115,10 @@ class TrackParameterPerformanceCollector { std::size_t nTotalFakeTracks = 0; std::size_t nTotalParticles = 0; std::size_t nTotalMatchedParticles = 0; + /// Track states skipped for lack of the requested parameters. + std::size_t nMissingStateParameters = 0; + /// Track states skipped for lack of truth hits. + std::size_t nMissingStateTruth = 0; }; /// Return accumulated event counts. @@ -111,6 +154,14 @@ class TrackParameterPerformanceCollector { private: const Acts::Logger& logger() const { return *m_logger; } + /// Fill the residuals of the selected measurement states of one track + /// against the truth on their own surfaces. + void fillTrackStates(const Acts::GeometryContext& geoContext, + const ConstTrackProxy& track, + const SimParticle& particle, + const SimHitContainer& simHits, + const MeasurementSimHitsMap& measurementSimHitsMap); + /// Fit every histogram in @p histMap and append the resulting mean/width /// profiles to @p out, warning on excessive fit failures. template @@ -126,6 +177,8 @@ class TrackParameterPerformanceCollector { EffPlotTool m_effPlotTool; TrackSummaryPlotTool m_trackSummaryPlotTool; + Acts::GeometryHierarchyMap m_geometrySelection; + Stats m_stats; }; diff --git a/Examples/Framework/src/Validation/ParametersOnSurface.cpp b/Examples/Framework/src/Validation/ParametersOnSurface.cpp new file mode 100644 index 00000000000..0afe4af967e --- /dev/null +++ b/Examples/Framework/src/Validation/ParametersOnSurface.cpp @@ -0,0 +1,108 @@ +// This file is part of the ACTS project. +// +// Copyright (C) 2016 CERN for the benefit of the ACTS project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at https://mozilla.org/MPL/2.0/. + +#include "ActsExamples/Validation/ParametersOnSurface.hpp" + +#include "Acts/Utilities/TrackHelpers.hpp" +#include "Acts/Utilities/VectorHelpers.hpp" +#include "ActsExamples/EventData/AverageSimHits.hpp" +#include "ActsExamples/Utilities/Range.hpp" + +#include + +std::optional +ActsExamples::truthParametersOnSurface( + const Acts::GeometryContext& gctx, const Acts::Surface& surface, + Index measurementIndex, const SimParticle& particle, + const SimHitContainer& simHits, + const MeasurementSimHitsMap& measurementSimHitsMap, + const Acts::Logger& logger) { + using Acts::VectorHelpers::phi; + using Acts::VectorHelpers::theta; + + using enum Acts::BoundIndices; + + const auto indices = + makeRange(measurementSimHitsMap.equal_range(measurementIndex)); + if (indices.empty()) { + ACTS_WARNING("No truth hits associated to measurement " << measurementIndex + << " found"); + return std::nullopt; + } + + const auto [truthLocal, truthPos4, truthUnitDir] = + averageSimHits(gctx, surface, simHits, indices, logger); + + // position, direction, and time are averaged over the hits above. the + // momentum is not: an average over hits of different particles is not the + // momentum of any of them. take the first hit instead, which exists because + // the range was checked to be non-empty. + const auto simHitIdx0 = indices.begin()->second; + const auto& simHit0 = *simHits.nth(simHitIdx0); + const auto momentum = simHit0.momentum4Before().segment<3>(Acts::eMom0); + + Acts::BoundVector params = Acts::BoundVector::Zero(); + params[eBoundLoc0] = truthLocal[Acts::ePos0]; + params[eBoundLoc1] = truthLocal[Acts::ePos1]; + params[eBoundPhi] = phi(truthUnitDir); + params[eBoundTheta] = theta(truthUnitDir); + params[eBoundQOverP] = + particle.hypothesis().qOverP(momentum.norm(), particle.charge()); + params[eBoundTime] = truthPos4[Acts::eTime]; + + return Acts::BoundTrackParameters(surface.getSharedPtr(), params, + std::nullopt, particle.hypothesis()); +} + +std::optional ActsExamples::recoParametersOnSurface( + const ConstTrackStateProxy& state, + std::optional parameterType, + const Acts::ParticleHypothesis& hypothesis) { + using enum TrackParameterType; + + if (!state.hasReferenceSurface()) { + return std::nullopt; + } + + const auto stateParameters = + [&]() -> std::optional> { + if (!parameterType.has_value()) { + if (!state.hasSmoothed() && !state.hasFiltered() && + !state.hasPredicted()) { + return std::nullopt; + } + // the choice is the proxy's: `parameters()` returns smoothed, else + // filtered, else predicted. `Unbiased` is not among them and has to be + // requested explicitly. + return std::pair(state.parameters(), state.covariance()); + } + if (parameterType == Predicted && state.hasPredicted()) { + return std::pair(state.predicted(), state.predictedCovariance()); + } + if (parameterType == Filtered && state.hasFiltered()) { + return std::pair(state.filtered(), state.filteredCovariance()); + } + if (parameterType == Smoothed && state.hasSmoothed()) { + return std::pair(state.smoothed(), state.smoothedCovariance()); + } + if (parameterType == Unbiased && state.hasSmoothed() && + state.hasProjector() && state.hasCalibrated()) { + return Acts::calculateUnbiasedParametersCovariance( + Acts::AnyConstTrackStateProxy{state}); + } + return std::nullopt; + }(); + + if (!stateParameters.has_value()) { + return std::nullopt; + } + + return Acts::BoundTrackParameters(state.referenceSurface().getSharedPtr(), + stateParameters->first, + stateParameters->second, hypothesis); +} diff --git a/Examples/Framework/src/Validation/ResPlotTool.cpp b/Examples/Framework/src/Validation/ResPlotTool.cpp index befa88a6fc3..d34f9d69fb5 100644 --- a/Examples/Framework/src/Validation/ResPlotTool.cpp +++ b/Examples/Framework/src/Validation/ResPlotTool.cpp @@ -16,6 +16,7 @@ #include #include +#include namespace ActsExamples { @@ -24,10 +25,25 @@ static constexpr double nan = std::numeric_limits::quiet_NaN(); ResPlotTool::ResPlotTool(const ResPlotTool::Config& cfg, Acts::Logging::Level lvl) : m_cfg(cfg), m_logger(Acts::getDefaultLogger("ResPlotTool", lvl)) { - const auto& etaAxis = m_cfg.varBinning.at("Eta"); - const auto& phiAxis = m_cfg.varBinning.at("Phi"); - const auto& ptAxis = m_cfg.varBinning.at("Pt"); - const auto& pullAxis = m_cfg.varBinning.at("Pull"); + // `varBinning.at` would only report the key type, not the missing key + const auto binning = [this](const std::string& key) -> const AxisVariant& { + const auto it = m_cfg.varBinning.find(key); + if (it == m_cfg.varBinning.end()) { + throw std::invalid_argument("ResPlotTool: missing binning for '" + key + + "'"); + } + return it->second; + }; + + const auto& etaAxis = binning("Eta"); + const auto& phiAxis = binning("Phi"); + const auto& ptAxis = binning("Pt"); + const auto& pullAxis = binning("Pull"); + + if (m_cfg.paramNames.size() != Acts::eBoundSize) { + throw std::invalid_argument( + "ResPlotTool: expected one name per bound parameter"); + } ACTS_DEBUG("Initialize the histograms for residual and pull plots"); @@ -36,8 +52,7 @@ ResPlotTool::ResPlotTool(const ResPlotTool::Config& cfg, allParamNames.push_back(m_cfg.relQoverPtName); for (const std::string& parName : allParamNames) { - const std::string parResidual = "Residual_" + parName; - const auto& residualAxis = m_cfg.varBinning.at(parResidual); + const auto& residualAxis = binning("Residual_" + parName); // residual distributions m_res.emplace(parName, Acts::Experimental::Histogram1( @@ -116,11 +131,6 @@ void ResPlotTool::fill(const Acts::GeometryContext& gctx, using enum Acts::BoundIndices; - // get the fitted parameter (at perigee surface) and its error - const Acts::BoundVector& trackParameters = fittedParamters.parameters(); - const Acts::BoundMatrix& trackCovariance = - fittedParamters.covariance().value_or(Acts::BoundMatrix::Zero()); - // get the perigee surface const Acts::Surface& pSurface = fittedParamters.referenceSurface(); @@ -145,16 +155,56 @@ void ResPlotTool::fill(const Acts::GeometryContext& gctx, truthParameters[eBoundQOverP] = truthParticle.qOverP(); truthParameters[eBoundTime] = truthParticle.time(); - // get the truth eta and pT - const double truthEta = eta(truthParticle.direction()); - const double truthPhi = phi(truthParticle.direction()); - const double truthPt = truthParticle.transverseMomentum(); + // bin on the particle, not the bound parameters, which would round-trip the + // direction through phi/theta + fill(truthParameters, + TruthBinning{eta(truthParticle.direction()), + phi(truthParticle.direction()), + truthParticle.transverseMomentum(), truthParticle.charge(), + truthParticle.absoluteCharge()}, + fittedParamters); +} + +void ResPlotTool::fill(const Acts::BoundTrackParameters& truthParameters, + const Acts::BoundTrackParameters& fittedParameters) { + using Acts::VectorHelpers::eta; + using Acts::VectorHelpers::phi; + + if (truthParameters.referenceSurface() != + fittedParameters.referenceSurface()) { + throw std::invalid_argument( + "ResPlotTool: truth and fitted parameters are expressed on different " + "reference surfaces"); + } + + const double truthCharge = truthParameters.charge(); + fill(truthParameters.parameters(), + TruthBinning{eta(truthParameters.direction()), + phi(truthParameters.direction()), + truthParameters.transverseMomentum(), truthCharge, + std::abs(truthCharge)}, + fittedParameters); +} + +void ResPlotTool::fill(const Acts::BoundVector& truthVector, + const TruthBinning& truthBinning, + const Acts::BoundTrackParameters& fittedParameters) { + using enum Acts::BoundIndices; + + const double truthEta = truthBinning.eta; + const double truthPhi = truthBinning.phi; + const double truthPt = truthBinning.pt; + + // get the fitted parameter and its error + const Acts::BoundVector& trackParameters = fittedParameters.parameters(); + const Acts::BoundMatrix& trackCovariance = + fittedParameters.covariance().value_or(Acts::BoundMatrix::Zero()); // fill the histograms for residual and pull for (unsigned int paramId = 0; paramId < Acts::eBoundSize; paramId++) { const std::string& parName = m_cfg.paramNames.at(paramId); - const double residual = trackParameters[paramId] - truthParameters[paramId]; + const double residual = trackParameters[paramId] - truthVector[paramId]; fillResidual(parName, residual, truthEta, truthPhi, truthPt); const double var = trackCovariance(paramId, paramId); @@ -165,8 +215,8 @@ void ResPlotTool::fill(const Acts::GeometryContext& gctx, // `reco(q/pT)` and `true(pT/q) * reco(q/pT)` residual and pull { - const double truthQoverPt = truthParticle.charge() / truthPt; - const double truthPtOverAbsQ = truthPt / truthParticle.absoluteCharge(); + const double truthQoverPt = truthBinning.charge / truthPt; + const double truthPtOverAbsQ = truthPt / truthBinning.absCharge; const double recoQoverPt = trackParameters[eBoundQOverP] / std::sin(trackParameters[eBoundTheta]); const double residualQoverPt = recoQoverPt - truthQoverPt; diff --git a/Examples/Framework/src/Validation/TrackParameterPerformanceCollector.cpp b/Examples/Framework/src/Validation/TrackParameterPerformanceCollector.cpp index 70d5fd2e6cf..8aca4f88c54 100644 --- a/Examples/Framework/src/Validation/TrackParameterPerformanceCollector.cpp +++ b/Examples/Framework/src/Validation/TrackParameterPerformanceCollector.cpp @@ -8,9 +8,12 @@ #include "ActsExamples/Validation/TrackParameterPerformanceCollector.hpp" +#include "Acts/Surfaces/Surface.hpp" #include "Acts/Utilities/Logger.hpp" #include "Acts/Utilities/VectorHelpers.hpp" +#include "ActsExamples/EventData/IndexSourceLink.hpp" +#include #include namespace ActsExamples { @@ -22,12 +25,40 @@ TrackParameterPerformanceCollector::TrackParameterPerformanceCollector( m_resPlotTool(m_cfg.resPlotToolConfig, m_logger->level()), m_effPlotTool(m_cfg.effPlotToolConfig, m_logger->level()), m_trackSummaryPlotTool(m_cfg.trackSummaryPlotToolConfig, - m_logger->level()) {} + m_logger->level()) { + if (m_cfg.parameterSource == TrackParameterSource::Track && + (m_cfg.parameterType.has_value() || !m_cfg.geometrySelection.empty())) { + throw std::invalid_argument( + "Parameter type and geometry selection only apply to the TrackState " + "parameter source"); + } + + std::vector::InputElement> elements; + elements.reserve(m_cfg.geometrySelection.size()); + for (const Acts::GeometryIdentifier& geoId : m_cfg.geometrySelection) { + elements.emplace_back(geoId, 0u); + } + m_geometrySelection = + Acts::GeometryHierarchyMap(std::move(elements)); +} void TrackParameterPerformanceCollector::fill( const Acts::GeometryContext& geoContext, const ConstTrackContainer& tracks, const SimParticleContainer& particles, - const TrackParticleMatching& trackParticleMatching) { + const TrackParticleMatching& trackParticleMatching, + const SimHitContainer* simHits, + const MeasurementSimHitsMap* measurementSimHitsMap) { + // with the track-state source the comparison happens per measurement state + // on the surface that state sits on, so the track itself needs no reference + // surface, but the simulated hits behind the measurements are required + const bool fromTrackStates = + m_cfg.parameterSource == TrackParameterSource::TrackState; + if (fromTrackStates && + (simHits == nullptr || measurementSimHitsMap == nullptr)) { + throw std::invalid_argument( + "Missing simulated hits for the TrackState parameter source"); + } + // Truth particles with corresponding reconstructed tracks std::vector reconParticleIds; reconParticleIds.reserve(tracks.size()); @@ -37,12 +68,10 @@ void TrackParameterPerformanceCollector::fill( ++m_stats.nTotalTracks; // Select reco track with fitted parameters - if (!track.hasReferenceSurface()) { + if (!fromTrackStates && !track.hasReferenceSurface()) { ACTS_DEBUG("No fitted track parameters for track " << track.index()); continue; } - Acts::BoundTrackParameters fittedParameters = - track.createParametersAtReference(); // Get the truth-matched particle auto imatched = trackParticleMatching.find(track.index()); @@ -70,6 +99,14 @@ void TrackParameterPerformanceCollector::fill( // Record this majority particle ID reconParticleIds.push_back(ip->particleId()); + if (fromTrackStates) { + fillTrackStates(geoContext, track, *ip, *simHits, *measurementSimHitsMap); + continue; + } + + Acts::BoundTrackParameters fittedParameters = + track.createParametersAtReference(); + // Fill residual plots m_resPlotTool.fill(geoContext, ip->initialState(), fittedParameters); @@ -107,6 +144,54 @@ void TrackParameterPerformanceCollector::fill( } } +void TrackParameterPerformanceCollector::fillTrackStates( + const Acts::GeometryContext& geoContext, const ConstTrackProxy& track, + const SimParticle& particle, const SimHitContainer& simHits, + const MeasurementSimHitsMap& measurementSimHitsMap) { + for (const auto& state : track.trackStatesReversed()) { + if (!state.typeFlags().isMeasurement() || state.typeFlags().isOutlier()) { + continue; + } + if (!state.hasReferenceSurface()) { + continue; + } + const Acts::Surface& surface = state.referenceSurface(); + + if (!m_geometrySelection.empty() && + m_geometrySelection.find(surface.geometryId()) == + m_geometrySelection.end()) { + continue; + } + + const std::optional reco = + recoParametersOnSurface(state, m_cfg.parameterType, + track.particleHypothesis()); + if (!reco.has_value()) { + ++m_stats.nMissingStateParameters; + continue; + } + + // the source link must outlive the pointer into it + const Acts::SourceLink sourceLink = state.getUncalibratedSourceLink(); + const auto* indexSourceLink = sourceLink.getPtr(); + if (indexSourceLink == nullptr) { + ++m_stats.nMissingStateTruth; + continue; + } + + const std::optional truth = + truthParametersOnSurface(geoContext, surface, indexSourceLink->index(), + particle, simHits, measurementSimHitsMap, + logger()); + if (!truth.has_value()) { + ++m_stats.nMissingStateTruth; + continue; + } + + m_resPlotTool.fill(truth.value(), reco.value()); + } +} + void TrackParameterPerformanceCollector::logSummary() const { ACTS_INFO("=== Track Parameter Performance Summary ==="); ACTS_INFO("Total tracks: " << m_stats.nTotalTracks); @@ -114,6 +199,16 @@ void TrackParameterPerformanceCollector::logSummary() const { ACTS_INFO("Total particles: " << m_stats.nTotalParticles); ACTS_INFO("Total matched particles: " << m_stats.nTotalMatchedParticles); + if (m_cfg.parameterSource == TrackParameterSource::TrackState) { + // a state counts here when it does not carry the requested parameters at + // all, which is not a failure per se: an input that stores its estimate on + // a single state, e.g. seeding output, skips every other state of a track + ACTS_INFO("Skipped states without the requested parameters: " + << m_stats.nMissingStateParameters); + ACTS_INFO( + "Skipped states without truth hits: " << m_stats.nMissingStateTruth); + } + if (m_stats.nTotalTracks > 0) { double efficiency = static_cast(m_stats.nTotalMatchedTracks) / m_stats.nTotalTracks; diff --git a/Examples/Io/Root/include/ActsExamples/Io/Root/RootTrackParameterPerformanceWriter.hpp b/Examples/Io/Root/include/ActsExamples/Io/Root/RootTrackParameterPerformanceWriter.hpp index 132c04e1163..9470778a61c 100644 --- a/Examples/Io/Root/include/ActsExamples/Io/Root/RootTrackParameterPerformanceWriter.hpp +++ b/Examples/Io/Root/include/ActsExamples/Io/Root/RootTrackParameterPerformanceWriter.hpp @@ -9,6 +9,7 @@ #pragma once #include "Acts/Utilities/Logger.hpp" +#include "ActsExamples/EventData/SimHit.hpp" #include "ActsExamples/EventData/SimParticle.hpp" #include "ActsExamples/EventData/Track.hpp" #include "ActsExamples/EventData/TruthMatching.hpp" @@ -21,7 +22,9 @@ #include "ActsExamples/Validation/TrackSummaryPlotTool.hpp" #include +#include #include +#include class TFile; class TTree; @@ -32,6 +35,11 @@ namespace ActsExamples { /// /// Efficiency here is the fraction of smoothed tracks compared to all tracks. /// +/// With `parameterSource = TrackState` the residuals are taken from the +/// individual measurement states on their own surfaces instead of the track +/// reference surface, optionally restricted to a geometry region. That needs +/// `inputSimHits` and `inputMeasurementSimHitsMap` for the truth. +/// /// A common file can be provided for the writer to attach his TTree, /// this is done by setting the Config::rootFile pointer to an existing file /// @@ -46,6 +54,10 @@ class RootTrackParameterPerformanceWriter final std::string inputParticles; /// Input track-particle matching. std::string inputTrackParticleMatching; + /// Input simulated hits collection. `TrackState` source only. + std::string inputSimHits; + /// Input measurement to simulated hits map. `TrackState` source only. + std::string inputMeasurementSimHitsMap; /// Output filename. std::string filePath = "performance_track_parameters.root"; /// Plot tool configurations. @@ -53,6 +65,15 @@ class RootTrackParameterPerformanceWriter final EffPlotTool::Config effPlotToolConfig; TrackSummaryPlotTool::Config trackSummaryPlotToolConfig; + /// Where to take the reconstructed parameters from. + TrackParameterSource parameterSource = TrackParameterSource::Track; + /// Which track-state parameters to use. If not set, the best available + /// ones (smoothed, filtered, or predicted). `TrackState` source only. + std::optional parameterType; + /// If non-empty, only track states in these geometry regions are used. + /// `TrackState` source only. + std::vector geometrySelection; + /// Minimum number of entries in a bin for it to be included in the /// mean/width fit. int fitMinEntries = 10; @@ -87,6 +108,9 @@ class RootTrackParameterPerformanceWriter final ReadDataHandle m_inputParticles{this, "InputParticles"}; ReadDataHandle m_inputTrackParticleMatching{ this, "InputTrackParticleMatching"}; + ReadDataHandle m_inputSimHits{this, "InputSimHits"}; + ReadDataHandle m_inputMeasurementSimHitsMap{ + this, "InputMeasurementSimHitsMap"}; /// Mutex used to protect multi-threaded writes. std::mutex m_writeMutex; diff --git a/Examples/Io/Root/src/RootTrackParameterPerformanceWriter.cpp b/Examples/Io/Root/src/RootTrackParameterPerformanceWriter.cpp index c047960cd1f..0746e5974c5 100644 --- a/Examples/Io/Root/src/RootTrackParameterPerformanceWriter.cpp +++ b/Examples/Io/Root/src/RootTrackParameterPerformanceWriter.cpp @@ -28,18 +28,35 @@ using ActsPlugins::toRoot; namespace ActsExamples { +namespace { + +/// Translate the writer configuration into the collector configuration. +TrackParameterPerformanceCollector::Config collectorConfig( + const RootTrackParameterPerformanceWriter::Config& cfg) { + TrackParameterPerformanceCollector::Config collectorCfg; + collectorCfg.resPlotToolConfig = cfg.resPlotToolConfig; + collectorCfg.effPlotToolConfig = cfg.effPlotToolConfig; + collectorCfg.trackSummaryPlotToolConfig = cfg.trackSummaryPlotToolConfig; + collectorCfg.parameterSource = cfg.parameterSource; + collectorCfg.parameterType = cfg.parameterType; + collectorCfg.geometrySelection = cfg.geometrySelection; + collectorCfg.fitFunction = ActsPlugins::RootHistogramFit(); + collectorCfg.fitMinEntries = cfg.fitMinEntries; + collectorCfg.fitSigmaRange = cfg.fitSigmaRange; + collectorCfg.fitIterations = cfg.fitIterations; + collectorCfg.warningThresholdFitFailureFraction = + cfg.warningThresholdFitFailureFraction; + return collectorCfg; +} + +} // namespace + RootTrackParameterPerformanceWriter::RootTrackParameterPerformanceWriter( RootTrackParameterPerformanceWriter::Config config, Acts::Logging::Level level) : WriterT(config.inputTracks, "RootTrackParameterPerformanceWriter", level), m_cfg(std::move(config)), - m_collector( - TrackParameterPerformanceCollector::Config{ - m_cfg.resPlotToolConfig, m_cfg.effPlotToolConfig, - m_cfg.trackSummaryPlotToolConfig, ActsPlugins::RootHistogramFit(), - m_cfg.fitMinEntries, m_cfg.fitSigmaRange, m_cfg.fitIterations, - m_cfg.warningThresholdFitFailureFraction}, - logger().clone()) { + m_collector(collectorConfig(m_cfg), logger().clone()) { // trajectories collection name is already checked by base ctor if (m_cfg.inputParticles.empty()) { throw std::invalid_argument("Missing particles input collection"); @@ -54,6 +71,18 @@ RootTrackParameterPerformanceWriter::RootTrackParameterPerformanceWriter( m_inputParticles.initialize(m_cfg.inputParticles); m_inputTrackParticleMatching.initialize(m_cfg.inputTrackParticleMatching); + if (m_cfg.parameterSource == TrackParameterSource::TrackState) { + if (m_cfg.inputSimHits.empty()) { + throw std::invalid_argument("Missing simulated hits input collection"); + } + if (m_cfg.inputMeasurementSimHitsMap.empty()) { + throw std::invalid_argument("Missing measurement to simulated hits map"); + } + + m_inputSimHits.initialize(m_cfg.inputSimHits); + m_inputMeasurementSimHitsMap.initialize(m_cfg.inputMeasurementSimHitsMap); + } + // the output file can not be given externally since TFile accesses to the // same file from multiple threads are unsafe. // must always be opened internally @@ -158,11 +187,18 @@ ProcessCode RootTrackParameterPerformanceWriter::writeT( const auto& particles = m_inputParticles(ctx); const auto& trackParticleMatching = m_inputTrackParticleMatching(ctx); + const SimHitContainer* simHits = nullptr; + const MeasurementSimHitsMap* measurementSimHitsMap = nullptr; + if (m_cfg.parameterSource == TrackParameterSource::TrackState) { + simHits = &m_inputSimHits(ctx); + measurementSimHitsMap = &m_inputMeasurementSimHitsMap(ctx); + } + // Exclusive access to the histograms while filling std::lock_guard lock(m_writeMutex); - m_collector.fill(ctx.recoGeoContext, tracks, particles, - trackParticleMatching); + m_collector.fill(ctx.recoGeoContext, tracks, particles, trackParticleMatching, + simHits, measurementSimHitsMap); return ProcessCode::SUCCESS; } diff --git a/Python/Examples/src/Framework.cpp b/Python/Examples/src/Framework.cpp index c0f5e94e051..1d2ef0daf1e 100644 --- a/Python/Examples/src/Framework.cpp +++ b/Python/Examples/src/Framework.cpp @@ -21,6 +21,7 @@ #include "ActsExamples/Validation/EffPlotTool.hpp" #include "ActsExamples/Validation/FakePlotTool.hpp" #include "ActsExamples/Validation/ResPlotTool.hpp" +#include "ActsExamples/Validation/TrackParameterPerformanceCollector.hpp" #include "ActsExamples/Validation/TrackQualityPlotTool.hpp" #include "ActsExamples/Validation/TrackSummaryPlotTool.hpp" #include "ActsPython/Utilities/Macros.hpp" @@ -508,6 +509,16 @@ void addFramework(py::module& mex) { .def_readwrite("paramNames", &ResPlotTool::Config::paramNames) .def_readwrite("varBinning", &ResPlotTool::Config::varBinning); + py::enum_(mex, "TrackParameterSource") + .value("Track", TrackParameterSource::Track) + .value("TrackState", TrackParameterSource::TrackState); + + py::enum_(mex, "TrackParameterType") + .value("Predicted", TrackParameterType::Predicted) + .value("Filtered", TrackParameterType::Filtered) + .value("Smoothed", TrackParameterType::Smoothed) + .value("Unbiased", TrackParameterType::Unbiased); + py::class_(mex, "TrackQualityPlotToolConfig") .def(py::init<>()) .def_readwrite("varBinning", &TrackQualityPlotTool::Config::varBinning); diff --git a/Python/Examples/src/PythonSpecific.cpp b/Python/Examples/src/PythonSpecific.cpp index 26f1e3b39e5..7b3ae2aa35a 100644 --- a/Python/Examples/src/PythonSpecific.cpp +++ b/Python/Examples/src/PythonSpecific.cpp @@ -213,16 +213,26 @@ class PythonTrackParameterPerformanceWriter final double warningThresholdFitFailureFraction = 0.55; }; + /// Translate the writer configuration into the collector configuration. + static TrackParameterPerformanceCollector::Config collectorConfig( + const Config& cfg) { + TrackParameterPerformanceCollector::Config collectorCfg; + collectorCfg.resPlotToolConfig = cfg.resPlotToolConfig; + collectorCfg.effPlotToolConfig = cfg.effPlotToolConfig; + collectorCfg.trackSummaryPlotToolConfig = cfg.trackSummaryPlotToolConfig; + collectorCfg.fitFunction = cfg.fitFunction; + collectorCfg.fitMinEntries = cfg.fitMinEntries; + collectorCfg.fitSigmaRange = cfg.fitSigmaRange; + collectorCfg.fitIterations = cfg.fitIterations; + collectorCfg.warningThresholdFitFailureFraction = + cfg.warningThresholdFitFailureFraction; + return collectorCfg; + } + PythonTrackParameterPerformanceWriter(Config cfg, Acts::Logging::Level lvl) : WriterT(cfg.inputTracks, "PythonTrackParameterPerformanceWriter", lvl), m_cfg(std::move(cfg)), - m_collector( - TrackParameterPerformanceCollector::Config{ - m_cfg.resPlotToolConfig, m_cfg.effPlotToolConfig, - m_cfg.trackSummaryPlotToolConfig, m_cfg.fitFunction, - m_cfg.fitMinEntries, m_cfg.fitSigmaRange, m_cfg.fitIterations, - m_cfg.warningThresholdFitFailureFraction}, - logger().clone()) { + m_collector(collectorConfig(m_cfg), logger().clone()) { if (m_cfg.inputParticles.empty()) { throw std::invalid_argument("Missing particles input collection"); } diff --git a/Python/Examples/src/plugins/Root.cpp b/Python/Examples/src/plugins/Root.cpp index 2798fdf9006..274db73f38b 100644 --- a/Python/Examples/src/plugins/Root.cpp +++ b/Python/Examples/src/plugins/Root.cpp @@ -45,8 +45,6 @@ #include "ActsPlugins/Root/RootHistogramFit.hpp" #include "ActsPython/Utilities/Macros.hpp" -#include - #include #include #include @@ -136,9 +134,11 @@ PYBIND11_MODULE(ActsExamplesPythonBindingsRoot, root) { ACTS_PYTHON_DECLARE_WRITER( RootTrackParameterPerformanceWriter, root, "RootTrackParameterPerformanceWriter", inputTracks, inputParticles, - inputTrackParticleMatching, filePath, resPlotToolConfig, - effPlotToolConfig, trackSummaryPlotToolConfig, fitMinEntries, - fitSigmaRange, fitIterations, warningThresholdFitFailureFraction); + inputTrackParticleMatching, inputSimHits, inputMeasurementSimHitsMap, + filePath, resPlotToolConfig, effPlotToolConfig, + trackSummaryPlotToolConfig, parameterSource, parameterType, + geometrySelection, fitMinEntries, fitSigmaRange, fitIterations, + warningThresholdFitFailureFraction); ACTS_PYTHON_DECLARE_WRITER( RootTrackParameterWriter, root, "RootTrackParameterWriter",