refactor: Disentangle stepper state creation and initialization (#4069)

Disentangles the state creation and initialization for the steppers.
This is analog to what is done in the navigator. The idea is that a
state can be used for multiple propagations so its creation should be
decoupled from the initialization to allow the reuse of the existing
state. At the same time `resetState` is removed as it has no use
anymore. `initialize` will reset and reinitialize the state.

Pulled out of #4036

<!-- This is an auto-generated comment: release notes by coderabbit.ai
-->

## Summary by CodeRabbit

- **Refactor**
- Redesigned the state management workflow across propagation and
tracking features to separate object creation from parameter
initialization, resulting in a clearer, more streamlined API.
- Removed legacy reset mechanisms in favor of explicit initialization
methods that accommodate configurable parameters.
  
- **Tests**
- Updated unit tests to align with the new initialization process,
ensuring improved consistency and robustness across propagation-related
functionalities.

<!-- end of auto-generated comment: release notes by coderabbit.ai -->

Core/include/Acts/Propagator/AtlasStepper.hpp

@@ -12,6 +12,8 @@
 #include "Acts/Utilities/detail/ReferenceWrapperAnyCompat.hpp"
 
 #include "Acts/Definitions/Algebra.hpp"
+#include "Acts/Definitions/Common.hpp"
+#include "Acts/Definitions/TrackParametrization.hpp"
 #include "Acts/Definitions/Units.hpp"
 #include "Acts/EventData/TrackParameters.hpp"
 #include "Acts/EventData/TransformationHelpers.hpp"
@@ -137,33 +139,53 @@ class AtlasStepper {
 
   explicit AtlasStepper(const Config& config) : m_bField(config.bField) {}
 
-  State makeState(const Options& options,
-                  const BoundTrackParameters& par) const {
+  State makeState(const Options& options) const {
     State state{options, m_bField->makeCache(options.magFieldContext)};
+    return state;
+  }
+
+  void initialize(State& state, const BoundTrackParameters& par) const {
+    initialize(state, par.parameters(), par.covariance(),
+               par.particleHypothesis(), par.referenceSurface());
+  }
+
+  void initialize(State& state, const BoundVector& boundParams,
+                  const std::optional<BoundMatrix>& cov,
+                  ParticleHypothesis particleHypothesis,
+                  const Surface& surface) const {
+    state.particleHypothesis = particleHypothesis;
 
-    state.particleHypothesis = par.particleHypothesis();
+    state.pathAccumulated = 0;
+    state.nSteps = 0;
+    state.nStepTrials = 0;
+    state.stepSize = ConstrainedStep(state.options.maxStepSize);
+    state.previousStepSize = 0;
+    state.statistics = StepperStatistics();
 
     // The rest of this constructor is copy&paste of AtlasStepper::update() -
     // this is a nasty but working solution for the stepper state without
     // functions
 
-    const auto pos = par.position(options.geoContext);
-    const auto Vp = par.parameters();
+    const auto& Vp = boundParams;
 
     double Sf = std::sin(Vp[eBoundPhi]);
     double Cf = std::cos(Vp[eBoundPhi]);
     double Se = std::sin(Vp[eBoundTheta]);
     double Ce = std::cos(Vp[eBoundTheta]);
 
+    const Vector3 dir = {Cf * Se, Sf * Se, Ce};
+    const auto pos = surface.localToGlobal(
+        state.options.geoContext, boundParams.segment<2>(eBoundLoc0), dir);
+
     double* pVector = state.pVector.data();
 
     pVector[0] = pos[ePos0];
     pVector[1] = pos[ePos1];
     pVector[2] = pos[ePos2];
-    pVector[3] = par.time();
-    pVector[4] = Cf * Se;
-    pVector[5] = Sf * Se;
-    pVector[6] = Ce;
+    pVector[3] = boundParams[eBoundTime];
+    pVector[4] = dir[ePos0];
+    pVector[5] = dir[ePos1];
+    pVector[6] = dir[ePos2];
     pVector[7] = Vp[eBoundQOverP];
 
     // @todo: remove magic numbers - is that the charge ?
@@ -173,13 +195,13 @@ class AtlasStepper {
     }
 
     // prepare the jacobian if we have a covariance
-    if (par.covariance()) {
+    state.covTransport = cov.has_value();
+    if (state.covTransport) {
       // copy the covariance matrix
-      state.covariance = new BoundSquareMatrix(*par.covariance());
-      state.covTransport = true;
+      state.covariance = new BoundSquareMatrix(*cov);
       state.useJacobian = true;
-      const auto transform = par.referenceSurface().referenceFrame(
-          options.geoContext, pos, par.direction());
+      const auto transform =
+          surface.referenceFrame(state.options.geoContext, pos, dir);
 
       pVector[8] = transform(0, eBoundLoc0);
       pVector[16] = transform(0, eBoundLoc1);
@@ -243,7 +265,6 @@ class AtlasStepper {
       pVector[59] = 0.;
 
       // special treatment for surface types
-      const auto& surface = par.referenceSurface();
       // the disc needs polar coordinate adaptations
       if (surface.type() == Surface::Disc) {
         double lCf = std::cos(Vp[1]);
@@ -305,34 +326,8 @@ class AtlasStepper {
       }
     }
 
-    state.stepSize = ConstrainedStep(options.maxStepSize);
-
     // now declare the state as ready
     state.state_ready = true;
-
-    return state;
-  }
-
-  /// @brief Resets the state
-  ///
-  /// @param [in, out] state State of the stepper
-  /// @param [in] boundParams Parameters in bound parametrisation
-  /// @param [in] cov Covariance matrix
-  /// @param [in] surface Reset state will be on this surface
-  /// @param [in] stepSize Step size
-  void resetState(
-      State& state, const BoundVector& boundParams,
-      const BoundSquareMatrix& cov, const Surface& surface,
-      const double stepSize = std::numeric_limits<double>::max()) const {
-    // Update the stepping state
-    update(state,
-           transformBoundToFreeParameters(surface, state.options.geoContext,
-                                          boundParams),
-           boundParams, cov, surface);
-    state.stepSize = ConstrainedStep(stepSize);
-    state.pathAccumulated = 0.;
-
-    setIdentityJacobian(state);
   }
 
   /// Get the field for the stepping

Core/include/Acts/Propagator/EigenStepper.hpp

@@ -25,7 +25,6 @@
 #include "Acts/Utilities/Intersection.hpp"
 #include "Acts/Utilities/Result.hpp"
 
-#include <limits>
 #include <type_traits>
 
 namespace Acts {
@@ -150,20 +149,14 @@ class EigenStepper {
   /// @param [in] config The configuration of the stepper
   explicit EigenStepper(const Config& config) : m_bField(config.bField) {}
 
-  State makeState(const Options& options,
-                  const BoundTrackParameters& par) const;
+  State makeState(const Options& options) const;
 
-  /// @brief Resets the state
-  ///
-  /// @param [in, out] state State of the stepper
-  /// @param [in] boundParams Parameters in bound parametrisation
-  /// @param [in] cov Covariance matrix
-  /// @param [in] surface The reference surface of the bound parameters
-  /// @param [in] stepSize Step size
-  void resetState(
-      State& state, const BoundVector& boundParams,
-      const BoundSquareMatrix& cov, const Surface& surface,
-      const double stepSize = std::numeric_limits<double>::max()) const;
+  void initialize(State& state, const BoundTrackParameters& par) const;
+
+  void initialize(State& state, const BoundVector& boundParams,
+                  const std::optional<BoundMatrix>& cov,
+                  ParticleHypothesis particleHypothesis,
+                  const Surface& surface) const;
 
   /// Get the field for the stepping, it checks first if the access is still
   /// within the Cell, and updates the cell if necessary.

Core/include/Acts/Propagator/EigenStepper.ipp

@@ -20,57 +20,49 @@ Acts::EigenStepper<E>::EigenStepper(
     : m_bField(std::move(bField)) {}
 
 template <typename E>
-auto Acts::EigenStepper<E>::makeState(
-    const Options& options, const BoundTrackParameters& par) const -> State {
+auto Acts::EigenStepper<E>::makeState(const Options& options) const -> State {
   State state{options, m_bField->makeCache(options.magFieldContext)};
-
-  state.particleHypothesis = par.particleHypothesis();
-
-  Vector3 position = par.position(options.geoContext);
-  Vector3 direction = par.direction();
-  state.pars.template segment<3>(eFreePos0) = position;
-  state.pars.template segment<3>(eFreeDir0) = direction;
-  state.pars[eFreeTime] = par.time();
-  state.pars[eFreeQOverP] = par.parameters()[eBoundQOverP];
-
-  // Init the jacobian matrix if needed
-  if (par.covariance()) {
-    // Get the reference surface for navigation
-    const auto& surface = par.referenceSurface();
-    // set the covariance transport flag to true and copy
-    state.covTransport = true;
-    state.cov = BoundSquareMatrix(*par.covariance());
-    state.jacToGlobal =
-        surface.boundToFreeJacobian(options.geoContext, position, direction);
-  }
-
-  state.stepSize = ConstrainedStep(options.maxStepSize);
-
   return state;
 }
 
 template <typename E>
-void Acts::EigenStepper<E>::resetState(State& state,
+void Acts::EigenStepper<E>::initialize(State& state,
+                                       const BoundTrackParameters& par) const {
+  initialize(state, par.parameters(), par.covariance(),
+             par.particleHypothesis(), par.referenceSurface());
+}
+
+template <typename E>
+void Acts::EigenStepper<E>::initialize(State& state,
                                        const BoundVector& boundParams,
-                                       const BoundSquareMatrix& cov,
-                                       const Surface& surface,
-                                       const double stepSize) const {
+                                       const std::optional<BoundMatrix>& cov,
+                                       ParticleHypothesis particleHypothesis,
+                                       const Surface& surface) const {
   FreeVector freeParams = transformBoundToFreeParameters(
       surface, state.options.geoContext, boundParams);
 
-  // Update the stepping state
+  state.particleHypothesis = particleHypothesis;
+
+  state.pathAccumulated = 0;
+  state.nSteps = 0;
+  state.nStepTrials = 0;
+  state.stepSize = ConstrainedStep(state.options.maxStepSize);
+  state.previousStepSize = 0;
+  state.statistics = StepperStatistics();
+
   state.pars = freeParams;
-  state.cov = cov;
-  state.stepSize = ConstrainedStep(stepSize);
-  state.pathAccumulated = 0.;
 
-  // Reinitialize the stepping jacobian
-  state.jacToGlobal = surface.boundToFreeJacobian(
-      state.options.geoContext, freeParams.template segment<3>(eFreePos0),
-      freeParams.template segment<3>(eFreeDir0));
-  state.jacobian = BoundMatrix::Identity();
-  state.jacTransport = FreeMatrix::Identity();
-  state.derivative = FreeVector::Zero();
+  // Init the jacobian matrix if needed
+  state.covTransport = cov.has_value();
+  if (state.covTransport) {
+    state.cov = *cov;
+    state.jacToGlobal = surface.boundToFreeJacobian(
+        state.options.geoContext, freeParams.segment<3>(eFreePos0),
+        freeParams.segment<3>(eFreeDir0));
+    state.jacobian = BoundMatrix::Identity();
+    state.jacTransport = FreeMatrix::Identity();
+    state.derivative = FreeVector::Zero();
+  }
 }
 
 template <typename E>

Core/include/Acts/Propagator/MultiEigenStepperLoop.hpp

@@ -198,6 +198,9 @@ class MultiEigenStepperLoop : public EigenStepper<extension_t> {
       SingleState state;
       double weight;
       IntersectionStatus status;
+
+      Component(SingleState state_, double weight_, IntersectionStatus status_)
+          : state(std::move(state_)), weight(weight_), status(status_) {}
     };
 
     Options options;
@@ -240,49 +243,34 @@ class MultiEigenStepperLoop : public EigenStepper<extension_t> {
                             getDefaultLogger("GSF", Logging::INFO))
       : EigenStepper<extension_t>(config), m_logger(std::move(logger)) {}
 
-  /// Construct and initialize a state
-  State makeState(const Options& options,
+  State makeState(const Options& options) const {
+    State state(options);
+    return state;
+  }
+
+  void initialize(State& state,
                   const MultiComponentBoundTrackParameters& par) const {
     if (par.components().empty()) {
       throw std::invalid_argument(
           "Cannot construct MultiEigenStepperLoop::State with empty "
           "multi-component parameters");
     }
 
-    State state(options);
-
     state.particleHypothesis = par.particleHypothesis();
 
     const auto surface = par.referenceSurface().getSharedPtr();
 
     for (auto i = 0ul; i < par.components().size(); ++i) {
       const auto& [weight, singlePars] = par[i];
-      state.components.push_back({SingleStepper::makeState(options, singlePars),
-                                  weight, IntersectionStatus::onSurface});
+      auto& cmp =
+          state.components.emplace_back(SingleStepper::makeState(state.options),
+                                        weight, IntersectionStatus::onSurface);
+      SingleStepper::initialize(cmp.state, singlePars);
     }
 
     if (std::get<2>(par.components().front())) {
       state.covTransport = true;
     }
-
-    return state;
-  }
-
-  /// @brief Resets the state
-  ///
-  /// @param [in, out] state State of the stepper
-  /// @param [in] boundParams Parameters in bound parametrisation
-  /// @param [in] cov Covariance matrix
-  /// @param [in] surface The reference surface of the bound parameters
-  /// @param [in] stepSize Step size
-  void resetState(
-      State& state, const BoundVector& boundParams,
-      const BoundSquareMatrix& cov, const Surface& surface,
-      const double stepSize = std::numeric_limits<double>::max()) const {
-    for (auto& component : state.components) {
-      SingleStepper::resetState(component.state, boundParams, cov, surface,
-                                stepSize);
-    }
   }
 
   /// A proxy struct which allows access to a single component of the
@@ -420,8 +408,10 @@ class MultiEigenStepperLoop : public EigenStepper<extension_t> {
   Result<ComponentProxy> addComponent(State& state,
                                       const BoundTrackParameters& pars,
                                       double weight) const {
-    state.components.push_back({SingleStepper::makeState(state.options, pars),
-                                weight, IntersectionStatus::onSurface});
+    auto& cmp =
+        state.components.emplace_back(SingleStepper::makeState(state.options),
+                                      weight, IntersectionStatus::onSurface);
+    SingleStepper::initialize(cmp.state, pars);
 
     return ComponentProxy{state.components.back(), state};
   }

Core/include/Acts/Propagator/Propagator.hpp

@@ -386,8 +386,9 @@ class Propagator final
  private:
   const Logger& logger() const { return *m_logger; }
 
-  template <typename propagator_state_t, typename path_aborter_t>
-  void initialize(propagator_state_t& state) const;
+  template <typename propagator_state_t, typename parameters_t,
+            typename path_aborter_t>
+  void initialize(propagator_state_t& state, const parameters_t& start) const;
 
   template <typename propagator_state_t, typename result_t>
   void moveStateToResult(propagator_state_t& state, result_t& result) const;