Set EigenStepper components

Core/include/Acts/Propagator/EigenStepper.hpp

@@ -164,6 +164,8 @@ class EigenStepper {
       Vector3D B_first, B_middle, B_last;
       /// k_i of the RKN4 algorithm
       Vector3D k1, k2, k3, k4;
+      /// k_i elements of the momenta
+       std::array<double, 4> kQoP;
     } stepData;
   };
 

Core/include/Acts/Propagator/EigenStepper.ipp

@@ -216,7 +216,7 @@ Acts::Result<double> Acts::EigenStepper<B, E, A>::step(
   // First Runge-Kutta point (at current position)
   sd.B_first = getField(state.stepping, state.stepping.pos);
   if (!state.stepping.extension.validExtensionForStep(state, *this) ||
-      !state.stepping.extension.k1(state, *this, sd.k1, sd.B_first)) {
+      !state.stepping.extension.k1(state, *this, sd.k1, sd.B_first, sd.kQoP)) {
     return 0.;
   }
 
@@ -233,13 +233,13 @@ Acts::Result<double> Acts::EigenStepper<B, E, A>::step(
     const Vector3D pos1 =
         state.stepping.pos + half_h * state.stepping.dir + h2 * 0.125 * sd.k1;
     sd.B_middle = getField(state.stepping, pos1);
-    if (!state.stepping.extension.k2(state, *this, sd.k2, sd.B_middle, half_h,
+    if (!state.stepping.extension.k2(state, *this, sd.k2, sd.B_middle, sd.kQoP, half_h,
                                      sd.k1)) {
       return false;
     }
 
     // Third Runge-Kutta point
-    if (!state.stepping.extension.k3(state, *this, sd.k3, sd.B_middle, half_h,
+    if (!state.stepping.extension.k3(state, *this, sd.k3, sd.B_middle, sd.kQoP, half_h,
                                      sd.k2)) {
       return false;
     }
@@ -248,14 +248,14 @@ Acts::Result<double> Acts::EigenStepper<B, E, A>::step(
     const Vector3D pos2 =
         state.stepping.pos + h * state.stepping.dir + h2 * 0.5 * sd.k3;
     sd.B_last = getField(state.stepping, pos2);
-    if (!state.stepping.extension.k4(state, *this, sd.k4, sd.B_last, h,
+    if (!state.stepping.extension.k4(state, *this, sd.k4, sd.B_last, sd.kQoP, h,
                                      sd.k3)) {
       return false;
     }
 
     // Compute and check the local integration error estimate
     error_estimate = std::max(
-        h2 * (sd.k1 - sd.k2 - sd.k3 + sd.k4).template lpNorm<1>(), 1e-20);
+        h2 * ((sd.k1 - sd.k2 - sd.k3 + sd.k4).template lpNorm<1>() + std::abs(sd.kQoP[0] + sd.kQoP[1] + sd.kQoP[2] + sd.kQoP[3]), 1e-20);
     return (error_estimate <= state.options.tolerance);
   };