From 1ea1bb4c71563f668bc55a12b18738c0bbb43a8e Mon Sep 17 00:00:00 2001
From: Fabian Klimpel <fklimpel@cern.ch>
Date: Wed, 18 Mar 2020 13:07:24 +0100
Subject: [PATCH] Adapted extension interfaces

---
 .../include/Acts/Propagator/DefaultExtension.hpp |  3 ++-
 .../Propagator/DenseEnvironmentExtension.hpp     |  4 +++-
 Core/include/Acts/Propagator/EigenStepper.ipp    |  2 +-
 .../Acts/Propagator/StepperExtensionList.hpp     | 16 ++++++++--------
 .../stepper_extension_list_implementation.hpp    | 10 +++++-----
 5 files changed, 19 insertions(+), 16 deletions(-)

diff --git a/Core/include/Acts/Propagator/DefaultExtension.hpp b/Core/include/Acts/Propagator/DefaultExtension.hpp
index c11118537..9c8fda4e5 100644
--- a/Core/include/Acts/Propagator/DefaultExtension.hpp
+++ b/Core/include/Acts/Propagator/DefaultExtension.hpp
@@ -44,13 +44,14 @@ struct DefaultExtension {
   /// @return Boolean flag if the calculation is valid
   template <typename propagator_state_t, typename stepper_t>
   bool k(const propagator_state_t& state, const stepper_t& stepper,
-         Vector3D& knew, const Vector3D& bField, const int i = 0,
+         Vector3D& knew, const Vector3D& bField, std::array<double, 4>& kQoP, const int i = 0,
          const double h = 0., const Vector3D& kprev = Vector3D()) {
     auto qop =
         stepper.charge(state.stepping) / stepper.momentum(state.stepping);
     // First step does not rely on previous data
     if (i == 0) {
       knew = qop * stepper.direction(state.stepping).cross(bField);
+      kQoP = {0., 0., 0., 0.};
     } else {
       knew =
           qop * (stepper.direction(state.stepping) + h * kprev).cross(bField);
diff --git a/Core/include/Acts/Propagator/DenseEnvironmentExtension.hpp b/Core/include/Acts/Propagator/DenseEnvironmentExtension.hpp
index 3492fa815..350a8f4df 100644
--- a/Core/include/Acts/Propagator/DenseEnvironmentExtension.hpp
+++ b/Core/include/Acts/Propagator/DenseEnvironmentExtension.hpp
@@ -86,7 +86,7 @@ struct DenseEnvironmentExtension {
   /// @return Boolean flag if the calculation is valid
   template <typename propagator_state_t, typename stepper_t>
   bool k(const propagator_state_t& state, const stepper_t& stepper,
-         Vector3D& knew, const Vector3D& bField, const int i = 0,
+         Vector3D& knew, const Vector3D& bField, std::array<double, 4>& kQoP, const int i = 0,
          const double h = 0., const Vector3D& kprev = Vector3D()) {
     // i = 0 is used for setup and evaluation of k
     if (i == 0) {
@@ -106,6 +106,7 @@ struct DenseEnvironmentExtension {
           (stepper.charge(state.stepping) * stepper.charge(state.stepping));
       //~ tKi[0] = std::hypot(1, state.options.mass / initialMomentum);
       tKi[0] = std::hypot(1, state.options.mass * qop[0]);
+      kQoP[0] = Lambdappi[0];
     } else {
       // Update parameters and check for momentum condition
       updateEnergyLoss(state.options.mass, h, state.stepping, stepper, i);
@@ -122,6 +123,7 @@ struct DenseEnvironmentExtension {
           (stepper.charge(state.stepping) * stepper.charge(state.stepping) *
            UnitConstants::C * UnitConstants::C);
       tKi[i] = std::hypot(1, state.options.mass * qopNew);
+      kQoP[i] = Lambdappi[i];
     }
     return true;
   }
diff --git a/Core/include/Acts/Propagator/EigenStepper.ipp b/Core/include/Acts/Propagator/EigenStepper.ipp
index 6abf968a0..9ef674b6e 100644
--- a/Core/include/Acts/Propagator/EigenStepper.ipp
+++ b/Core/include/Acts/Propagator/EigenStepper.ipp
@@ -255,7 +255,7 @@ Acts::Result<double> Acts::EigenStepper<B, E, A>::step(
 
     // Compute and check the local integration error estimate
     error_estimate = std::max(
-        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);
+        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);
   };
 
diff --git a/Core/include/Acts/Propagator/StepperExtensionList.hpp b/Core/include/Acts/Propagator/StepperExtensionList.hpp
index d6de656d4..8ba023a02 100644
--- a/Core/include/Acts/Propagator/StepperExtensionList.hpp
+++ b/Core/include/Acts/Propagator/StepperExtensionList.hpp
@@ -79,8 +79,8 @@ struct StepperExtensionList : private detail::Extendable<extensions...> {
   /// the evaluation is valid.
   template <typename propagator_state_t, typename stepper_t>
   bool k1(const propagator_state_t& state, const stepper_t& stepper,
-          Vector3D& knew, const Vector3D& bField) {
-    return impl::k(tuple(), state, stepper, knew, bField, validExtensions);
+          Vector3D& knew, const Vector3D& bField, std::array<double, 4>& kQoP) {
+    return impl::k(tuple(), state, stepper, knew, bField, kQoP, validExtensions);
   }
 
   /// @brief This functions broadcasts the call for evaluating k2. It collects
@@ -88,9 +88,9 @@ struct StepperExtensionList : private detail::Extendable<extensions...> {
   /// returns a boolean as indicator if the evaluation is valid.
   template <typename propagator_state_t, typename stepper_t>
   bool k2(const propagator_state_t& state, const stepper_t& stepper,
-          Vector3D& knew, const Vector3D& bField, const double h,
+          Vector3D& knew, const Vector3D& bField, std::array<double, 4>& kQoP, const double h,
           const Vector3D& kprev) {
-    return impl::k(tuple(), state, stepper, knew, bField, validExtensions, 1, h,
+    return impl::k(tuple(), state, stepper, knew, bField, kQoP, validExtensions, 1, h,
                    kprev);
   }
 
@@ -99,9 +99,9 @@ struct StepperExtensionList : private detail::Extendable<extensions...> {
   /// returns a boolean as indicator if the evaluation is valid.
   template <typename propagator_state_t, typename stepper_t>
   bool k3(const propagator_state_t& state, const stepper_t& stepper,
-          Vector3D& knew, const Vector3D& bField, const double h,
+          Vector3D& knew, const Vector3D& bField, std::array<double, 4>& kQoP, const double h,
           const Vector3D& kprev) {
-    return impl::k(tuple(), state, stepper, knew, bField, validExtensions, 2, h,
+    return impl::k(tuple(), state, stepper, knew, bField, kQoP, validExtensions, 2, h,
                    kprev);
   }
 
@@ -110,9 +110,9 @@ struct StepperExtensionList : private detail::Extendable<extensions...> {
   /// returns a boolean as indicator if the evaluation is valid.
   template <typename propagator_state_t, typename stepper_t>
   bool k4(const propagator_state_t& state, const stepper_t& stepper,
-          Vector3D& knew, const Vector3D& bField, const double h,
+          Vector3D& knew, const Vector3D& bField, std::array<double, 4>& kQoP, const double h,
           const Vector3D& kprev) {
-    return impl::k(tuple(), state, stepper, knew, bField, validExtensions, 3, h,
+    return impl::k(tuple(), state, stepper, knew, bField, kQoP, validExtensions, 3, h,
                    kprev);
   }
 
diff --git a/Core/include/Acts/Propagator/detail/stepper_extension_list_implementation.hpp b/Core/include/Acts/Propagator/detail/stepper_extension_list_implementation.hpp
index 0345eaab4..ea4adcc36 100644
--- a/Core/include/Acts/Propagator/detail/stepper_extension_list_implementation.hpp
+++ b/Core/include/Acts/Propagator/detail/stepper_extension_list_implementation.hpp
@@ -55,22 +55,22 @@ struct stepper_extension_list_impl {
   template <typename propagator_state_t, typename stepper_t, typename... T>
   static bool k(std::tuple<T...>& obs_tuple, const propagator_state_t& state,
                 const stepper_t& stepper, Vector3D& knew,
-                const Vector3D& bField,
+                const Vector3D& bField, std::array<double, 4>& kQoP,
                 const std::array<bool, sizeof...(T)>& validExtensions,
                 const int i = 0, const double h = 0,
                 const Vector3D& kprev = Vector3D()) {
     // If element is invalid: continue
     if (!std::get<N - 1>(validExtensions)) {
       return stepper_extension_list_impl<N - 1>::k(
-          obs_tuple, state, stepper, knew, bField, validExtensions, i, h,
+          obs_tuple, state, stepper, knew, bField, kQoP, validExtensions, i, h,
           kprev);
     }
 
     // Continue as long as evaluations are 'true'
-    if (std::get<N - 1>(obs_tuple).k(state, stepper, knew, bField, i, h,
+    if (std::get<N - 1>(obs_tuple).k(state, stepper, knew, bField, kQoP, i, h,
                                      kprev)) {
       return stepper_extension_list_impl<N - 1>::k(
-          obs_tuple, state, stepper, knew, bField, validExtensions, i, h,
+          obs_tuple, state, stepper, knew, bField, kQoP, validExtensions, i, h,
           kprev);
     } else {
       // Break at false
@@ -142,7 +142,7 @@ struct stepper_extension_list_impl<0u> {
   static bool k(std::tuple<T...>& /*unused*/,
                 const propagator_state_t& /*unused*/,
                 const stepper_t& /*unused*/, Vector3D& /*unused*/,
-                const Vector3D& /*unused*/,
+                const Vector3D& /*unused*/, std::array<double, 4>& /*unused*/,
                 const std::array<bool, sizeof...(T)>& /*unused*/,
                 const int /*unused*/, const double /*unused*/,
                 const Vector3D& /*unused*/) {
-- 
GitLab