feat: Collect and write material statistics in GSF (#2695)

This refactors the collection of material statistics in the GSF, so it only accumulates values until the last measurement in the forward pass (not collecting material in the solenoid etc...)

Also writes the sum and the maximum material of the forward pass to the tracksummary.

Core/include/Acts/TrackFitting/GaussianSumFitter.hpp

@@ -305,7 +305,7 @@ struct GaussianSumFitter {
       return return_error_or_abort(fwdResult.error());
     }
 
-    auto& fwdGsfResult =
+    const auto& fwdGsfResult =
         fwdResult->template get<typename GsfActor::result_type>();
 
     if (!fwdGsfResult.result.ok()) {
@@ -321,8 +321,8 @@ struct GaussianSumFitter {
     ACTS_VERBOSE("- processed states: " << fwdGsfResult.processedStates);
     ACTS_VERBOSE("- measurement states: " << fwdGsfResult.measurementStates);
 
-    std::size_t nInvalidBetheHeitler = fwdGsfResult.nInvalidBetheHeitler;
-    double maxPathXOverX0 = fwdGsfResult.maxPathXOverX0;
+    std::size_t nInvalidBetheHeitler = fwdGsfResult.nInvalidBetheHeitler.val();
+    double maxPathXOverX0 = fwdGsfResult.maxPathXOverX0.val();
 
     //////////////////
     // Backward pass
@@ -408,8 +408,14 @@ struct GaussianSumFitter {
           GsfError::NoMeasurementStatesCreatedBackward);
     }
 
-    nInvalidBetheHeitler += bwdGsfResult.nInvalidBetheHeitler;
-    maxPathXOverX0 = std::max(maxPathXOverX0, bwdGsfResult.maxPathXOverX0);
+    // For the backward pass we want the counters at in end (= at the
+    // interaction point) and not at the last measurement surface
+    bwdGsfResult.nInvalidBetheHeitler.update();
+    bwdGsfResult.maxPathXOverX0.update();
+    bwdGsfResult.sumPathXOverX0.update();
+    nInvalidBetheHeitler += bwdGsfResult.nInvalidBetheHeitler.val();
+    maxPathXOverX0 =
+        std::max(maxPathXOverX0, bwdGsfResult.maxPathXOverX0.val());
 
     if (nInvalidBetheHeitler > 0) {
       ACTS_WARNING("Encountered " << nInvalidBetheHeitler
@@ -479,12 +485,23 @@ struct GaussianSumFitter {
 
       if (trackContainer.hasColumn(
               hashString(GsfConstants::kFinalMultiComponentStateColumn))) {
-        ACTS_DEBUG("Add final multi-component state to track")
+        ACTS_DEBUG("Add final multi-component state to track");
         track.template component<GsfConstants::FinalMultiComponentState>(
             GsfConstants::kFinalMultiComponentStateColumn) = std::move(params);
       }
     }
 
+    if (trackContainer.hasColumn(
+            hashString(GsfConstants::kFwdMaxMaterialXOverX0))) {
+      track.template component<double>(GsfConstants::kFwdMaxMaterialXOverX0) =
+          fwdGsfResult.maxPathXOverX0.val();
+    }
+    if (trackContainer.hasColumn(
+            hashString(GsfConstants::kFwdSumMaterialXOverX0))) {
+      track.template component<double>(GsfConstants::kFwdSumMaterialXOverX0) =
+          fwdGsfResult.sumPathXOverX0.val();
+    }
+
     calculateTrackQuantities(track);
 
     return track;

Core/include/Acts/TrackFitting/GsfOptions.hpp

@@ -39,6 +39,10 @@ constexpr std::string_view kFinalMultiComponentStateColumn =
     "gsf-final-multi-component-state";
 using FinalMultiComponentState =
     std::optional<Acts::MultiComponentBoundTrackParameters>;
+constexpr std::string_view kFwdSumMaterialXOverX0 =
+    "gsf-fwd-sum-material-x-over-x0";
+constexpr std::string_view kFwdMaxMaterialXOverX0 =
+    "gsf-fwd-max-material-x-over-x0";
 }  // namespace GsfConstants
 
 /// The extensions needed for the GSF

Core/include/Acts/TrackFitting/detail/GsfActor.hpp

@@ -56,9 +56,10 @@ struct GsfResult {
   std::vector<const Acts::Surface*> visitedSurfaces;
   std::vector<const Acts::Surface*> surfacesVisitedBwdAgain;
 
-  /// Statistics about encounterings of invalid bethe heitler approximations
-  std::size_t nInvalidBetheHeitler = 0;
-  double maxPathXOverX0 = 0.0;
+  /// Statistics about material encounterings
+  Updatable<std::size_t> nInvalidBetheHeitler;
+  Updatable<double> maxPathXOverX0;
+  Updatable<double> sumPathXOverX0;
 
   // Propagate potential errors to the outside
   Result<void> result{Result<void>::success()};
@@ -232,6 +233,15 @@ struct GsfActor {
       return;
     }
 
+    // Update the counters. Note that this should be done before potential
+    // material interactions, because if this is our last measurement this would
+    // not influence the fit anymore.
+    if (haveMeasurement) {
+      result.maxPathXOverX0.update();
+      result.sumPathXOverX0.update();
+      result.nInvalidBetheHeitler.update();
+    }
+
     for (auto cmp : stepper.componentIterable(state.stepping)) {
       auto singleState = cmp.singleState(state);
       cmp.singleStepper(stepper).transportCovarianceToBound(
@@ -332,25 +342,31 @@ struct GsfActor {
                            std::vector<ComponentCache>& componentCache,
                            result_type& result) const {
     auto cmps = stepper.componentIterable(state.stepping);
+    double pathXOverX0 = 0.0;
     for (auto [idx, cmp] : zip(tmpStates.tips, cmps)) {
       auto proxy = tmpStates.traj.getTrackState(idx);
 
       BoundTrackParameters bound(proxy.referenceSurface().getSharedPtr(),
                                  proxy.filtered(), proxy.filteredCovariance(),
                                  stepper.particleHypothesis(state.stepping));
 
-      applyBetheHeitler(state, navigator, bound, tmpStates.weights.at(idx),
-                        componentCache, result);
+      pathXOverX0 +=
+          applyBetheHeitler(state, navigator, bound, tmpStates.weights.at(idx),
+                            componentCache, result);
     }
+
+    // Store average material seen by the components
+    // Should not be too broadly distributed
+    result.sumPathXOverX0.tmp() += pathXOverX0 / tmpStates.tips.size();
   }
 
   template <typename propagator_state_t, typename navigator_t>
-  void applyBetheHeitler(const propagator_state_t& state,
-                         const navigator_t& navigator,
-                         const BoundTrackParameters& old_bound,
-                         const double old_weight,
-                         std::vector<ComponentCache>& componentCaches,
-                         result_type& result) const {
+  double applyBetheHeitler(const propagator_state_t& state,
+                           const navigator_t& navigator,
+                           const BoundTrackParameters& old_bound,
+                           const double old_weight,
+                           std::vector<ComponentCache>& componentCaches,
+                           result_type& result) const {
     const auto& surface = *navigator.currentSurface(state.navigation);
     const auto p_prev = old_bound.absoluteMomentum();
 
@@ -365,11 +381,12 @@ struct GsfActor {
     slab.scaleThickness(pathCorrection);
 
     const double pathXOverX0 = slab.thicknessInX0();
+    result.maxPathXOverX0.tmp() =
+        std::max(result.maxPathXOverX0.tmp(), pathXOverX0);
 
     // Emit a warning if the approximation is not valid for this x/x0
     if (!m_cfg.bethe_heitler_approx->validXOverX0(pathXOverX0)) {
-      ++result.nInvalidBetheHeitler;
-      result.maxPathXOverX0 = std::max(result.maxPathXOverX0, pathXOverX0);
+      ++result.nInvalidBetheHeitler.tmp();
       ACTS_DEBUG(
           "Bethe-Heitler approximation encountered invalid value for x/x0="
           << pathXOverX0 << " at surface " << surface.geometryId());
@@ -428,6 +445,8 @@ struct GsfActor {
       // Set the remaining things and push to vector
       componentCaches.push_back({new_weight, new_pars, new_cov});
     }
+
+    return pathXOverX0;
   }
 
   /// Remove components with low weights and renormalize from the component

Core/include/Acts/TrackFitting/detail/GsfUtils.hpp

@@ -239,5 +239,22 @@ struct MultiTrajectoryProjector {
   }
 };
 
+/// Small Helper class that allows to carry a temporary value until we decide to
+/// update the actual value. The temporary value is deliberately only accessible
+/// with a mutable reference
+template <typename T>
+class Updatable {
+  T m_tmp{};
+  T m_val{};
+
+ public:
+  Updatable() : m_tmp(0), m_val(0) {}
+
+  T &tmp() { return m_tmp; }
+  void update() { m_val = m_tmp; }
+
+  const T &val() const { return m_val; }
+};
+
 }  // namespace detail
 }  // namespace Acts

Examples/Algorithms/TrackFitting/src/GsfFitterFunction.cpp

@@ -78,6 +78,7 @@ struct GsfFitterFunctionImpl final : public ActsExamples::TrackFitterFunction {
 
   std::size_t maxComponents = 0;
   double weightCutoff = 0;
+  const double momentumCutoff = 0;  // 500_MeV;
   bool abortOnError = false;
   bool disableAllMaterialHandling = false;
   MixtureReductionAlgorithm reductionAlg =
@@ -146,6 +147,14 @@ struct GsfFitterFunctionImpl final : public ActsExamples::TrackFitterFunction {
       tracks.template addColumn<FinalMultiComponentState>(key);
     }
 
+    if (!tracks.hasColumn(Acts::hashString(kFwdMaxMaterialXOverX0))) {
+      tracks.template addColumn<double>(std::string(kFwdMaxMaterialXOverX0));
+    }
+
+    if (!tracks.hasColumn(Acts::hashString(kFwdSumMaterialXOverX0))) {
+      tracks.template addColumn<double>(std::string(kFwdSumMaterialXOverX0));
+    }
+
     return fitter.fit(sourceLinks.begin(), sourceLinks.end(), initialParameters,
                       gsfOptions, tracks);
   }

Examples/Io/Root/include/ActsExamples/Io/Root/RootTrackSummaryWriter.hpp

@@ -67,6 +67,8 @@ class RootTrackSummaryWriter final : public WriterT<ConstTrackContainer> {
     std::string fileMode = "RECREATE";
     /// Switch for adding full covariance matrix to output file.
     bool writeCovMat = false;
+    /// Write GSF specific things (for now only some material statistics)
+    bool writeGsfSpecific = false;
   };
 
   /// Constructor
@@ -226,6 +228,9 @@ class RootTrackSummaryWriter final : public WriterT<ConstTrackContainer> {
   std::vector<float> m_cov_eT_eTHETA;
   std::vector<float> m_cov_eT_eQOP;
   std::vector<float> m_cov_eT_eT;
+
+  std::vector<float> m_gsf_max_material_fwd;
+  std::vector<float> m_gsf_sum_material_fwd;
 };
 
 }  // namespace ActsExamples

Examples/Io/Root/src/RootTrackSummaryWriter.cpp

@@ -13,6 +13,7 @@
 #include "Acts/EventData/MultiTrajectoryHelpers.hpp"
 #include "Acts/EventData/VectorMultiTrajectory.hpp"
 #include "Acts/Surfaces/Surface.hpp"
+#include "Acts/TrackFitting/GsfOptions.hpp"
 #include "Acts/Utilities/Intersection.hpp"
 #include "Acts/Utilities/MultiIndex.hpp"
 #include "Acts/Utilities/Result.hpp"
@@ -137,6 +138,12 @@ ActsExamples::RootTrackSummaryWriter::RootTrackSummaryWriter(
     m_outputTree->Branch("pull_eTHETA_fit", &m_pull_eTHETA_fit);
     m_outputTree->Branch("pull_eQOP_fit", &m_pull_eQOP_fit);
     m_outputTree->Branch("pull_eT_fit", &m_pull_eT_fit);
+
+    if (m_cfg.writeGsfSpecific) {
+      m_outputTree->Branch("max_material_fwd", &m_gsf_max_material_fwd);
+      m_outputTree->Branch("sum_material_fwd", &m_gsf_sum_material_fwd);
+    }
+
     if (m_cfg.writeCovMat == true) {
       // create one branch for every entry of covariance matrix
       // one block for every row of the matrix, every entry gets own branch
@@ -440,6 +447,24 @@ ActsExamples::ProcessCode ActsExamples::RootTrackSummaryWriter::writeT(
     m_pull_eT_fit.push_back(pull[Acts::eBoundTime]);
 
     m_hasFittedParams.push_back(hasFittedParams);
+
+    if (m_cfg.writeGsfSpecific) {
+      using namespace Acts::GsfConstants;
+      if (tracks.hasColumn(Acts::hashString(kFwdMaxMaterialXOverX0))) {
+        m_gsf_max_material_fwd.push_back(
+            track.template component<double>(kFwdMaxMaterialXOverX0));
+      } else {
+        m_gsf_max_material_fwd.push_back(NaNfloat);
+      }
+
+      if (tracks.hasColumn(Acts::hashString(kFwdSumMaterialXOverX0))) {
+        m_gsf_sum_material_fwd.push_back(
+            track.template component<double>(kFwdSumMaterialXOverX0));
+      } else {
+        m_gsf_sum_material_fwd.push_back(NaNfloat);
+      }
+    }
+
     if (m_cfg.writeCovMat) {
       // write all entries of covariance matrix to output file
       // one branch for every entry of the matrix.
@@ -549,6 +574,9 @@ ActsExamples::ProcessCode ActsExamples::RootTrackSummaryWriter::writeT(
   m_pull_eQOP_fit.clear();
   m_pull_eT_fit.clear();
 
+  m_gsf_max_material_fwd.clear();
+  m_gsf_sum_material_fwd.clear();
+
   if (m_cfg.writeCovMat) {
     m_cov_eLOC0_eLOC0.clear();
     m_cov_eLOC0_eLOC1.clear();

Examples/Python/src/Output.cpp

@@ -323,10 +323,10 @@ void addOutput(Context& ctx) {
       inputTracks, inputParticles, inputSimHits, inputMeasurementParticlesMap,
       inputMeasurementSimHitsMap, filePath, treeName, fileMode);
 
-  ACTS_PYTHON_DECLARE_WRITER(ActsExamples::RootTrackSummaryWriter, mex,
-                             "RootTrackSummaryWriter", inputTracks,
-                             inputParticles, inputMeasurementParticlesMap,
-                             filePath, treeName, fileMode, writeCovMat);
+  ACTS_PYTHON_DECLARE_WRITER(
+      ActsExamples::RootTrackSummaryWriter, mex, "RootTrackSummaryWriter",
+      inputTracks, inputParticles, inputMeasurementParticlesMap, filePath,
+      treeName, fileMode, writeCovMat, writeGsfSpecific);
 
   ACTS_PYTHON_DECLARE_WRITER(
       ActsExamples::VertexPerformanceWriter, mex, "VertexPerformanceWriter",

Examples/Python/tests/root_file_hashes.txt

@@ -35,9 +35,9 @@ test_truth_tracking_kalman[odd-1000.0]__trackstates_fitter.root: f616b2234db239b
 test_truth_tracking_kalman[odd-1000.0]__tracksummary_fitter.root: a53253b509b5779a5d856f8c09d76a499b217b55ba4b0e52d9076ffad726463f
 test_truth_tracking_kalman[odd-1000.0]__performance_track_finder.root: 39aec6316cceb90e314e16b02947faa691c18f57c3a851a25e547a8fc05a4593
 test_truth_tracking_gsf[generic]__trackstates_gsf.root: 8422ad83843d69aac7d1e37d6d3e4019b1abfaf854d9246c1fb8a4fc5d6c3f13
-test_truth_tracking_gsf[generic]__tracksummary_gsf.root: fca19c92a21bec029325189d7a404e5da4851b6cca886170e14486863b86c29e
+test_truth_tracking_gsf[generic]__tracksummary_gsf.root: fe5bce60f13f6221290a68340772beefe5bb5389f424b23c14787a944aa5d84f
 test_truth_tracking_gsf[odd]__trackstates_gsf.root: 9e1c0acaf1aed2347992c0a4c076ca27ce71d770521b33680f30f9734dd8494e
-test_truth_tracking_gsf[odd]__tracksummary_gsf.root: ca78cba367449062e985ee5aacd4db88ef721cbe56a772b818f1586770ec05ba
+test_truth_tracking_gsf[odd]__tracksummary_gsf.root: 09e06b254b18e18cb3086038e2f8b7ac96a33d980ac14395bc9a02227c4433bd
 test_particle_gun__particles.root: 8549ba6e20338004ab8ba299fc65e1ee5071985b46df8f77f887cb6fef56a8ec
 test_material_mapping__material-map_tracks.root: 2deac7a48ff1185ba3889cfd4cc0bd0a6de57293048dfd6766f3c3907232e45e
 test_material_mapping__propagation-material.root: 84b04ebd5721550867f0f193b5eb4e9f3f205df542cb46090d320be6bff565c5

Examples/Scripts/Python/truth_tracking_gsf.py

@@ -139,6 +139,7 @@ def runTruthTrackingGsf(
             inputParticles="truth_seeds_selected",
             inputMeasurementParticlesMap="measurement_particles_map",
             filePath=str(outputDir / "tracksummary_gsf.root"),
+            writeGsfSpecific=True,
         )
     )