feat: Write track parameter performance from track states in a geometry region

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

Examples/Framework/include/ActsExamples/Validation/ParametersOnSurface.hpp

@@ -0,0 +1,79 @@
+// 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 <optional>
+
+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<Acts::BoundTrackParameters> 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)
+/// @param hypothesis The particle hypothesis for the parameters
+/// @return The parameters with covariance, or nullopt without a reference
+///         surface or matching parameters
+std::optional<Acts::BoundTrackParameters> recoParametersOnSurface(
+    const ConstTrackStateProxy& state,
+    std::optional<TrackParameterType> parameterType,
+    const Acts::ParticleHypothesis& hypothesis);
+
+}  // namespace ActsExamples

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<std::string, Histogram1>& res() const { return m_res; }
   const std::map<std::string, Histogram2>& 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<const Acts::Logger> m_logger;
@@ -123,6 +140,10 @@ class ResPlotTool {
   /// Pull vs eta-pT scatter plot
   std::map<std::string, Histogram3> 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,

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 <cstddef>
 #include <map>
+#include <optional>
 #include <string>
 #include <vector>
 
 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<TrackParameterType> parameterType;
+    /// If non-empty, only track states in these geometry regions are used.
+    /// `TrackState` source only.
+    std::vector<Acts::GeometryIdentifier> 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 <std::size_t Dim>
@@ -126,6 +177,8 @@ class TrackParameterPerformanceCollector {
   EffPlotTool m_effPlotTool;
   TrackSummaryPlotTool m_trackSummaryPlotTool;
 
+  Acts::GeometryHierarchyMap<unsigned int> m_geometrySelection;
+
   Stats m_stats;
 };
 

Examples/Framework/src/Validation/ParametersOnSurface.cpp

@@ -0,0 +1,104 @@
+// 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 <utility>
+
+std::optional<Acts::BoundTrackParameters>
+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);
+
+  // averaging the momentum makes even less sense than position and direction,
+  // so take the first one. the indices are known to be valid.
+  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<Acts::BoundTrackParameters> ActsExamples::recoParametersOnSurface(
+    const ConstTrackStateProxy& state,
+    std::optional<TrackParameterType> parameterType,
+    const Acts::ParticleHypothesis& hypothesis) {
+  using enum TrackParameterType;
+
+  if (!state.hasReferenceSurface()) {
+    return std::nullopt;
+  }
+
+  const auto stateParameters =
+      [&]() -> std::optional<std::pair<Acts::BoundVector, Acts::BoundMatrix>> {
+    if (!parameterType.has_value()) {
+      if (!state.hasSmoothed() && !state.hasFiltered() &&
+          !state.hasPredicted()) {
+        return std::nullopt;
+      }
+      // best available parameters, i.e. smoothed, filtered, or predicted
+      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);
+}