diff --git a/Core/include/Acts/Propagator/DefaultExtension.hpp b/Core/include/Acts/Propagator/DefaultExtension.hpp
index c111185379253ed246851c52e9cea699dc90faf9..c1cddc578fa76f157ef9218c24fb353b33a8a660 100644
--- a/Core/include/Acts/Propagator/DefaultExtension.hpp
+++ b/Core/include/Acts/Propagator/DefaultExtension.hpp
@@ -38,19 +38,22 @@ struct DefaultExtension {
/// @param [in] stepper Stepper of the propagation
/// @param [out] knew Next k_i that is evaluated
/// @param [in] bField B-Field at the evaluation position
+ /// @param [out] kQoP k_i elements of the momenta
/// @param [in] i Index of the k_i, i = [0, 3]
/// @param [in] h Step size (= 0. ^ 0.5 * StepSize ^ StepSize)
/// @param [in] kprev Evaluated k_{i - 1}
/// @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,
- const double h = 0., const Vector3D& kprev = Vector3D()) {
+ 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 3492fa815deae7ba76b3bb3fa8b5d30d9df5ece3..9608a7dcd7da622fb9d6cb7a0e391cc31cde5036 100644
--- a/Core/include/Acts/Propagator/DenseEnvironmentExtension.hpp
+++ b/Core/include/Acts/Propagator/DenseEnvironmentExtension.hpp
@@ -79,6 +79,7 @@ struct DenseEnvironmentExtension {
/// @tparam stepper_t Type of the stepper
/// @param [in] state State of the propagator
/// @param [out] knew Next k_i that is evaluated
+ /// @param [out] kQoP k_i elements of the momenta
/// @param [in] bField B-Field at the evaluation position
/// @param [in] i Index of the k_i, i = [0, 3]
/// @param [in] h Step size (= 0. ^ 0.5 * StepSize ^ StepSize)
@@ -86,8 +87,9 @@ 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,
- const double h = 0., const Vector3D& kprev = Vector3D()) {
+ 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) {
// Set up container for energy loss
@@ -106,6 +108,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 +125,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;
}
@@ -473,4 +477,4 @@ struct DenseStepperPropagatorOptions
}
};
-} // namespace Acts
+} // namespace Acts
\ No newline at end of file
diff --git a/Core/include/Acts/Propagator/EigenStepper.hpp b/Core/include/Acts/Propagator/EigenStepper.hpp
index d3f493ed1bf17eb78b187cdfe873157c8421f6be..16483628967ede95fa285b67cd30795437099364 100644
--- a/Core/include/Acts/Propagator/EigenStepper.hpp
+++ b/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;
};
@@ -366,4 +368,4 @@ class EigenStepper {
};
} // namespace Acts
-#include "Acts/Propagator/EigenStepper.ipp"
+#include "Acts/Propagator/EigenStepper.ipp"
\ No newline at end of file
diff --git a/Core/include/Acts/Propagator/EigenStepper.ipp b/Core/include/Acts/Propagator/EigenStepper.ipp
index a07d0cc9bb391d3de1024248303249592a9ce6f2..97329488c49c30612e67af39cf67e63e27de2e96 100644
--- a/Core/include/Acts/Propagator/EigenStepper.ipp
+++ b/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,14 +233,14 @@ 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,
- sd.k1)) {
+ 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,
- sd.k2)) {
+ if (!state.stepping.extension.k3(state, *this, sd.k3, sd.B_middle, sd.kQoP,
+ half_h, sd.k2)) {
return false;
}
@@ -248,14 +248,16 @@ 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);
};
@@ -321,4 +323,4 @@ Acts::Result<double> Acts::EigenStepper<B, E, A>::step(
}
state.stepping.pathAccumulated += h;
return h;
-}
+}
\ No newline at end of file
diff --git a/Core/include/Acts/Propagator/StepperExtensionList.hpp b/Core/include/Acts/Propagator/StepperExtensionList.hpp
index d6de656d48c9831e83e54bb705e70f60ff5845e5..cbbc133cb52b0f38974a9c1d7a5459de8397e8af 100644
--- a/Core/include/Acts/Propagator/StepperExtensionList.hpp
+++ b/Core/include/Acts/Propagator/StepperExtensionList.hpp
@@ -79,8 +79,9 @@ 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,10 +89,10 @@ 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,
- const Vector3D& kprev) {
- return impl::k(tuple(), state, stepper, knew, bField, validExtensions, 1, h,
- kprev);
+ Vector3D& knew, const Vector3D& bField, std::array<double, 4>& kQoP,
+ const double h, const Vector3D& kprev) {
+ return impl::k(tuple(), state, stepper, knew, bField, kQoP, validExtensions,
+ 1, h, kprev);
}
/// @brief This functions broadcasts the call for evaluating k3. It collects
@@ -99,10 +100,10 @@ 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,
- const Vector3D& kprev) {
- return impl::k(tuple(), state, stepper, knew, bField, validExtensions, 2, h,
- kprev);
+ Vector3D& knew, const Vector3D& bField, std::array<double, 4>& kQoP,
+ const double h, const Vector3D& kprev) {
+ return impl::k(tuple(), state, stepper, knew, bField, kQoP, validExtensions,
+ 2, h, kprev);
}
/// @brief This functions broadcasts the call for evaluating k4. It collects
@@ -110,10 +111,10 @@ 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,
- const Vector3D& kprev) {
- return impl::k(tuple(), state, stepper, knew, bField, validExtensions, 3, h,
- kprev);
+ Vector3D& knew, const Vector3D& bField, std::array<double, 4>& kQoP,
+ const double h, const Vector3D& kprev) {
+ return impl::k(tuple(), state, stepper, knew, bField, kQoP, validExtensions,
+ 3, h, kprev);
}
/// @brief This functions broadcasts the call of the method finalize(). It
@@ -135,4 +136,4 @@ struct StepperExtensionList : private detail::Extendable<extensions...> {
}
};
-} // namespace Acts
+} // namespace Acts
\ No newline at end of file
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 0345eaab468d15cca5fa9e03bc2f4d36a89ab2b7..ea4adcc36aac8465a8833a2e8b967e6c8d42a491 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*/) {
diff --git a/Tests/UnitTests/Core/Propagator/StepperTests.cpp b/Tests/UnitTests/Core/Propagator/StepperTests.cpp
index 2836e6178e0b217ff249269215e5d712e3853702..a464f136367b45582b8ab99a6549823594cfe454 100644
--- a/Tests/UnitTests/Core/Propagator/StepperTests.cpp
+++ b/Tests/UnitTests/Core/Propagator/StepperTests.cpp
@@ -576,7 +576,7 @@ BOOST_AUTO_TEST_CASE(step_extension_vacmatvac_test) {
propOptsDense(tgContext, mfContext);
abortList.get<EndOfWorld>().maxX = 2_m;
propOptsDense.abortList = abortList;
- propOptsDense.maxSteps = 100;
+ propOptsDense.maxSteps = 1000;
propOptsDense.maxStepSize = 1.5_m;
propOptsDense.tolerance = 1e-8;