refactor!: Revise stepper state creation

Core/include/Acts/Propagator/EigenStepper.hpp

@@ -15,8 +15,6 @@
 #include "Acts/Definitions/Tolerance.hpp"
 #include "Acts/EventData/TrackParameters.hpp"
 #include "Acts/EventData/detail/CorrectedTransformationFreeToBound.hpp"
-#include "Acts/Geometry/GeometryContext.hpp"
-#include "Acts/MagneticField/MagneticFieldContext.hpp"
 #include "Acts/MagneticField/MagneticFieldProvider.hpp"
 #include "Acts/Propagator/ConstrainedStep.hpp"
 #include "Acts/Propagator/EigenStepperDefaultExtension.hpp"
@@ -27,7 +25,6 @@
 #include "Acts/Utilities/Intersection.hpp"
 #include "Acts/Utilities/Result.hpp"
 
-#include <functional>
 #include <limits>
 #include <type_traits>
 
@@ -60,6 +57,9 @@ class EigenStepper {
   };
 
   struct Options : public StepperPlainOptions {
+    Options(const GeometryContext& gctx, const MagneticFieldContext& mctx)
+        : StepperPlainOptions(gctx, mctx) {}
+
     void setPlainOptions(const StepperPlainOptions& options) {
       static_cast<StepperPlainOptions&>(*this) = options;
     }
@@ -70,41 +70,16 @@ class EigenStepper {
   /// It contains the stepping information and is provided thread local
   /// by the propagator
   struct State {
-    State() = delete;
-
     /// Constructor from the initial bound track parameters
     ///
-    /// @param [in] gctx is the context object for the geometry
+    /// @param [in] optionsIn The stepper options
     /// @param [in] fieldCacheIn is the cache object for the magnetic field
-    /// @param [in] par The track parameters at start
-    /// @param [in] ssize is the maximum step size
     ///
     /// @note the covariance matrix is copied when needed
-    explicit State(const GeometryContext& gctx,
-                   MagneticFieldProvider::Cache fieldCacheIn,
-                   const BoundTrackParameters& par,
-                   double ssize = std::numeric_limits<double>::max())
-        : particleHypothesis(par.particleHypothesis()),
-          stepSize(ssize),
-          fieldCache(std::move(fieldCacheIn)),
-          geoContext(gctx) {
-      Vector3 position = par.position(gctx);
-      Vector3 direction = par.direction();
-      pars.template segment<3>(eFreePos0) = position;
-      pars.template segment<3>(eFreeDir0) = direction;
-      pars[eFreeTime] = par.time();
-      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
-        covTransport = true;
-        cov = BoundSquareMatrix(*par.covariance());
-        jacToGlobal = surface.boundToFreeJacobian(gctx, position, direction);
-      }
-    }
+    State(const Options& optionsIn, MagneticFieldProvider::Cache fieldCacheIn)
+        : options(optionsIn), fieldCache(std::move(fieldCacheIn)) {}
+
+    Options options;
 
     /// Internal free vector parameters
     FreeVector pars = FreeVector::Zero();
@@ -149,9 +124,6 @@ class EigenStepper {
     /// See step() code for details.
     MagneticFieldProvider::Cache fieldCache;
 
-    /// The geometry context
-    std::reference_wrapper<const GeometryContext> geoContext;
-
     /// Algorithmic extension
     extension_t extension;
 
@@ -178,10 +150,8 @@ class EigenStepper {
   /// @param [in] config The configuration of the stepper
   explicit EigenStepper(const Config& config) : m_bField(config.bField) {}
 
-  State makeState(std::reference_wrapper<const GeometryContext> gctx,
-                  std::reference_wrapper<const MagneticFieldContext> mctx,
-                  const BoundTrackParameters& par,
-                  double ssize = std::numeric_limits<double>::max()) const;
+  State makeState(const Options& options,
+                  const BoundTrackParameters& par) const;
 
   /// @brief Resets the state
   ///

Core/include/Acts/Propagator/EigenStepper.ipp

@@ -21,10 +21,32 @@ Acts::EigenStepper<E>::EigenStepper(
 
 template <typename E>
 auto Acts::EigenStepper<E>::makeState(
-    std::reference_wrapper<const GeometryContext> gctx,
-    std::reference_wrapper<const MagneticFieldContext> mctx,
-    const BoundTrackParameters& par, double ssize) const -> State {
-  return State{gctx, m_bField->makeCache(mctx), par, ssize};
+    const Options& options, const BoundTrackParameters& par) 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>
@@ -33,8 +55,8 @@ void Acts::EigenStepper<E>::resetState(State& state,
                                        const BoundSquareMatrix& cov,
                                        const Surface& surface,
                                        const double stepSize) const {
-  FreeVector freeParams =
-      transformBoundToFreeParameters(surface, state.geoContext, boundParams);
+  FreeVector freeParams = transformBoundToFreeParameters(
+      surface, state.options.geoContext, boundParams);
 
   // Update the stepping state
   state.pars = freeParams;
@@ -44,7 +66,7 @@ void Acts::EigenStepper<E>::resetState(State& state,
 
   // Reinitialize the stepping jacobian
   state.jacToGlobal = surface.boundToFreeJacobian(
-      state.geoContext, freeParams.template segment<3>(eFreePos0),
+      state.options.geoContext, freeParams.template segment<3>(eFreePos0),
       freeParams.template segment<3>(eFreeDir0));
   state.jacobian = BoundMatrix::Identity();
   state.jacTransport = FreeMatrix::Identity();
@@ -57,10 +79,10 @@ auto Acts::EigenStepper<E>::boundState(
     const FreeToBoundCorrection& freeToBoundCorrection) const
     -> Result<BoundState> {
   return detail::boundState(
-      state.geoContext, surface, state.cov, state.jacobian, state.jacTransport,
-      state.derivative, state.jacToGlobal, state.pars, state.particleHypothesis,
-      state.covTransport && transportCov, state.pathAccumulated,
-      freeToBoundCorrection);
+      state.options.geoContext, surface, state.cov, state.jacobian,
+      state.jacTransport, state.derivative, state.jacToGlobal, state.pars,
+      state.particleHypothesis, state.covTransport && transportCov,
+      state.pathAccumulated, freeToBoundCorrection);
 }
 
 template <typename E>
@@ -113,7 +135,7 @@ void Acts::EigenStepper<E>::update(State& state, const FreeVector& freeParams,
   state.pars = freeParams;
   state.cov = covariance;
   state.jacToGlobal = surface.boundToFreeJacobian(
-      state.geoContext, freeParams.template segment<3>(eFreePos0),
+      state.options.geoContext, freeParams.template segment<3>(eFreePos0),
       freeParams.template segment<3>(eFreeDir0));
 }
 
@@ -139,10 +161,10 @@ template <typename E>
 void Acts::EigenStepper<E>::transportCovarianceToBound(
     State& state, const Surface& surface,
     const FreeToBoundCorrection& freeToBoundCorrection) const {
-  detail::transportCovarianceToBound(state.geoContext.get(), surface, state.cov,
-                                     state.jacobian, state.jacTransport,
-                                     state.derivative, state.jacToGlobal,
-                                     state.pars, freeToBoundCorrection);
+  detail::transportCovarianceToBound(
+      state.options.geoContext, surface, state.cov, state.jacobian,
+      state.jacTransport, state.derivative, state.jacToGlobal, state.pars,
+      freeToBoundCorrection);
 }
 
 template <typename E>

Core/include/Acts/Propagator/MultiEigenStepperLoop.hpp

@@ -30,7 +30,6 @@
 #include <algorithm>
 #include <cmath>
 #include <cstddef>
-#include <functional>
 #include <limits>
 #include <sstream>
 #include <vector>
@@ -156,6 +155,9 @@ class MultiEigenStepperLoop : public EigenStepper<extension_t> {
   /// @brief Typedef to the Single-Component Eigen Stepper
   using SingleStepper = EigenStepper<extension_t>;
 
+  /// @brief Typedef to the Single-Component Stepper Options
+  using SingleOptions = typename SingleStepper::Options;
+
   /// @brief Typedef to the State of the single component Stepper
   using SingleState = typename SingleStepper::State;
 
@@ -181,7 +183,10 @@ class MultiEigenStepperLoop : public EigenStepper<extension_t> {
     std::shared_ptr<const MagneticFieldProvider> bField;
   };
 
-  struct Options : public StepperPlainOptions {
+  struct Options : public SingleOptions {
+    Options(const GeometryContext& gctx, const MagneticFieldContext& mctx)
+        : SingleOptions(gctx, mctx) {}
+
     void setPlainOptions(const StepperPlainOptions& options) {
       static_cast<StepperPlainOptions&>(*this) = options;
     }
@@ -195,6 +200,8 @@ class MultiEigenStepperLoop : public EigenStepper<extension_t> {
       IntersectionStatus status;
     };
 
+    Options options;
+
     /// Particle hypothesis
     ParticleHypothesis particleHypothesis = ParticleHypothesis::pion();
 
@@ -205,58 +212,19 @@ class MultiEigenStepperLoop : public EigenStepper<extension_t> {
     double pathAccumulated = 0.;
     std::size_t steps = 0;
 
-    /// geoContext
-    std::reference_wrapper<const GeometryContext> geoContext;
-
-    /// MagneticFieldContext
-    std::reference_wrapper<const MagneticFieldContext> magContext;
-
     /// Step-limit counter which limits the number of steps when one component
     /// reached a surface
     std::optional<std::size_t> stepCounterAfterFirstComponentOnSurface;
 
     /// The stepper statistics
     StepperStatistics statistics;
 
-    /// No default constructor is provided
-    State() = delete;
-
     /// Constructor from the initial bound track parameters
     ///
-    /// @param [in] gctx is the context object for the geometry
-    /// @param [in] mctx is the context object for the magnetic field
-    /// @param [in] bfield the shared magnetic filed provider
-    /// @param [in] multipars The track multi-component track-parameters at start
-    /// @param [in] ssize is the maximum step size
+    /// @param [in] optionsIn The options for the stepper
     ///
     /// @note the covariance matrix is copied when needed
-    explicit State(const GeometryContext& gctx,
-                   const MagneticFieldContext& mctx,
-                   const std::shared_ptr<const MagneticFieldProvider>& bfield,
-                   const MultiComponentBoundTrackParameters& multipars,
-                   double ssize = std::numeric_limits<double>::max())
-        : particleHypothesis(multipars.particleHypothesis()),
-          geoContext(gctx),
-          magContext(mctx) {
-      if (multipars.components().empty()) {
-        throw std::invalid_argument(
-            "Cannot construct MultiEigenStepperLoop::State with empty "
-            "multi-component parameters");
-      }
-
-      const auto surface = multipars.referenceSurface().getSharedPtr();
-
-      for (auto i = 0ul; i < multipars.components().size(); ++i) {
-        const auto& [weight, singlePars] = multipars[i];
-        components.push_back(
-            {SingleState(gctx, bfield->makeCache(mctx), singlePars, ssize),
-             weight, IntersectionStatus::onSurface});
-      }
-
-      if (std::get<2>(multipars.components().front())) {
-        covTransport = true;
-      }
-    }
+    explicit State(const Options& optionsIn) : options(optionsIn) {}
   };
 
   /// Constructor from a magnetic field and a optionally provided Logger
@@ -273,11 +241,31 @@ class MultiEigenStepperLoop : public EigenStepper<extension_t> {
       : EigenStepper<extension_t>(config), m_logger(std::move(logger)) {}
 
   /// Construct and initialize a state
-  State makeState(std::reference_wrapper<const GeometryContext> gctx,
-                  std::reference_wrapper<const MagneticFieldContext> mctx,
-                  const MultiComponentBoundTrackParameters& par,
-                  double ssize = std::numeric_limits<double>::max()) const {
-    return State(gctx, mctx, SingleStepper::m_bField, par, ssize);
+  State makeState(const Options& options,
+                  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});
+    }
+
+    if (std::get<2>(par.components().front())) {
+      state.covTransport = true;
+    }
+
+    return state;
   }
 
   /// @brief Resets the state
@@ -432,11 +420,8 @@ class MultiEigenStepperLoop : public EigenStepper<extension_t> {
   Result<ComponentProxy> addComponent(State& state,
                                       const BoundTrackParameters& pars,
                                       double weight) const {
-    state.components.push_back(
-        {SingleState(state.geoContext,
-                     SingleStepper::m_bField->makeCache(state.magContext),
-                     pars),
-         weight, IntersectionStatus::onSurface});
+    state.components.push_back({SingleStepper::makeState(state.options, pars),
+                                weight, IntersectionStatus::onSurface});
 
     return ComponentProxy{state.components.back(), state};
   }
@@ -541,8 +526,8 @@ class MultiEigenStepperLoop : public EigenStepper<extension_t> {
 
     ACTS_VERBOSE("Component status wrt "
                  << surface.geometryId() << " at {"
-                 << surface.center(state.geoContext).transpose() << "}:\t"
-                 << [&]() {
+                 << surface.center(state.options.geoContext).transpose()
+                 << "}:\t" << [&]() {
                       std::stringstream ss;
                       for (auto& component : state.components) {
                         ss << component.status << "\t";
@@ -599,7 +584,8 @@ class MultiEigenStepperLoop : public EigenStepper<extension_t> {
 
     for (auto& component : state.components) {
       auto intersection = surface.intersect(
-          component.state.geoContext, SingleStepper::position(component.state),
+          component.state.options.geoContext,
+          SingleStepper::position(component.state),
           direction * SingleStepper::direction(component.state),
           BoundaryTolerance::None())[oIntersection.index()];
 

Core/include/Acts/Propagator/MultiEigenStepperLoop.ipp

@@ -34,7 +34,7 @@ auto MultiEigenStepperLoop<E, R>::boundState(
     // onSurface.
     cmpState.pars.template segment<3>(eFreePos0) =
         surface
-            .intersect(state.geoContext,
+            .intersect(state.options.geoContext,
                        cmpState.pars.template segment<3>(eFreePos0),
                        cmpState.pars.template segment<3>(eFreeDir0),
                        BoundaryTolerance::Infinite())
@@ -78,7 +78,7 @@ auto MultiEigenStepperLoop<E, R>::curvilinearState(
         state.components[i].state, transportCov);
 
     cmps.emplace_back(state.components[i].weight,
-                      cp.fourPosition(state.geoContext), cp.direction(),
+                      cp.fourPosition(state.options.geoContext), cp.direction(),
                       cp.qOverP(),
                       cp.covariance().value_or(BoundSquareMatrix::Zero()));
     accumulatedPathLength += state.components[i].weight * pl;