Merge branch 'main' into RecomputeParametersFromSeed

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/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepUtil.hpp

@@ -118,5 +118,19 @@ void writeTrajectory(const Acts::GeometryContext& gctx, double Bz,
                      const Acts::ParticleHypothesis& particleHypothesis,
                      const IndexMultimap<ActsFatras::Barcode>& hitParticlesMap);
 
+/// Helper function to either return an id as is, or unpack an index from it
+/// if it is a podio::ObjectID.
+/// @tparam T The type of the id.
+/// @param o The id to convert.
+/// @return The id as an unsigned integer.
+template <typename T>
+uint64_t podioObjectIDToInteger(T&& o) {
+  if constexpr (!std::is_same_v<T, podio::ObjectID>) {
+    return o;
+  } else {
+    return o.index;
+  }
+}
+
 }  // namespace EDM4hepUtil
 }  // namespace ActsExamples

Examples/Io/EDM4hep/src/EDM4hepParticleReader.cpp

@@ -55,8 +55,7 @@ ProcessCode EDM4hepParticleReader::read(const AlgorithmContext& ctx) {
     auto particle =
         EDM4hepUtil::readParticle(mcParticle, [](const edm4hep::MCParticle& p) {
           ActsFatras::Barcode result;
-          // TODO dont use podio internal id
-          result.setParticle(p.id());
+          result.setParticle(EDM4hepUtil::podioObjectIDToInteger(p.id()));
           return result;
         });
     unordered.push_back(particle);

Examples/Io/EDM4hep/src/EDM4hepSimHitReader.cpp

@@ -54,8 +54,7 @@ ProcessCode EDM4hepSimHitReader::read(const AlgorithmContext& ctx) {
       auto particle = EDM4hepUtil::readParticle(
           mcParticle, [](const edm4hep::MCParticle& p) {
             ActsFatras::Barcode result;
-            // TODO dont use podio internal id
-            result.setParticle(p.id());
+            result.setParticle(EDM4hepUtil::podioObjectIDToInteger(p.id()));
             return result;
           });
       unordered.push_back(particle);
@@ -78,8 +77,8 @@ ProcessCode EDM4hepSimHitReader::read(const AlgorithmContext& ctx) {
           simTrackerHit,
           [](const edm4hep::MCParticle& particle) {
             ActsFatras::Barcode result;
-            // TODO dont use podio internal id
-            result.setParticle(particle.id());
+            result.setParticle(
+                EDM4hepUtil::podioObjectIDToInteger(particle.id()));
             return result;
           },
           [&](std::uint64_t cellId) {

Examples/Io/EDM4hep/src/EDM4hepUtil.cpp

@@ -152,7 +152,8 @@ Measurement EDM4hepUtil::readMeasurement(
   // no need for digitization as we only want to identify the sensor
   Acts::GeometryIdentifier geometryId = geometryMapper(from.getCellID());
 
-  IndexSourceLink sourceLink{geometryId, from.id()};
+  IndexSourceLink sourceLink{
+      geometryId, static_cast<Index>(podioObjectIDToInteger(from.id()))};
 
   auto pos = from.getPosition();
   auto cov = from.getCovMatrix();

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,
         )
     )