refactor: Remove direction from stepper (#2312)

The propagation direction is currently stored in multiple places which makes it tedious to change it consistently during propagation which can easily result in bugs.

In this PR I propose to remove the `navDir` from the stepper interface and replace it consistently with the propagation option for the direction.

This should not effect any physics performance as it is a code refactor of how we pass down the propagation direction to the stepper.

In a follow up PR I want to do the same thing for the navigation.

Core/include/Acts/Material/MaterialCollector.hpp

@@ -89,14 +89,14 @@ struct MaterialCollector {
             ACTS_VERBOSE("Update on start surface: post-update mode.");
             prepofu = navigator.currentSurface(state.navigation)
                           ->surfaceMaterial()
-                          ->factor(state.stepping.navDir,
+                          ->factor(state.options.direction,
                                    MaterialUpdateStage::PostUpdate);
           } else if (navigator.targetSurface(state.navigation) ==
                      navigator.currentSurface(state.navigation)) {
             ACTS_VERBOSE("Update on target surface: pre-update mode");
             prepofu = navigator.currentSurface(state.navigation)
                           ->surfaceMaterial()
-                          ->factor(state.stepping.navDir,
+                          ->factor(state.options.direction,
                                    MaterialUpdateStage::PreUpdate);
           } else {
             ACTS_VERBOSE("Update while pass through: full mode.");

Core/include/Acts/Navigation/DetectorNavigator.hpp

@@ -222,8 +222,9 @@ class DetectorNavigator {
                                   << surface.geometryId());
       // Estimate the surface status
       bool boundaryCheck = c.boundaryCheck;
-      auto surfaceStatus = stepper.updateSurfaceStatus(state.stepping, surface,
-                                                       boundaryCheck, logger());
+      auto surfaceStatus = stepper.updateSurfaceStatus(
+          state.stepping, surface, state.options.direction, boundaryCheck,
+          state.options.targetTolerance, logger());
       if (surfaceStatus == Intersection3D::Status::reachable) {
         ACTS_VERBOSE(volInfo(state)
                      << posInfo(state, stepper)
@@ -290,7 +291,8 @@ class DetectorNavigator {
 
     // TODO not sure about the boundary check
     auto surfaceStatus = stepper.updateSurfaceStatus(
-        state.stepping, *nextSurface, boundaryCheck, logger());
+        state.stepping, *nextSurface, state.options.direction, boundaryCheck,
+        state.options.targetTolerance, logger());
 
     // Check if we are at a surface
     if (surfaceStatus == Intersection3D::Status::onSurface) {

Core/include/Acts/Propagator/AtlasStepper.hpp

@@ -52,17 +52,14 @@ class AtlasStepper {
     /// @param [in] gctx The geometry context tof this call
     /// @param [in] fieldCacheIn The magnetic field cache for this call
     /// @param [in] pars Input parameters
-    /// @param [in] ndir The navigation direction w.r.t. parameters
     /// @param [in] ssize the steps size limitation
     /// @param [in] stolerance is the stepping tolerance
     template <typename Parameters>
     State(const GeometryContext& gctx,
           MagneticFieldProvider::Cache fieldCacheIn, const Parameters& pars,
-          Direction ndir = Direction::Forward,
           double ssize = std::numeric_limits<double>::max(),
           double stolerance = s_onSurfaceTolerance)
-        : navDir(ndir),
-          field(0., 0., 0.),
+        : field(0., 0., 0.),
           stepSize(ssize),
           tolerance(stolerance),
           fieldCache(std::move(fieldCacheIn)),
@@ -233,7 +230,6 @@ class AtlasStepper {
     // optimisation that init is not called twice
     bool state_ready = false;
     // configuration
-    Direction navDir = Direction::Forward;
     bool useJacobian = false;
     double step = 0;
     double maxPathLength = 0;
@@ -298,11 +294,9 @@ class AtlasStepper {
   State makeState(std::reference_wrapper<const GeometryContext> gctx,
                   std::reference_wrapper<const MagneticFieldContext> mctx,
                   const GenericBoundTrackParameters<charge_t>& par,
-                  Direction navDir = Direction::Forward,
                   double ssize = std::numeric_limits<double>::max(),
                   double stolerance = s_onSurfaceTolerance) const {
-    return State{gctx,      m_bField->makeCache(mctx), par, navDir, ssize,
-                 stolerance};
+    return State{gctx, m_bField->makeCache(mctx), par, ssize, stolerance};
   }
 
   /// @brief Resets the state
@@ -311,18 +305,16 @@ class AtlasStepper {
   /// @param [in] boundParams Parameters in bound parametrisation
   /// @param [in] cov Covariance matrix
   /// @param [in] surface Reset state will be on this surface
-  /// @param [in] navDir Navigation direction
   /// @param [in] stepSize Step size
   void resetState(
       State& state, const BoundVector& boundParams, const BoundSymMatrix& cov,
-      const Surface& surface, const Direction navDir = Direction::Forward,
+      const Surface& surface,
       const double stepSize = std::numeric_limits<double>::max()) const {
     // Update the stepping state
     update(state,
            detail::transformBoundToFreeParameters(surface, state.geoContext,
                                                   boundParams),
            boundParams, cov, surface);
-    state.navDir = navDir;
     state.stepSize = ConstrainedStep(stepSize);
     state.pathAccumulated = 0.;
 
@@ -383,15 +375,17 @@ class AtlasStepper {
   ///
   /// @param [in,out] state The stepping state (thread-local cache)
   /// @param [in] surface The surface provided
+  /// @param [in] navDir The navigation direction
   /// @param [in] bcheck The boundary check for this status update
-  /// @param [in] logger Logger instance to use
   /// @param [in] surfaceTolerance Surface tolerance used for intersection
+  /// @param [in] logger Logger instance to use
   Intersection3D::Status updateSurfaceStatus(
-      State& state, const Surface& surface, const BoundaryCheck& bcheck,
-      const Logger& logger = getDummyLogger(),
-      ActsScalar surfaceTolerance = s_onSurfaceTolerance) const {
+      State& state, const Surface& surface, Direction navDir,
+      const BoundaryCheck& bcheck,
+      ActsScalar surfaceTolerance = s_onSurfaceTolerance,
+      const Logger& logger = getDummyLogger()) const {
     return detail::updateSingleSurfaceStatus<AtlasStepper>(
-        *this, state, surface, bcheck, logger, surfaceTolerance);
+        *this, state, surface, navDir, bcheck, surfaceTolerance, logger);
   }
 
   /// Update step size
@@ -1113,7 +1107,7 @@ class AtlasStepper {
   Result<double> step(propagator_state_t& state,
                       const navigator_t& /*navigator*/) const {
     // we use h for keeping the nominclature with the original atlas code
-    auto h = state.stepping.stepSize.value() * state.stepping.navDir;
+    auto h = state.stepping.stepSize.value() * state.options.direction;
     bool Jac = state.stepping.useJacobian;
 
     double* R = &(state.stepping.pVector[0]);  // Coordinates
@@ -1225,7 +1219,7 @@ class AtlasStepper {
       if (std::abs(EST) > std::abs(state.options.tolerance)) {
         h = h * .5;
         // neutralize the sign of h again
-        state.stepping.stepSize.setValue(h * state.stepping.navDir);
+        state.stepping.stepSize.setValue(h * state.options.direction);
         //        dltm = 0.;
         nStepTrials++;
         continue;

Core/include/Acts/Propagator/DirectNavigator.hpp

@@ -221,8 +221,9 @@ class DirectNavigator {
     if (state.navigation.navSurfaceIter != state.navigation.navSurfaces.end()) {
       // Establish & update the surface status
       auto surfaceStatus = stepper.updateSurfaceStatus(
-          state.stepping, **state.navigation.navSurfaceIter, false, *m_logger,
-          state.options.targetTolerance);
+          state.stepping, **state.navigation.navSurfaceIter,
+          state.options.direction, false, state.options.targetTolerance,
+          *m_logger);
       if (surfaceStatus == Intersection3D::Status::unreachable) {
         ACTS_VERBOSE(
             "Surface not reachable anymore, switching to next one in "
@@ -269,8 +270,9 @@ class DirectNavigator {
     if (state.navigation.navSurfaceIter != state.navigation.navSurfaces.end()) {
       // Establish the surface status
       auto surfaceStatus = stepper.updateSurfaceStatus(
-          state.stepping, **state.navigation.navSurfaceIter, false, *m_logger,
-          state.options.targetTolerance);
+          state.stepping, **state.navigation.navSurfaceIter,
+          state.options.direction, false, state.options.targetTolerance,
+          *m_logger);
       if (surfaceStatus == Intersection3D::Status::onSurface) {
         // Set the current surface
         state.navigation.currentSurface = *state.navigation.navSurfaceIter;

Core/include/Acts/Propagator/EigenStepper.hpp

@@ -69,18 +69,15 @@ class EigenStepper {
     /// @param [in] gctx is the context object for the geometry
     /// @param [in] fieldCacheIn is the cache object for the magnetic field
     /// @param [in] par The track parameters at start
-    /// @param [in] ndir The navigation direction w.r.t momentum
     /// @param [in] ssize is the maximum step size
     ///
     /// @note the covariance matrix is copied when needed
     template <typename charge_t>
     explicit State(const GeometryContext& gctx,
                    MagneticFieldProvider::Cache fieldCacheIn,
                    const GenericBoundTrackParameters<charge_t>& par,
-                   Direction ndir = Direction::Forward,
                    double ssize = std::numeric_limits<double>::max())
         : absCharge(std::abs(par.charge())),
-          navDir(ndir),
           stepSize(ssize),
           fieldCache(std::move(fieldCacheIn)),
           geoContext(gctx) {
@@ -111,9 +108,6 @@ class EigenStepper {
     bool covTransport = false;
     Covariance cov = Covariance::Zero();
 
-    /// Navigation direction, this is needed for searching
-    Direction navDir;
-
     /// The full jacobian of the transport entire transport
     Jacobian jacobian = Jacobian::Identity();
 
@@ -168,7 +162,6 @@ class EigenStepper {
   State makeState(std::reference_wrapper<const GeometryContext> gctx,
                   std::reference_wrapper<const MagneticFieldContext> mctx,
                   const GenericBoundTrackParameters<charge_t>& par,
-                  Direction navDir = Direction::Forward,
                   double ssize = std::numeric_limits<double>::max()) const;
 
   /// @brief Resets the state
@@ -177,11 +170,10 @@ class EigenStepper {
   /// @param [in] boundParams Parameters in bound parametrisation
   /// @param [in] cov Covariance matrix
   /// @param [in] surface The reference surface of the bound parameters
-  /// @param [in] navDir Navigation direction
   /// @param [in] stepSize Step size
   void resetState(
       State& state, const BoundVector& boundParams, const BoundSymMatrix& cov,
-      const Surface& surface, const Direction navDir = Direction::Forward,
+      const Surface& surface,
       const double stepSize = std::numeric_limits<double>::max()) const;
 
   /// Get the field for the stepping, it checks first if the access is still
@@ -248,15 +240,17 @@ class EigenStepper {
   ///
   /// @param [in,out] state The stepping state (thread-local cache)
   /// @param [in] surface The surface provided
+  /// @param [in] navDir The navigation direction
   /// @param [in] bcheck The boundary check for this status update
-  /// @param [in] logger A @c Logger instance
   /// @param [in] surfaceTolerance Surface tolerance used for intersection
+  /// @param [in] logger A @c Logger instance
   Intersection3D::Status updateSurfaceStatus(
-      State& state, const Surface& surface, const BoundaryCheck& bcheck,
-      const Logger& logger = getDummyLogger(),
-      ActsScalar surfaceTolerance = s_onSurfaceTolerance) const {
+      State& state, const Surface& surface, Direction navDir,
+      const BoundaryCheck& bcheck,
+      ActsScalar surfaceTolerance = s_onSurfaceTolerance,
+      const Logger& logger = getDummyLogger()) const {
     return detail::updateSingleSurfaceStatus<EigenStepper>(
-        *this, state, surface, bcheck, logger, surfaceTolerance);
+        *this, state, surface, navDir, bcheck, surfaceTolerance, logger);
   }
 
   /// Update step size

Core/include/Acts/Propagator/EigenStepper.ipp

@@ -19,24 +19,22 @@ template <typename charge_t>
 auto Acts::EigenStepper<E, A>::makeState(
     std::reference_wrapper<const GeometryContext> gctx,
     std::reference_wrapper<const MagneticFieldContext> mctx,
-    const GenericBoundTrackParameters<charge_t>& par, Direction navDir,
-    double ssize) const -> State {
-  return State{gctx, m_bField->makeCache(mctx), par, navDir, ssize};
+    const GenericBoundTrackParameters<charge_t>& par, double ssize) const
+    -> State {
+  return State{gctx, m_bField->makeCache(mctx), par, ssize};
 }
 
 template <typename E, typename A>
 void Acts::EigenStepper<E, A>::resetState(State& state,
                                           const BoundVector& boundParams,
                                           const BoundSymMatrix& cov,
                                           const Surface& surface,
-                                          const Direction navDir,
                                           const double stepSize) const {
   // Update the stepping state
   update(state,
          detail::transformBoundToFreeParameters(surface, state.geoContext,
                                                 boundParams),
          boundParams, cov, surface);
-  state.navDir = navDir;
   state.stepSize = ConstrainedStep(stepSize);
   state.pathAccumulated = 0.;
 
@@ -147,7 +145,7 @@ Acts::Result<double> Acts::EigenStepper<E, A>::step(
     constexpr auto success = &Result<bool>::success;
     constexpr auto failure = &Result<bool>::failure;
 
-    const double h = step.value() * state.stepping.navDir;
+    const double h = step.value() * state.options.direction;
     // State the square and half of the step size
     h2 = h * h;
     half_h = h * 0.5;
@@ -230,7 +228,7 @@ Acts::Result<double> Acts::EigenStepper<E, A>::step(
   }
 
   // use the adjusted step size
-  const double h = state.stepping.stepSize.value() * state.stepping.navDir;
+  const double h = state.stepping.stepSize.value() * state.options.direction;
 
   // When doing error propagation, update the associated Jacobian matrix
   if (state.stepping.covTransport) {

Core/include/Acts/Propagator/MaterialInteractor.hpp

@@ -106,7 +106,7 @@ struct MaterialInteractor {
         stepper.transportCovarianceToCurvilinear(state.stepping);
       }
       // Change the noise updater depending on the navigation direction
-      NoiseUpdateMode mode = (state.stepping.navDir == Direction::Forward)
+      NoiseUpdateMode mode = (state.options.direction == Direction::Forward)
                                  ? addNoise
                                  : removeNoise;
       // Apply the material interactions