refactor: Remove missed components as early as possible & more (#1921)

This refactors the `Acts::MultiEigenStepperLoop` code in a way, that missed components are removed as early as possible. Now the code should also be more consistent, so that some checks and workarounds in the `Acts::GsfActor` could be removed.

Core/include/Acts/Propagator/MultiEigenStepperLoop.hpp

@@ -286,7 +286,7 @@ class MultiEigenStepperLoop
         components.push_back(
             {SingleState(gctx, bfield->makeCache(mctx), std::move(singlePars),
                          ndir, ssize, stolerance),
-             weight, Intersection3D::Status::reachable});
+             weight, Intersection3D::Status::onSurface});
       }
 
       if (std::get<2>(multipars.components().front())) {
@@ -532,11 +532,15 @@ class MultiEigenStepperLoop
   }
 
   /// Get the number of components
+  ///
+  /// @param state [in,out] The stepping state (thread-local cache)
   std::size_t numberComponents(const State& state) const {
     return state.components.size();
   }
 
   /// Remove missed components from the component state
+  ///
+  /// @param state [in,out] The stepping state (thread-local cache)
   void removeMissedComponents(State& state) const {
     auto new_end = std::remove_if(
         state.components.begin(), state.components.end(), [](const auto& cmp) {
@@ -546,10 +550,31 @@ class MultiEigenStepperLoop
     state.components.erase(new_end, state.components.end());
   }
 
+  /// Reweight the components
+  ///
+  /// @param [in,out] state The stepping state (thread-local cache)
+  void reweightComponents(State& state) const {
+    ActsScalar sumOfWeights = 0.0;
+    for (const auto& cmp : state.components) {
+      sumOfWeights += cmp.weight;
+    }
+    for (auto& cmp : state.components) {
+      cmp.weight /= sumOfWeights;
+    }
+  }
+
   /// Reset the number of components
+  ///
+  /// @param [in,out] state  The stepping state (thread-local cache)
   void clearComponents(State& state) const { state.components.clear(); }
 
   /// Add a component to the Multistepper
+  ///
+  /// @param [in,out] state  The stepping state (thread-local cache)
+  /// @param [in] pars Parameters of the component to add
+  /// @param [in] weight Weight of the component to add
+  ///
+  /// @note: It is not ensured that the weights are normalized afterwards
   /// @note This function makes no garantuees about how new components are
   /// initialized, it is up to the caller to ensure that all components are
   /// valid in the end.
@@ -631,6 +656,14 @@ class MultiEigenStepperLoop
       ++counts[static_cast<std::size_t>(component.status)];
     }
 
+    // If at least one component is on a surface, we can remove all missed
+    // components before the step. If not, we must keep them for the case that
+    // all components miss and we need to retarget
+    if (counts[static_cast<std::size_t>(Status::onSurface)] > 0) {
+      removeMissedComponents(state);
+      reweightComponents(state);
+    }
+
     ACTS_VERBOSE("Component status wrt "
                  << surface.geometryId() << " at {"
                  << surface.center(state.geoContext).transpose() << "}:\t"

Core/include/Acts/Propagator/MultiEigenStepperLoop.ipp

@@ -129,22 +129,14 @@ template <typename E, typename R, typename A>
 template <typename propagator_state_t>
 Result<double> MultiEigenStepperLoop<E, R, A>::step(
     propagator_state_t& state) const {
+  using Status = Acts::Intersection3D::Status;
+
   State& stepping = state.stepping;
+  auto& components = stepping.components;
 
   // @TODO: This needs to be a real logger
   const Logger& logger = getDummyLogger();
 
-  // 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++;
 
@@ -156,14 +148,14 @@ Result<double> MultiEigenStepperLoop<E, R, A>::step(
     // surface, reweight the components, perform no step and return 0
     if (*stepping.stepCounterAfterFirstComponentOnSurface >=
         m_stepLimitAfterFirstComponentOnSurface) {
-      for (auto& cmp : stepping.components) {
-        if (cmp.status != Intersection3D::Status::onSurface) {
-          cmp.status = Intersection3D::Status::missed;
+      for (auto& cmp : components) {
+        if (cmp.status != Status::onSurface) {
+          cmp.status = Status::missed;
         }
       }
 
       removeMissedComponents(stepping);
-      reweight(stepping.components);
+      reweightComponents(stepping);
 
       ACTS_VERBOSE("Stepper performed "
                    << m_stepLimitAfterFirstComponentOnSurface
@@ -177,43 +169,66 @@ Result<double> MultiEigenStepperLoop<E, R, A>::step(
     }
   }
 
+  // Flag indicating if we need to reweight in the end
+  bool reweightNecessary = false;
+
+  // If at least one component is on a surface, we can remove all missed
+  // components before the step. If not, we must keep them for the case that all
+  // components miss and we need to retarget
+  const auto cmpsOnSurface =
+      std::count_if(components.cbegin(), components.cend(), [&](auto& cmp) {
+        return cmp.status == Intersection3D::Status::onSurface;
+      });
+
+  if (cmpsOnSurface > 0) {
+    removeMissedComponents(stepping);
+    reweightNecessary = true;
+  }
+
   // Loop over all components and collect results in vector, write some
   // summary information to a stringstream
   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) {
+  // Type of the proxy single propagation2 state
+  using ThisSinglePropState =
+      SinglePropState<SingleState, decltype(state.navigation),
+                      decltype(state.options), decltype(state.geoContext)>;
+
+  // Lambda that performs the step for a component and returns false if the step
+  // went ok and true if there was an error
+  auto errorInStep = [&](auto& component) {
+    if (component.status == Status::onSurface) {
       // 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;
+      return false;
     }
 
-    using ThisSinglePropState =
-        SinglePropState<SingleState, decltype(state.navigation),
-                        decltype(state.options), decltype(state.geoContext)>;
-
     ThisSinglePropState single_state(component.state, state.navigation,
                                      state.options, state.geoContext);
 
     results.emplace_back(SingleStepper::step(single_state));
 
     if (results.back()->ok()) {
       accumulatedPathLength += component.weight * results.back()->value();
+      return false;
     } else {
       ++errorSteps;
-      component.status = Intersection3D::Status::missed;
+      reweightNecessary = true;
+      return true;
     }
-  }
+  };
 
-  // Since we have invalidated some components, we need to reweight
-  if (errorSteps > 0) {
-    removeMissedComponents(stepping);
-    reweight(stepping.components);
+  // Loop over components and remove errorous components
+  stepping.components.erase(
+      std::remove_if(components.begin(), components.end(), errorInStep),
+      components.end());
+
+  // Reweight if necessary
+  if (reweightNecessary) {
+    reweightComponents(stepping);
   }
 
   // Print the result vector to a string so we can log it

Core/include/Acts/Propagator/MultiStepperAborters.hpp

@@ -63,6 +63,8 @@ struct MultiStepperSurfaceReached {
       if (!SurfaceReached{}(singleState, singleStepper, targetSurface,
                             logger)) {
         reached = false;
+      } else {
+        cmp.status() = Acts::Intersection3D::Status::onSurface;
       }
     }
 
@@ -78,6 +80,11 @@ struct MultiStepperSurfaceReached {
         ACTS_VERBOSE("Reached target in average mode");
         state.navigation.currentSurface = &targetSurface;
         state.navigation.targetReached = true;
+
+        for (auto cmp : stepper.componentIterable(state.stepping)) {
+          cmp.status() = Intersection3D::Status::onSurface;
+        }
+
         return true;
       }
     }