Skip to content
Snippets Groups Projects
Commit 1ea1bb4c authored by Fabian Klimpel's avatar Fabian Klimpel Committed by Fabian Klimpel
Browse files

Adapted extension interfaces

parent b2f73901
No related branches found
No related tags found
1 merge request!791Include momentum in error estimation of the EigenStepper
...@@ -44,13 +44,14 @@ struct DefaultExtension { ...@@ -44,13 +44,14 @@ struct DefaultExtension {
/// @return Boolean flag if the calculation is valid /// @return Boolean flag if the calculation is valid
template <typename propagator_state_t, typename stepper_t> template <typename propagator_state_t, typename stepper_t>
bool k(const propagator_state_t& state, const stepper_t& stepper, 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()) { const double h = 0., const Vector3D& kprev = Vector3D()) {
auto qop = auto qop =
stepper.charge(state.stepping) / stepper.momentum(state.stepping); stepper.charge(state.stepping) / stepper.momentum(state.stepping);
// First step does not rely on previous data // First step does not rely on previous data
if (i == 0) { if (i == 0) {
knew = qop * stepper.direction(state.stepping).cross(bField); knew = qop * stepper.direction(state.stepping).cross(bField);
kQoP = {0., 0., 0., 0.};
} else { } else {
knew = knew =
qop * (stepper.direction(state.stepping) + h * kprev).cross(bField); qop * (stepper.direction(state.stepping) + h * kprev).cross(bField);
......
...@@ -86,7 +86,7 @@ struct DenseEnvironmentExtension { ...@@ -86,7 +86,7 @@ struct DenseEnvironmentExtension {
/// @return Boolean flag if the calculation is valid /// @return Boolean flag if the calculation is valid
template <typename propagator_state_t, typename stepper_t> template <typename propagator_state_t, typename stepper_t>
bool k(const propagator_state_t& state, const stepper_t& stepper, 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()) { const double h = 0., const Vector3D& kprev = Vector3D()) {
// i = 0 is used for setup and evaluation of k // i = 0 is used for setup and evaluation of k
if (i == 0) { if (i == 0) {
...@@ -106,6 +106,7 @@ struct DenseEnvironmentExtension { ...@@ -106,6 +106,7 @@ struct DenseEnvironmentExtension {
(stepper.charge(state.stepping) * stepper.charge(state.stepping)); (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 / initialMomentum);
tKi[0] = std::hypot(1, state.options.mass * qop[0]); tKi[0] = std::hypot(1, state.options.mass * qop[0]);
kQoP[0] = Lambdappi[0];
} else { } else {
// Update parameters and check for momentum condition // Update parameters and check for momentum condition
updateEnergyLoss(state.options.mass, h, state.stepping, stepper, i); updateEnergyLoss(state.options.mass, h, state.stepping, stepper, i);
...@@ -122,6 +123,7 @@ struct DenseEnvironmentExtension { ...@@ -122,6 +123,7 @@ struct DenseEnvironmentExtension {
(stepper.charge(state.stepping) * stepper.charge(state.stepping) * (stepper.charge(state.stepping) * stepper.charge(state.stepping) *
UnitConstants::C * UnitConstants::C); UnitConstants::C * UnitConstants::C);
tKi[i] = std::hypot(1, state.options.mass * qopNew); tKi[i] = std::hypot(1, state.options.mass * qopNew);
kQoP[i] = Lambdappi[i];
} }
return true; return true;
} }
......
...@@ -255,7 +255,7 @@ Acts::Result<double> Acts::EigenStepper<B, E, A>::step( ...@@ -255,7 +255,7 @@ Acts::Result<double> Acts::EigenStepper<B, E, A>::step(
// Compute and check the local integration error estimate // Compute and check the local integration error estimate
error_estimate = std::max( 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); return (error_estimate <= state.options.tolerance);
}; };
......
...@@ -79,8 +79,8 @@ struct StepperExtensionList : private detail::Extendable<extensions...> { ...@@ -79,8 +79,8 @@ struct StepperExtensionList : private detail::Extendable<extensions...> {
/// the evaluation is valid. /// the evaluation is valid.
template <typename propagator_state_t, typename stepper_t> template <typename propagator_state_t, typename stepper_t>
bool k1(const propagator_state_t& state, const stepper_t& stepper, bool k1(const propagator_state_t& state, const stepper_t& stepper,
Vector3D& knew, const Vector3D& bField) { Vector3D& knew, const Vector3D& bField, std::array<double, 4>& kQoP) {
return impl::k(tuple(), state, stepper, knew, bField, validExtensions); return impl::k(tuple(), state, stepper, knew, bField, kQoP, validExtensions);
} }
/// @brief This functions broadcasts the call for evaluating k2. It collects /// @brief This functions broadcasts the call for evaluating k2. It collects
...@@ -88,9 +88,9 @@ struct StepperExtensionList : private detail::Extendable<extensions...> { ...@@ -88,9 +88,9 @@ struct StepperExtensionList : private detail::Extendable<extensions...> {
/// returns a boolean as indicator if the evaluation is valid. /// returns a boolean as indicator if the evaluation is valid.
template <typename propagator_state_t, typename stepper_t> template <typename propagator_state_t, typename stepper_t>
bool k2(const propagator_state_t& state, const stepper_t& stepper, 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) { 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); kprev);
} }
...@@ -99,9 +99,9 @@ struct StepperExtensionList : private detail::Extendable<extensions...> { ...@@ -99,9 +99,9 @@ struct StepperExtensionList : private detail::Extendable<extensions...> {
/// returns a boolean as indicator if the evaluation is valid. /// returns a boolean as indicator if the evaluation is valid.
template <typename propagator_state_t, typename stepper_t> template <typename propagator_state_t, typename stepper_t>
bool k3(const propagator_state_t& state, const stepper_t& stepper, 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) { 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); kprev);
} }
...@@ -110,9 +110,9 @@ struct StepperExtensionList : private detail::Extendable<extensions...> { ...@@ -110,9 +110,9 @@ struct StepperExtensionList : private detail::Extendable<extensions...> {
/// returns a boolean as indicator if the evaluation is valid. /// returns a boolean as indicator if the evaluation is valid.
template <typename propagator_state_t, typename stepper_t> template <typename propagator_state_t, typename stepper_t>
bool k4(const propagator_state_t& state, const stepper_t& stepper, 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) { 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); kprev);
} }
......
...@@ -55,22 +55,22 @@ struct stepper_extension_list_impl { ...@@ -55,22 +55,22 @@ struct stepper_extension_list_impl {
template <typename propagator_state_t, typename stepper_t, typename... T> template <typename propagator_state_t, typename stepper_t, typename... T>
static bool k(std::tuple<T...>& obs_tuple, const propagator_state_t& state, static bool k(std::tuple<T...>& obs_tuple, const propagator_state_t& state,
const stepper_t& stepper, Vector3D& knew, 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 std::array<bool, sizeof...(T)>& validExtensions,
const int i = 0, const double h = 0, const int i = 0, const double h = 0,
const Vector3D& kprev = Vector3D()) { const Vector3D& kprev = Vector3D()) {
// If element is invalid: continue // If element is invalid: continue
if (!std::get<N - 1>(validExtensions)) { if (!std::get<N - 1>(validExtensions)) {
return stepper_extension_list_impl<N - 1>::k( 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); kprev);
} }
// Continue as long as evaluations are 'true' // 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)) { kprev)) {
return stepper_extension_list_impl<N - 1>::k( 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); kprev);
} else { } else {
// Break at false // Break at false
...@@ -142,7 +142,7 @@ struct stepper_extension_list_impl<0u> { ...@@ -142,7 +142,7 @@ struct stepper_extension_list_impl<0u> {
static bool k(std::tuple<T...>& /*unused*/, static bool k(std::tuple<T...>& /*unused*/,
const propagator_state_t& /*unused*/, const propagator_state_t& /*unused*/,
const stepper_t& /*unused*/, Vector3D& /*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 std::array<bool, sizeof...(T)>& /*unused*/,
const int /*unused*/, const double /*unused*/, const int /*unused*/, const double /*unused*/,
const Vector3D& /*unused*/) { const Vector3D& /*unused*/) {
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment