fix: Multi-Stepper stepping error in multiple components + refactoring [backport #1339 to develop/v19.x]

Core/include/Acts/Propagator/MultiEigenStepperLoop.hpp

@@ -177,9 +177,6 @@ template <typename extensionlist_t = StepperExtensionList<DefaultExtension>,
           typename auctioneer_t = detail::VoidAuctioneer>
 class MultiEigenStepperLoop
     : public EigenStepper<extensionlist_t, auctioneer_t> {
-  /// Allows logging in the member functions
-  const LoggerWrapper logger;
-
   /// Limits the number of steps after at least one component reached the
   /// surface
   std::size_t m_stepLimitAfterFirstComponentOnSurface = 50;
@@ -209,9 +206,8 @@ class MultiEigenStepperLoop
   static constexpr int maxComponents = std::numeric_limits<int>::max();
 
   /// Constructor from a magnetic field and a optionally provided Logger
-  MultiEigenStepperLoop(std::shared_ptr<const MagneticFieldProvider> bField,
-                        LoggerWrapper l = getDummyLogger())
-      : EigenStepper<extensionlist_t, auctioneer_t>(bField), logger(l) {}
+  MultiEigenStepperLoop(std::shared_ptr<const MagneticFieldProvider> bField)
+      : EigenStepper<extensionlist_t, auctioneer_t>(bField) {}
 
   struct State {
     /// The struct that stores the individual components
@@ -608,9 +604,10 @@ class MultiEigenStepperLoop
   /// @param state [in,out] The stepping state (thread-local cache)
   /// @param surface [in] The surface provided
   /// @param bcheck [in] The boundary check for this status update
+  /// @param logger [in] A @c LoggerWrapper instance
   Intersection3D::Status updateSurfaceStatus(
       State& state, const Surface& surface, const BoundaryCheck& bcheck,
-      LoggerWrapper /*extLogger*/ = getDummyLogger()) const {
+      LoggerWrapper logger = getDummyLogger()) const {
     using Status = Intersection3D::Status;
 
     std::array<int, 4> counts = {0, 0, 0, 0};

Core/include/Acts/Propagator/MultiEigenStepperLoop.ipp

@@ -35,16 +35,14 @@ auto MultiEigenStepperLoop<E, R, A>::boundState(
       }
     }
 
-    if (failedBoundTransforms > 0) {
-      ACTS_ERROR("Multi component bound state: "
-                 << failedBoundTransforms << " of " << numberComponents(state)
-                 << " transforms failed");
-    }
-
     if (states.size() == 0) {
       return MultiStepperError::AllComponentsConversionToBoundFailed;
     }
 
+    if (failedBoundTransforms > 0) {
+      return MultiStepperError::SomeComponentsConversionToBoundFailed;
+    }
+
     // TODO At ATLAS, the final parameters seem to be computed with the mode of
     // the mixture. At the moment, we use the mean of the mixture here, but
     // there should be done a comparison sometimes in the future. This could
@@ -100,8 +98,29 @@ template <typename E, typename R, typename A>
 template <typename propagator_state_t>
 Result<double> MultiEigenStepperLoop<E, R, A>::step(
     propagator_state_t& state) const {
+  const auto& logger = state.options.logger;
   State& stepping = state.stepping;
 
+  // It is not possible to remove components from the vector, since the
+  // GSF actor relies on the fact that the ordering and number of
+  // components does not change
+  auto invalidateComponent = [](auto& cmp) {
+    cmp.status = Intersection3D::Status::missed;
+    cmp.weight = 0.0;
+    cmp.state.pars.template segment<3>(eFreeDir0) = Vector3::Zero();
+  };
+
+  // Lambda for reweighting the components
+  auto reweight = [](auto& cmps) {
+    ActsScalar sumOfWeights = 0.0;
+    for (const auto& cmp : cmps) {
+      sumOfWeights += cmp.weight;
+    }
+    for (auto& cmp : cmps) {
+      cmp.weight /= sumOfWeights;
+    }
+  };
+
   // Update step count
   stepping.steps++;
 
@@ -113,30 +132,19 @@ Result<double> MultiEigenStepperLoop<E, R, A>::step(
     // surface, reweight the components, perform no step and return 0
     if (*stepping.stepCounterAfterFirstComponentOnSurface >=
         m_stepLimitAfterFirstComponentOnSurface) {
-      auto& cmps = stepping.components;
-
-      // It is not possible to remove components from the vector, since the
-      // GSF actor relies on the fact that the ordering and number of
-      // components does not change
-      for (auto& cmp : cmps) {
+      for (auto& cmp : stepping.components) {
         if (cmp.status != Intersection3D::Status::onSurface) {
-          cmp.status = Intersection3D::Status::missed;
-          cmp.weight = 0.0;
-          cmp.state.pars.template segment<3>(eFreeDir0) = Vector3::Zero();
+          invalidateComponent(cmp);
         }
       }
 
-      // Reweight
-      const auto sum_of_weights = std::accumulate(
-          cmps.begin(), cmps.end(), ActsScalar{0},
-          [](auto sum, const auto& cmp) { return sum + cmp.weight; });
-      for (auto& cmp : cmps) {
-        cmp.weight /= sum_of_weights;
-      }
+      reweight(stepping.components);
 
+      ACTS_VERBOSE("Stepper performed "
+                   << m_stepLimitAfterFirstComponentOnSurface
+                   << " after the first component hit a surface.");
       ACTS_VERBOSE(
-          "hit m_stepLimitAfterFirstComponentOnSurface, "
-          "perform no step");
+          "-> remove all components not on a surface, perform no step");
 
       stepping.stepCounterAfterFirstComponentOnSurface.reset();
 
@@ -146,15 +154,17 @@ Result<double> MultiEigenStepperLoop<E, R, A>::step(
 
   // Loop over all components and collect results in vector, write some
   // summary information to a stringstream
-  SmallVector<Result<double>> results;
-  std::stringstream ss;
+  SmallVector<std::optional<Result<double>>> results;
   double accumulatedPathLength = 0.0;
+  std::size_t errorSteps = 0;
 
   for (auto& component : stepping.components) {
     // We must also propagate missed components for the case that all
     // components miss the target and we need to re-target
     if (component.status == Intersection3D::Status::onSurface) {
-      ss << "cmp skipped\t";
+      // We need to add these, so the propagation does not fail if we have only
+      // components on surfaces and failing states
+      results.emplace_back(std::nullopt);
       continue;
     }
 
@@ -165,41 +175,50 @@ Result<double> MultiEigenStepperLoop<E, R, A>::step(
     ThisSinglePropState single_state(component.state, state.navigation,
                                      state.options, state.geoContext);
 
-    results.push_back(SingleStepper::step(single_state));
+    results.emplace_back(SingleStepper::step(single_state));
 
-    if (results.back().ok()) {
-      accumulatedPathLength += component.weight * *results.back();
-      ss << *results.back() << "\t";
+    if (results.back()->ok()) {
+      accumulatedPathLength += component.weight * results.back()->value();
     } else {
-      ss << "step error: " << results.back().error() << "\t";
+      ++errorSteps;
+      invalidateComponent(component);
     }
   }
 
-  // Return no component was updated
-  if (results.empty()) {
-    return 0.0;
+  // Since we have invalidated some components, we need to reweight
+  if (errorSteps > 0) {
+    reweight(stepping.components);
   }
 
-  // Collect pointers to results which are ok, since Result is not copyable
-  SmallVector<Result<double>*> ok_results;
-  for (auto& res : results) {
-    if (res.ok()) {
-      ok_results.push_back(&res);
+  // Print the result vector to a string so we can log it
+  auto summary = [](auto& result_vec) {
+    std::stringstream ss;
+    for (auto& optRes : result_vec) {
+      if (not optRes) {
+        ss << "on surface | ";
+      } else if (optRes->ok()) {
+        ss << optRes->value() << " | ";
+      } else {
+        ss << optRes->error() << " | ";
+      }
     }
+    auto str = ss.str();
+    str.resize(str.size() - 3);
+    return str;
+  };
+
+  // Print the summary
+  if (errorSteps == 0) {
+    ACTS_VERBOSE("Performed steps: " << summary(results));
+  } else {
+    ACTS_WARNING("Performed steps with errors: " << summary(results));
   }
 
   // Return error if there is no ok result
-  if (ok_results.empty()) {
+  if (errorSteps == results.size()) {
     return MultiStepperError::AllComponentsSteppingError;
   }
 
-  // Print the summary
-  if (ok_results.size() == results.size()) {
-    ACTS_VERBOSE("Performed steps: " << ss.str());
-  } else {
-    ACTS_WARNING("Performed steps with errors: " << ss.str());
-  }
-
   // Return the weighted accumulated path length of all successful steps
   stepping.pathAccumulated += accumulatedPathLength;
   return accumulatedPathLength;

Core/include/Acts/Propagator/MultiStepperError.hpp

@@ -18,7 +18,8 @@ enum class MultiStepperError {
   StateOfMultipleComponentsRequested = 2,
   AverageTrackLeftCurrentVolume = 3,
   AllComponentsSteppingError = 4,
-  AllComponentsConversionToBoundFailed = 5
+  AllComponentsConversionToBoundFailed = 5,
+  SomeComponentsConversionToBoundFailed = 6
 };
 
 std::error_code make_error_code(Acts::MultiStepperError e);

Core/src/Propagator/MultiStepperError.cpp

@@ -31,6 +31,8 @@ class MultiStepperErrorCategory : public std::error_category {
         return "Stepping error occured in all components";
       case MultiStepperError::AllComponentsConversionToBoundFailed:
         return "The conversion to the bound state failed for all components";
+      case MultiStepperError::SomeComponentsConversionToBoundFailed:
+        return "The conversion to the bound state failed for some components";
       default:
         return "unknown";
     }

Tests/UnitTests/Core/Propagator/MultiStepperTests.cpp

@@ -12,7 +12,7 @@
 #include "Acts/MagneticField/NullBField.hpp"
 #include "Acts/Propagator/MultiEigenStepperLoop.hpp"
 #include "Acts/Propagator/MultiStepperAborters.hpp"
-#include <Acts/Propagator/Navigator.hpp>
+#include "Acts/Propagator/Navigator.hpp"
 
 using namespace Acts;
 using namespace Acts::VectorHelpers;
@@ -40,6 +40,7 @@ struct Options {
   double stepSizeCutOff = 0.0;
   std::size_t maxRungeKuttaStepTrials = 10;
   double mass = 1.0;
+  LoggerWrapper logger = Acts::getDummyLogger();
 };
 
 struct Navigation {};