refactor: Remove direction from constrained step (#2073)

With this PR I want to remove the propagation direction sign inside the `ConstrainedStep` and put it into the stepper as it feels like additional burden and logic that should not be handled by the `ConstrainedStep`. This will not change anything about the overstepping mechanics.

Core/include/Acts/Navigation/NextNavigator.hpp

@@ -16,7 +16,6 @@
 #include "Acts/Geometry/GeometryIdentifier.hpp"
 #include "Acts/Geometry/Layer.hpp"
 #include "Acts/Navigation/NavigationState.hpp"
-#include "Acts/Propagator/ConstrainedStep.hpp"
 #include "Acts/Propagator/Propagator.hpp"
 #include "Acts/Surfaces/Surface.hpp"
 #include "Acts/Utilities/Logger.hpp"

Core/include/Acts/Propagator/AtlasStepper.hpp

@@ -63,7 +63,7 @@ class AtlasStepper {
           double stolerance = s_onSurfaceTolerance)
         : navDir(ndir),
           field(0., 0., 0.),
-          stepSize(ndir * std::abs(ssize)),
+          stepSize(ssize),
           tolerance(stolerance),
           fieldCache(std::move(fieldCacheIn)),
           geoContext(gctx) {
@@ -1113,7 +1113,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();
+    auto h = state.stepping.stepSize.value() * state.stepping.navDir;
     bool Jac = state.stepping.useJacobian;
 
     double* R = &(state.stepping.pVector[0]);  // Coordinates
@@ -1224,7 +1224,8 @@ class AtlasStepper {
            std::abs((C1 + C6) - (C3 + C4)));
       if (std::abs(EST) > std::abs(state.options.tolerance)) {
         h = h * .5;
-        state.stepping.stepSize.setValue(h);
+        // neutralize the sign of h again
+        state.stepping.stepSize.setValue(h * state.stepping.navDir);
         //        dltm = 0.;
         nStepTrials++;
         continue;

Core/include/Acts/Propagator/ConstrainedStep.hpp

@@ -9,7 +9,6 @@
 #pragma once
 
 #include "Acts/Definitions/Algebra.hpp"
-#include "Acts/Definitions/Direction.hpp"
 
 #include <algorithm>
 #include <array>
@@ -21,7 +20,11 @@
 
 namespace Acts {
 
-/// A constrained step class for the steppers
+/// A constrained step class for the steppers.
+///
+/// This class is symmetrical for forward and backward propagation. The sign of
+/// the propagation direction should not enter here but rather be applied the
+/// step is actually taken.
 ///
 /// As simple as this class looks it hides a few very important details:
 /// - Overstepping handling. The step size sign will flip if we happened to pass
@@ -30,9 +33,8 @@ namespace Acts {
 /// order to converge on a target.
 ///
 /// Because of the points mentioned above, the update function will always
-/// prefer step sizes that point opposite the navigation direction. A side
-/// effect of this is that we will propagate in the opposite direction if the
-/// target is "behind us".
+/// prefer negative step sizes. A side effect of this is that we will propagate
+/// in the opposite direction if the target is "behind us".
 ///
 /// The hierarchy is:
 /// - Overstepping resolution / backpropagation
@@ -56,12 +58,8 @@ class ConstrainedStep {
   constexpr ConstrainedStep() = default;
 
   /// constructor from Scalar
-  /// navigation direction is inferred by the sign of the step size
   /// @param value is the user given initial value
-  constexpr explicit ConstrainedStep(Scalar value) {
-    m_values[user] = std::abs(value);
-    m_direction = Direction::fromScalar(value);
-  }
+  constexpr explicit ConstrainedStep(Scalar value) { m_values[user] = value; }
 
   /// set accuracy by one Scalar
   ///
@@ -71,7 +69,7 @@ class ConstrainedStep {
   /// @param value is the new accuracy value
   constexpr void setValue(Scalar value) {
     /// set the accuracy value
-    m_values[accuracy] = value * m_direction;
+    m_values[accuracy] = value;
   }
 
   /// returns the min step size
@@ -80,9 +78,7 @@ class ConstrainedStep {
   /// Access a specific value
   ///
   /// @param type is the requested parameter type
-  constexpr Scalar value(Type type) const {
-    return m_values[type] * m_direction;
-  }
+  constexpr Scalar value(Type type) const { return m_values[type]; }
 
   /// Access the currently leading type
   constexpr Type currentType() const {
@@ -109,14 +105,14 @@ class ConstrainedStep {
     }
     // check the current value and set it if appropriate
     // this will also allow signed values due to overstepping
-    if (std::abs(value) <= std::abs(m_values[type])) {
-      m_values[type] = value * m_direction;
+    if (std::abs(value) < std::abs(m_values[type])) {
+      m_values[type] = value;
     }
   }
 
   constexpr void scale(Scalar factor) {
     assert(factor > 0 && "ConstrainedStep scale factor was zero or negative.");
-    m_values[accuracy] = value() * factor * m_direction;
+    m_values[accuracy] = value() * factor;
   }
 
   std::ostream& toStream(std::ostream& os) const {
@@ -154,11 +150,7 @@ class ConstrainedStep {
   inline static constexpr auto kNotSet = std::numeric_limits<Scalar>::max();
 
   /// the step size tuple
-  /// all values point in the `m_direction`
   std::array<Scalar, 4> m_values = {kNotSet, kNotSet, kNotSet, kNotSet};
-  /// the navigation direction
-  /// the direction is invariant after initialization
-  Direction m_direction = Direction::Forward;
 };
 
 inline std::ostream& operator<<(std::ostream& os, const ConstrainedStep& step) {

Core/include/Acts/Propagator/DirectNavigator.hpp

@@ -12,7 +12,6 @@
 #include "Acts/Geometry/Layer.hpp"
 #include "Acts/Geometry/TrackingGeometry.hpp"
 #include "Acts/Geometry/TrackingVolume.hpp"
-#include "Acts/Propagator/ConstrainedStep.hpp"
 #include "Acts/Propagator/Propagator.hpp"
 #include "Acts/Surfaces/Surface.hpp"
 

Core/include/Acts/Propagator/EigenStepper.hpp

@@ -16,6 +16,7 @@
 #include "Acts/Definitions/Units.hpp"
 #include "Acts/EventData/TrackParameters.hpp"
 #include "Acts/MagneticField/MagneticFieldProvider.hpp"
+#include "Acts/Propagator/ConstrainedStep.hpp"
 #include "Acts/Propagator/DefaultExtension.hpp"
 #include "Acts/Propagator/DenseEnvironmentExtension.hpp"
 #include "Acts/Propagator/EigenStepperError.hpp"
@@ -80,7 +81,7 @@ class EigenStepper {
                    double ssize = std::numeric_limits<double>::max())
         : absCharge(std::abs(par.charge())),
           navDir(ndir),
-          stepSize(ndir * std::abs(ssize)),
+          stepSize(ssize),
           fieldCache(std::move(fieldCacheIn)),
           geoContext(gctx) {
       pars.template segment<3>(eFreePos0) = par.position(gctx);

Core/include/Acts/Propagator/EigenStepper.ipp

@@ -147,7 +147,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();
+    const double h = step.value() * state.stepping.navDir;
     // State the square and half of the step size
     h2 = h * h;
     half_h = h * 0.5;
@@ -230,7 +230,7 @@ Acts::Result<double> Acts::EigenStepper<E, A>::step(
   }
 
   // use the adjusted step size
-  const double h = state.stepping.stepSize.value();
+  const double h = state.stepping.stepSize.value() * state.stepping.navDir;
 
   // When doing error propagation, update the associated Jacobian matrix
   if (state.stepping.covTransport) {

Core/include/Acts/Propagator/MultiEigenStepperLoop.hpp

@@ -18,6 +18,7 @@
 #include "Acts/EventData/TrackParameters.hpp"
 #include "Acts/EventData/detail/CorrectedTransformationFreeToBound.hpp"
 #include "Acts/MagneticField/MagneticFieldProvider.hpp"
+#include "Acts/Propagator/ConstrainedStep.hpp"
 #include "Acts/Propagator/EigenStepper.hpp"
 #include "Acts/Propagator/EigenStepperError.hpp"
 #include "Acts/Propagator/Propagator.hpp"

Core/include/Acts/Propagator/Propagator.hpp

@@ -136,7 +136,7 @@ struct PropagatorOptions : public PropagatorPlainOptions {
     eoptions.mass = mass;
     eoptions.maxSteps = maxSteps;
     eoptions.maxRungeKuttaStepTrials = maxRungeKuttaStepTrials;
-    eoptions.maxStepSize = direction * std::abs(maxStepSize);
+    eoptions.maxStepSize = maxStepSize;
     eoptions.targetTolerance = targetTolerance;
     eoptions.pathLimit = direction * std::abs(pathLimit);
     eoptions.loopProtection = loopProtection;
@@ -162,7 +162,7 @@ struct PropagatorOptions : public PropagatorPlainOptions {
     mass = pOptions.mass;
     maxSteps = pOptions.maxSteps;
     maxRungeKuttaStepTrials = pOptions.maxRungeKuttaStepTrials;
-    maxStepSize = direction * std::abs(pOptions.maxStepSize);
+    maxStepSize = pOptions.maxStepSize;
     targetTolerance = pOptions.targetTolerance;
     pathLimit = direction * std::abs(pOptions.pathLimit);
     loopProtection = pOptions.loopProtection;

Core/include/Acts/Propagator/StandardAborters.hpp

@@ -63,8 +63,8 @@ struct PathLimitReached {
       return true;
     }
     // Check if the maximum allowed step size has to be updated
-    double distance = state.stepping.navDir * std::abs(internalLimit) -
-                      state.stepping.pathAccumulated;
+    double distance =
+        std::abs(internalLimit) - std::abs(state.stepping.pathAccumulated);
     double tolerance = state.options.targetTolerance;
     stepper.setStepSize(state.stepping, distance, ConstrainedStep::aborter,
                         false);
@@ -169,8 +169,8 @@ struct SurfaceReached {
         // Update the distance to the alternative solution
         distance = sIntersection.alternative.pathLength;
       }
-      stepper.setStepSize(state.stepping, state.stepping.navDir * distance,
-                          ConstrainedStep::aborter, false);
+      stepper.setStepSize(state.stepping, distance, ConstrainedStep::aborter,
+                          false);
 
       ACTS_VERBOSE("Target: 0 | "
                    << "Target stepSize (surface) updated to "

Core/include/Acts/Propagator/StraightLineStepper.hpp

@@ -79,7 +79,7 @@ class StraightLineStepper {
                    double stolerance = s_onSurfaceTolerance)
         : absCharge(std::abs(par.charge())),
           navDir(ndir),
-          stepSize(ndir * std::abs(ssize)),
+          stepSize(ssize),
           tolerance(stolerance),
           geoContext(gctx) {
       (void)mctx;
@@ -395,7 +395,7 @@ class StraightLineStepper {
   Result<double> step(propagator_state_t& state,
                       const navigator_t& /*navigator*/) const {
     // use the adjusted step size
-    const auto h = state.stepping.stepSize.value();
+    const auto h = state.stepping.stepSize.value() * state.stepping.navDir;
     const double p = absoluteMomentum(state.stepping);
     // time propagates along distance as 1/b = sqrt(1 + m²/p²)
     const auto dtds = std::hypot(1., state.options.mass / p);

Core/include/Acts/Propagator/detail/SteppingHelper.hpp

@@ -57,17 +57,15 @@ Acts::Intersection3D::Status updateSingleSurfaceStatus(
     if (detail::checkIntersection(sIntersection.intersection, pLimit, oLimit,
                                   surfaceTolerance, logger)) {
       ACTS_VERBOSE("Surface is reachable");
-      stepper.setStepSize(state,
-                          state.navDir * sIntersection.intersection.pathLength);
+      stepper.setStepSize(state, sIntersection.intersection.pathLength);
       return Intersection3D::Status::reachable;
     }
 
     if (sIntersection.alternative and
         detail::checkIntersection(sIntersection.alternative, pLimit, oLimit,
                                   surfaceTolerance, logger)) {
       ACTS_VERBOSE("Surface is reachable");
-      stepper.setStepSize(state,
-                          state.navDir * sIntersection.alternative.pathLength);
+      stepper.setStepSize(state, sIntersection.alternative.pathLength);
       return Intersection3D::Status::reachable;
     }
   }

Core/include/Acts/Propagator/detail/SteppingLogger.hpp

@@ -10,6 +10,7 @@
 
 #include "Acts/Definitions/Algebra.hpp"
 #include "Acts/Definitions/Common.hpp"
+#include "Acts/Definitions/Direction.hpp"
 #include "Acts/Geometry/GeometryIdentifier.hpp"
 #include "Acts/Propagator/ConstrainedStep.hpp"
 #include "Acts/Utilities/Logger.hpp"
@@ -28,6 +29,7 @@ namespace detail {
 /// @brief the step information for
 struct Step {
   ConstrainedStep stepSize;
+  Direction navDir;
   Vector3 position = Vector3(0., 0., 0.);
   Vector3 momentum = Vector3(0., 0., 0.);
   std::shared_ptr<const Surface> surface = nullptr;
@@ -70,6 +72,7 @@ struct SteppingLogger {
     // record the propagation state
     Step step;
     step.stepSize = state.stepping.stepSize;
+    step.navDir = state.stepping.navDir;
     step.position = stepper.position(state.stepping);
     step.momentum = stepper.momentum(state.stepping);
 

Core/include/Acts/TrackFinding/CombinatorialKalmanFilter.hpp

@@ -1226,8 +1226,7 @@ class CombinatorialKalmanFilter {
       state.stepping.jacTransport = FreeMatrix::Identity();
       state.stepping.derivative = FreeVector::Zero();
       // Reset the step size
-      state.stepping.stepSize = ConstrainedStep(
-          state.stepping.navDir * std::abs(state.options.maxStepSize));
+      state.stepping.stepSize = ConstrainedStep(state.options.maxStepSize);
       // Set accumulatd path to zero before targeting surface
       state.stepping.pathAccumulated = 0.;
 

Core/include/Acts/TrackFitting/KalmanFitter.hpp

@@ -515,9 +515,7 @@ class KalmanFitter {
       state.stepping.navDir = state.stepping.navDir.invert();
 
       // Reset propagator options
-      state.options.maxStepSize =
-          state.stepping.navDir * std::abs(state.options.maxStepSize);
-      // Not sure if reset of pathLimit during propagation makes any sense
+      // TODO Not sure if reset of pathLimit during propagation makes any sense
       state.options.pathLimit =
           state.stepping.navDir * std::abs(state.options.pathLimit);
 
@@ -978,8 +976,7 @@ class KalmanFitter {
         state.stepping.navDir = state.stepping.navDir.invert();
       }
       // Reset the step size
-      state.stepping.stepSize = ConstrainedStep(
-          state.stepping.navDir * std::abs(state.options.maxStepSize));
+      state.stepping.stepSize = ConstrainedStep(state.options.maxStepSize);
       // Set accumulatd path to zero before targeting surface
       state.stepping.pathAccumulated = 0.;
 

Core/src/Geometry/Layer.cpp

@@ -289,15 +289,12 @@ Acts::SurfaceIntersection Acts::Layer::surfaceOnApproach(
 
     if (detail::checkIntersection(isection.intersection, pLimit, oLimit,
                                   s_onSurfaceTolerance)) {
-      isection.intersection.pathLength *= options.navDir;
       return isection;
     }
 
     if (isection.alternative and
         detail::checkIntersection(isection.alternative, pLimit, oLimit,
                                   s_onSurfaceTolerance)) {
-      // Set the right sign for the path length
-      isection.alternative.pathLength *= options.navDir;
       return SurfaceIntersection(isection.alternative, isection.object);
     }