fix: use correct path length derivatives when computing the curviline…

…ar covariance for a zero step propagation (acts-project#2910)

The propagate can use curvilinear parametrisation for the returned parameters and covariance, However if the propagation step size is below the limit, the propagator does not compute the path length derivatives which leads to an incorrect covariance matrix when enabling curvilinear parameterization.

Core/include/Acts/Propagator/AtlasStepper.hpp

@@ -510,6 +510,21 @@ class AtlasStepper {
                       state.pathAccumulated);
   }
 
+  /// @brief If necessary fill additional members needed for curvilinearState
+  ///
+  /// Compute path length derivatives in case they have not been computed
+  /// yet, which is the case if no step has been executed yet.
+  ///
+  /// @param [in, out] prop_state State that will be presented as @c BoundState
+  /// @param [in] navigator the navigator of the propagation
+  /// @return true if nothing is missing after this call, false otherwise.
+  template <typename propagator_state_t, typename navigator_t>
+  bool prepareCurvilinearState(
+      [[maybe_unused]] propagator_state_t& prop_state,
+      [[maybe_unused]] const navigator_t& navigator) const {
+    return true;
+  }
+
   /// Create and return a curvilinear state at the current position
   ///
   ///

Core/include/Acts/Propagator/EigenStepper.hpp

@@ -339,6 +339,18 @@ class EigenStepper {
       const FreeToBoundCorrection& freeToBoundCorrection =
           FreeToBoundCorrection(false)) const;
 
+  /// @brief If necessary fill additional members needed for curvilinearState
+  ///
+  /// Compute path length derivatives in case they have not been computed
+  /// yet, which is the case if no step has been executed yet.
+  ///
+  /// @param [in, out] prop_state State that will be presented as @c BoundState
+  /// @param [in] navigator the navigator of the propagation
+  /// @return true if nothing is missing after this call, false otherwise.
+  template <typename propagator_state_t, typename navigator_t>
+  bool prepareCurvilinearState(propagator_state_t& prop_state,
+                               const navigator_t& navigator) const;
+
   /// Create and return a curvilinear state at the current position
   ///
   /// @brief This transports (if necessary) the covariance

Core/include/Acts/Propagator/EigenStepper.ipp

@@ -57,6 +57,43 @@ auto Acts::EigenStepper<E, A>::boundState(
       freeToBoundCorrection);
 }
 
+template <typename E, typename A>
+template <typename propagator_state_t, typename navigator_t>
+bool Acts::EigenStepper<E, A>::prepareCurvilinearState(
+    propagator_state_t& prop_state, const navigator_t& navigator) const {
+  // test whether the accumulated path has still its initial value.
+  if (prop_state.stepping.pathAccumulated == 0.) {
+    // if no step was executed the path length derivates have not been
+    // computed but are needed to compute the curvilinear covariance. The
+    // derivates are given by k1 for a zero step width.
+    if (prop_state.stepping.extension.validExtensionForStep(prop_state, *this,
+                                                            navigator)) {
+      // First Runge-Kutta point (at current position)
+      auto& sd = prop_state.stepping.stepData;
+      auto pos = position(prop_state.stepping);
+      auto fieldRes = getField(prop_state.stepping, pos);
+      if (fieldRes.ok()) {
+        sd.B_first = *fieldRes;
+        if (prop_state.stepping.extension.k1(prop_state, *this, navigator,
+                                             sd.k1, sd.B_first, sd.kQoP)) {
+          // dr/ds :
+          prop_state.stepping.derivative.template head<3>() =
+              prop_state.stepping.pars.template segment<3>(eFreeDir0);
+          // d (dr/ds) / ds :
+          prop_state.stepping.derivative.template segment<3>(4) = sd.k1;
+          // to set dt/ds :
+          prop_state.stepping.extension.finalize(
+              prop_state, *this, navigator,
+              prop_state.stepping.pathAccumulated);
+          return true;
+        }
+      }
+    }
+    return false;
+  }
+  return true;
+}
+
 template <typename E, typename A>
 auto Acts::EigenStepper<E, A>::curvilinearState(State& state,
                                                 bool transportCov) const

Core/include/Acts/Propagator/MultiEigenStepperLoop.hpp

@@ -759,6 +759,21 @@ class MultiEigenStepperLoop
       const FreeToBoundCorrection& freeToBoundCorrection =
           FreeToBoundCorrection(false)) const;
 
+  /// @brief If necessary fill additional members needed for curvilinearState
+  ///
+  /// Compute path length derivatives in case they have not been computed
+  /// yet, which is the case if no step has been executed yet.
+  ///
+  /// @param [in, out] prop_state State that will be presented as @c BoundState
+  /// @param [in] navigator the navigator of the propagation
+  /// @return true if nothing is missing after this call, false otherwise.
+  template <typename propagator_state_t, typename navigator_t>
+  bool prepareCurvilinearState(
+      [[maybe_unused]] propagator_state_t& prop_state,
+      [[maybe_unused]] const navigator_t& navigator) const {
+    return true;
+  }
+
   /// Create and return a curvilinear state at the current position
   ///
   /// @brief This transports (if necessary) the covariance

Core/include/Acts/Propagator/Propagator.ipp

@@ -258,6 +258,10 @@ auto Acts::Propagator<S, N>::makeResult(propagator_state_t state,
   moveStateToResult(state, result);
 
   if (makeCurvilinear) {
+    if (!m_stepper.prepareCurvilinearState(state, m_navigator)) {
+      // information to compute curvilinearState is incomplete.
+      return propagationResult.error();
+    }
     /// Convert into return type and fill the result object
     auto curvState = m_stepper.curvilinearState(state.stepping);
     // Fill the end parameters

Core/include/Acts/Propagator/StraightLineStepper.hpp

@@ -319,6 +319,36 @@ class StraightLineStepper {
       const FreeToBoundCorrection& freeToBoundCorrection =
           FreeToBoundCorrection(false)) const;
 
+  /// @brief If necessary fill additional members needed for curvilinearState
+  ///
+  /// Compute path length derivatives in case they have not been computed
+  /// yet, which is the case if no step has been executed yet.
+  ///
+  /// @param [in, out] prop_state State that will be presented as @c BoundState
+  /// @param [in] navigator the navigator of the propagation
+  /// @return true if nothing is missing after this call, false otherwise.
+  template <typename propagator_state_t, typename navigator_t>
+  bool prepareCurvilinearState(
+      [[maybe_unused]] propagator_state_t& prop_state,
+      [[maybe_unused]] const navigator_t& navigator) const {
+    // test whether the accumulated path has still its initial value.
+    if (prop_state.stepping.pathAccumulated == 0.) {
+      // dr/ds :
+      prop_state.stepping.derivative.template head<3>() =
+          direction(prop_state.stepping);
+      // dt / ds
+      prop_state.stepping.derivative(eFreeTime) =
+          std::hypot(1., prop_state.stepping.particleHypothesis.mass() /
+                             absoluteMomentum(prop_state.stepping));
+      // d (dr/ds) / ds : == 0
+      prop_state.stepping.derivative.template segment<3>(4) =
+          Acts::Vector3::Zero().transpose();
+      // d qop / ds  == 0
+      prop_state.stepping.derivative(eFreeQOverP) = 0.;
+    }
+    return true;
+  }
+
   /// Create and return a curvilinear state at the current position
   ///
   /// @brief This creates a curvilinear state.