From b78e945db06feec6d8d20b21d3817d1833e420eb Mon Sep 17 00:00:00 2001 From: Andreas Stefl Date: Mon, 20 Jan 2025 18:21:03 +0100 Subject: [PATCH 01/11] refactor!: Decouple stepping from navigation --- Core/include/Acts/Propagator/AtlasStepper.hpp | 78 ++++---- Core/include/Acts/Propagator/EigenStepper.hpp | 28 ++- Core/include/Acts/Propagator/EigenStepper.ipp | 168 +++++++++--------- .../EigenStepperDefaultExtension.hpp | 86 +++++---- .../Propagator/EigenStepperDenseExtension.hpp | 128 ++++++------- .../Acts/Propagator/MultiEigenStepperLoop.hpp | 18 +- .../Acts/Propagator/MultiEigenStepperLoop.ipp | 43 ++--- Core/include/Acts/Propagator/Propagator.hpp | 2 +- Core/include/Acts/Propagator/Propagator.ipp | 24 +-- .../Acts/Propagator/PropagatorState.hpp | 2 + .../Acts/Propagator/PropagatorTraits.hpp | 9 +- .../Acts/Propagator/StraightLineStepper.hpp | 92 +++++----- Core/include/Acts/Propagator/SympyStepper.hpp | 31 ++-- Core/include/Acts/Propagator/SympyStepper.ipp | 22 --- .../include/Acts/Propagator/VoidNavigator.hpp | 10 ++ .../Propagator/detail/LoopStepperUtils.hpp | 8 + .../Acts/TrackFitting/KalmanFitter.hpp | 2 +- .../Acts/TrackFitting/detail/GsfActor.hpp | 10 +- .../detail/KalmanUpdateHelpers.hpp | 15 +- Core/src/Propagator/SympyStepper.cpp | 22 ++- .../Core/Propagator/AtlasStepperTests.cpp | 103 +++++------ .../Core/Propagator/EigenStepperTests.cpp | 109 +++++------- .../Core/Propagator/MultiStepperTests.cpp | 58 +----- .../Propagator/StraightLineStepperTests.cpp | 98 ++++------ .../Core/Propagator/SympyStepperTests.cpp | 95 ++++------ 25 files changed, 526 insertions(+), 735 deletions(-) delete mode 100644 Core/include/Acts/Propagator/SympyStepper.ipp diff --git a/Core/include/Acts/Propagator/AtlasStepper.hpp b/Core/include/Acts/Propagator/AtlasStepper.hpp index 4a103bcaae5..c698e17dd52 100644 --- a/Core/include/Acts/Propagator/AtlasStepper.hpp +++ b/Core/include/Acts/Propagator/AtlasStepper.hpp @@ -30,6 +30,8 @@ namespace Acts { +class IVolumeMaterial; + /// @brief the AtlasStepper implementation for the /// /// This is based original stepper code from the ATLAS RungeKuttaPropagator @@ -528,13 +530,10 @@ class AtlasStepper { /// 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 + /// @param [in, out] state The stepping state (thread-local cache) /// @return true if nothing is missing after this call, false otherwise. - template - bool prepareCurvilinearState( - [[maybe_unused]] propagator_state_t& prop_state, - [[maybe_unused]] const navigator_t& navigator) const { + bool prepareCurvilinearState(State& state) const { + (void)state; return true; } @@ -1154,32 +1153,33 @@ class AtlasStepper { /// Perform the actual step on the state /// /// @param state is the provided stepper state (caller keeps thread locality) - template - Result step(propagator_state_t& state, - const navigator_t& /*navigator*/) const { + Result step(State& state, Direction propDir, + const IVolumeMaterial* material) const { + (void)material; + // we use h for keeping the nominclature with the original atlas code - auto h = state.stepping.stepSize.value() * state.options.direction; - bool Jac = state.stepping.useJacobian; + auto h = state.stepSize.value() * propDir; + bool Jac = state.useJacobian; - double* R = &(state.stepping.pVector[0]); // Coordinates - double* A = &(state.stepping.pVector[4]); // Directions - double* sA = &(state.stepping.pVector[56]); + double* R = &(state.pVector[0]); // Coordinates + double* A = &(state.pVector[4]); // Directions + double* sA = &(state.pVector[56]); // Invert mometum/2. - double Pi = 0.5 * state.stepping.pVector[7]; + double Pi = 0.5 * state.pVector[7]; // double dltm = 0.0002 * .03; Vector3 f0, f; // if new field is required get it - if (state.stepping.newfield) { + if (state.newfield) { const Vector3 pos(R[0], R[1], R[2]); // This is sd.B_first in EigenStepper - auto fRes = getField(state.stepping, pos); + auto fRes = getField(state, pos); if (!fRes.ok()) { return fRes.error(); } f0 = *fRes; } else { - f0 = state.stepping.field; + f0 = state.field; } bool Helix = false; @@ -1215,7 +1215,7 @@ class AtlasStepper { // This is pos1 in EigenStepper const Vector3 pos(R[0] + A1 * S4, R[1] + B1 * S4, R[2] + C1 * S4); // This is sd.B_middle in EigenStepper - auto fRes = getField(state.stepping, pos); + auto fRes = getField(state, pos); if (!fRes.ok()) { return fRes.error(); } @@ -1245,7 +1245,7 @@ class AtlasStepper { // This is pos2 in EigenStepper const Vector3 pos(R[0] + h * A4, R[1] + h * B4, R[2] + h * C4); // This is sd.B_last in Eigen stepper - auto fRes = getField(state.stepping, pos); + auto fRes = getField(state, pos); if (!fRes.ok()) { return fRes.error(); } @@ -1269,10 +1269,10 @@ class AtlasStepper { 2. * h * (std::abs((A1 + A6) - (A3 + A4)) + std::abs((B1 + B6) - (B3 + B4)) + std::abs((C1 + C6) - (C3 + C4))); - if (std::abs(EST) > std::abs(state.options.stepping.stepTolerance)) { + if (std::abs(EST) > std::abs(state.options.stepTolerance)) { h = h * .5; // neutralize the sign of h again - state.stepping.stepSize.setAccuracy(h * state.options.direction); + state.stepSize.setAccuracy(h * propDir); // dltm = 0.; continue; } @@ -1301,25 +1301,25 @@ class AtlasStepper { sA[1] = B6 * Sl; sA[2] = C6 * Sl; - double mass = particleHypothesis(state.stepping).mass(); - double momentum = absoluteMomentum(state.stepping); + double mass = particleHypothesis(state).mass(); + double momentum = absoluteMomentum(state); // Evaluate the time propagation double dtds = std::sqrt(1 + mass * mass / (momentum * momentum)); - state.stepping.pVector[3] += h * dtds; - state.stepping.pVector[59] = dtds; - state.stepping.field = f; - state.stepping.newfield = false; + state.pVector[3] += h * dtds; + state.pVector[59] = dtds; + state.field = f; + state.newfield = false; if (Jac) { - double dtdl = h * mass * mass * qOverP(state.stepping) / dtds; - state.stepping.pVector[43] += dtdl; + double dtdl = h * mass * mass * qOverP(state) / dtds; + state.pVector[43] += dtdl; // Jacobian calculation // - double* d2A = &state.stepping.pVector[28]; - double* d3A = &state.stepping.pVector[36]; - double* d4A = &state.stepping.pVector[44]; + double* d2A = &state.pVector[28]; + double* d3A = &state.pVector[36]; + double* d4A = &state.pVector[44]; double d2A0 = H0[2] * d2A[1] - H0[1] * d2A[2]; double d2B0 = H0[0] * d2A[2] - H0[2] * d2A[0]; double d2C0 = H0[1] * d2A[0] - H0[0] * d2A[1]; @@ -1378,7 +1378,7 @@ class AtlasStepper { double d4B6 = d4C5 * H2[0] - d4A5 * H2[2]; double d4C6 = d4A5 * H2[1] - d4B5 * H2[0]; - double* dR = &state.stepping.pVector[24]; + double* dR = &state.pVector[24]; dR[0] += (d2A2 + d2A3 + d2A4) * S3; dR[1] += (d2B2 + d2B3 + d2B4) * S3; dR[2] += (d2C2 + d2C3 + d2C4) * S3; @@ -1386,7 +1386,7 @@ class AtlasStepper { d2A[1] = ((d2B0 + 2. * d2B3) + (d2B5 + d2B6)) * (1. / 3.); d2A[2] = ((d2C0 + 2. * d2C3) + (d2C5 + d2C6)) * (1. / 3.); - dR = &state.stepping.pVector[32]; + dR = &state.pVector[32]; dR[0] += (d3A2 + d3A3 + d3A4) * S3; dR[1] += (d3B2 + d3B3 + d3B4) * S3; dR[2] += (d3C2 + d3C3 + d3C4) * S3; @@ -1394,7 +1394,7 @@ class AtlasStepper { d3A[1] = ((d3B0 + 2. * d3B3) + (d3B5 + d3B6)) * (1. / 3.); d3A[2] = ((d3C0 + 2. * d3C3) + (d3C5 + d3C6)) * (1. / 3.); - dR = &state.stepping.pVector[40]; + dR = &state.pVector[40]; dR[0] += (d4A2 + d4A3 + d4A4) * S3; dR[1] += (d4B2 + d4B3 + d4B4) * S3; dR[2] += (d4C2 + d4C3 + d4C4) * S3; @@ -1406,9 +1406,9 @@ class AtlasStepper { break; } - state.stepping.pathAccumulated += h; - ++state.stepping.nSteps; - state.stepping.nStepTrials += nStepTrials; + state.pathAccumulated += h; + ++state.nSteps; + state.nStepTrials += nStepTrials; return h; } diff --git a/Core/include/Acts/Propagator/EigenStepper.hpp b/Core/include/Acts/Propagator/EigenStepper.hpp index 5dd61626f3c..e28f45d3063 100644 --- a/Core/include/Acts/Propagator/EigenStepper.hpp +++ b/Core/include/Acts/Propagator/EigenStepper.hpp @@ -30,6 +30,8 @@ namespace Acts { +class IVolumeMaterial; + /// @brief Runge-Kutta-Nystroem stepper based on Eigen implementation /// for the following ODE: /// @@ -236,18 +238,18 @@ class EigenStepper { /// @param [in,out] state The stepping state (thread-local cache) /// @param [in] surface The surface provided /// @param [in] index The surface intersection index - /// @param [in] navDir The navigation direction + /// @param [in] propDir The propagation direction /// @param [in] boundaryTolerance The boundary check for this status update /// @param [in] surfaceTolerance Surface tolerance used for intersection /// @param [in] stype The step size type to be set /// @param [in] logger A @c Logger instance IntersectionStatus updateSurfaceStatus( State& state, const Surface& surface, std::uint8_t index, - Direction navDir, const BoundaryTolerance& boundaryTolerance, + Direction propDir, const BoundaryTolerance& boundaryTolerance, double surfaceTolerance, ConstrainedStep::Type stype, const Logger& logger = getDummyLogger()) const { return detail::updateSingleSurfaceStatus( - *this, state, surface, index, navDir, boundaryTolerance, + *this, state, surface, index, propDir, boundaryTolerance, surfaceTolerance, stype, logger); } @@ -330,12 +332,9 @@ class EigenStepper { /// 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 + /// @param [in, out] state The state of the stepper /// @return true if nothing is missing after this call, false otherwise. - template - bool prepareCurvilinearState(propagator_state_t& prop_state, - const navigator_t& navigator) const; + bool prepareCurvilinearState(State& state) const; /// Create and return a curvilinear state at the current position /// @@ -397,15 +396,13 @@ class EigenStepper { /// Perform a Runge-Kutta track parameter propagation step /// - /// @param [in,out] state the propagation state - /// @param [in] navigator the navigator of the propagation + /// @param [in,out] state The state of the stepper /// @note The state contains the desired step size. It can be negative during /// backwards track propagation, and since we're using an adaptive /// algorithm, it can be modified by the stepper class during /// propagation. - template - Result step(propagator_state_t& state, - const navigator_t& navigator) const; + Result step(State& state, Direction propDir, + const IVolumeMaterial* material) const; /// Method that reset the Jacobian to the Identity for when no bound state are /// available @@ -418,9 +415,8 @@ class EigenStepper { std::shared_ptr m_bField; }; -template -struct SupportsBoundParameters, navigator_t> - : public std::true_type {}; +template <> +struct SupportsBoundParameters> : public std::true_type {}; } // namespace Acts diff --git a/Core/include/Acts/Propagator/EigenStepper.ipp b/Core/include/Acts/Propagator/EigenStepper.ipp index a9e4b25d02b..6f66ea66302 100644 --- a/Core/include/Acts/Propagator/EigenStepper.ipp +++ b/Core/include/Acts/Propagator/EigenStepper.ipp @@ -20,8 +20,9 @@ Acts::EigenStepper::EigenStepper( : m_bField(std::move(bField)) {} template -auto Acts::EigenStepper::makeState( - const Options& options, const BoundTrackParameters& par) const -> State { +auto Acts::EigenStepper::makeState(const Options& options, + const BoundTrackParameters& par) const + -> State { State state{options, m_bField->makeCache(options.magFieldContext)}; state.particleHypothesis = par.particleHypothesis(); @@ -35,7 +36,6 @@ auto Acts::EigenStepper::makeState( // Init the jacobian matrix if needed if (par.covariance()) { - // Get the reference surface for navigation const auto& surface = par.referenceSurface(); // set the covariance transport flag to true and copy state.covTransport = true; @@ -86,41 +86,43 @@ auto Acts::EigenStepper::boundState( } template -template -bool Acts::EigenStepper::prepareCurvilinearState( - propagator_state_t& prop_state, const navigator_t& navigator) const { +bool Acts::EigenStepper::prepareCurvilinearState(State& state) 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. - // 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.template k<0>( - 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; - } - } + if (state.pathAccumulated != 0) { + return true; + } + + // 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. + // First Runge-Kutta point (at current position) + auto& sd = state.stepData; + auto pos = position(state); + auto fieldRes = getField(state, pos); + if (!fieldRes.ok()) { return false; } + + sd.B_first = *fieldRes; + if (!state.extension.template k<0>(state, *this, nullptr, sd.k1, sd.B_first, + sd.kQoP)) { + return false; + } + + // dr/ds : + state.derivative.template head<3>() = + state.pars.template segment<3>(eFreeDir0); + // d (dr/ds) / ds : + state.derivative.template segment<3>(4) = sd.k1; + // to set dt/ds : + state.extension.finalize(state, *this, nullptr, state.pathAccumulated); return true; } template -auto Acts::EigenStepper::curvilinearState( - State& state, bool transportCov) const -> CurvilinearState { +auto Acts::EigenStepper::curvilinearState(State& state, + bool transportCov) const + -> CurvilinearState { return detail::curvilinearState( state.cov, state.jacobian, state.jacTransport, state.derivative, state.jacToGlobal, state.pars, state.particleHypothesis, @@ -168,27 +170,26 @@ void Acts::EigenStepper::transportCovarianceToBound( } template -template Acts::Result Acts::EigenStepper::step( - propagator_state_t& state, const navigator_t& navigator) const { + State& state, Direction propDir, const IVolumeMaterial* material) const { // Runge-Kutta integrator state - auto& sd = state.stepping.stepData; + auto& sd = state.stepData; double errorEstimate = 0; double h2 = 0; double half_h = 0; - auto pos = position(state.stepping); - auto dir = direction(state.stepping); + auto pos = position(state); + auto dir = direction(state); // First Runge-Kutta point (at current position) - auto fieldRes = getField(state.stepping, pos); + auto fieldRes = getField(state, pos); if (!fieldRes.ok()) { return fieldRes.error(); } sd.B_first = *fieldRes; - if (!state.stepping.extension.template k<0>(state, *this, navigator, sd.k1, - sd.B_first, sd.kQoP)) { + if (!state.extension.template k<0>(state, *this, material, sd.k1, sd.B_first, + sd.kQoP)) { return 0.; } @@ -199,7 +200,7 @@ Acts::Result Acts::EigenStepper::step( // This is given by the order of the Runge-Kutta method constexpr double exponent = 0.25; - double x = state.options.stepping.stepTolerance / errorEstimate_; + double x = state.options.stepTolerance / errorEstimate_; if constexpr (exponent == 0.25) { // This is 3x faster than std::pow @@ -215,8 +216,7 @@ Acts::Result Acts::EigenStepper::step( // For details about these values see ATL-SOFT-PUB-2009-001 constexpr double marginFactor = 4.0; - return errorEstimate_ <= - marginFactor * state.options.stepping.stepTolerance; + return errorEstimate_ <= marginFactor * state.options.stepTolerance; }; // The following functor starts to perform a Runge-Kutta step of a certain @@ -234,34 +234,32 @@ Acts::Result Acts::EigenStepper::step( // Second Runge-Kutta point const Vector3 pos1 = pos + half_h * dir + h2 * 0.125 * sd.k1; - auto field = getField(state.stepping, pos1); + auto field = getField(state, pos1); if (!field.ok()) { return failure(field.error()); } sd.B_middle = *field; - if (!state.stepping.extension.template k<1>(state, *this, navigator, sd.k2, - sd.B_middle, sd.kQoP, half_h, - sd.k1)) { + if (!state.extension.template k<1>(state, *this, material, sd.k2, + sd.B_middle, sd.kQoP, half_h, sd.k1)) { return success(false); } // Third Runge-Kutta point - if (!state.stepping.extension.template k<2>(state, *this, navigator, sd.k3, - sd.B_middle, sd.kQoP, half_h, - sd.k2)) { + if (!state.extension.template k<2>(state, *this, material, sd.k3, + sd.B_middle, sd.kQoP, half_h, sd.k2)) { return success(false); } // Last Runge-Kutta point const Vector3 pos2 = pos + h * dir + h2 * 0.5 * sd.k3; - field = getField(state.stepping, pos2); + field = getField(state, pos2); if (!field.ok()) { return failure(field.error()); } sd.B_last = *field; - if (!state.stepping.extension.template k<3>(state, *this, navigator, sd.k4, - sd.B_last, sd.kQoP, h, sd.k3)) { + if (!state.extension.template k<3>(state, *this, material, sd.k4, sd.B_last, + sd.kQoP, h, sd.k3)) { return success(false); } @@ -275,15 +273,14 @@ Acts::Result Acts::EigenStepper::step( return success(isErrorTolerable(errorEstimate)); }; - const double initialH = - state.stepping.stepSize.value() * state.options.direction; + const double initialH = state.stepSize.value() * propDir; double h = initialH; std::size_t nStepTrials = 0; // Select and adjust the appropriate Runge-Kutta step size as given // ATL-SOFT-PUB-2009-001 while (true) { ++nStepTrials; - ++state.stepping.statistics.nAttemptedSteps; + ++state.statistics.nAttemptedSteps; auto res = tryRungeKuttaStep(h); if (!res.ok()) { @@ -293,33 +290,33 @@ Acts::Result Acts::EigenStepper::step( break; } - ++state.stepping.statistics.nRejectedSteps; + ++state.statistics.nRejectedSteps; const double stepSizeScaling = calcStepSizeScaling(errorEstimate); h *= stepSizeScaling; // If step size becomes too small the particle remains at the initial // place - if (std::abs(h) < std::abs(state.options.stepping.stepSizeCutOff)) { + if (std::abs(h) < std::abs(state.options.stepSizeCutOff)) { // Not moving due to too low momentum needs an aborter return EigenStepperError::StepSizeStalled; } // If the parameter is off track too much or given stepSize is not // appropriate - if (nStepTrials > state.options.stepping.maxRungeKuttaStepTrials) { + if (nStepTrials > state.options.maxRungeKuttaStepTrials) { // Too many trials, have to abort return EigenStepperError::StepSizeAdjustmentFailed; } } // When doing error propagation, update the associated Jacobian matrix - if (state.stepping.covTransport) { + if (state.covTransport) { // using the direction before updated below // The step transport matrix in global coordinates FreeMatrix D; - if (!state.stepping.extension.finalize(state, *this, navigator, h, D)) { + if (!state.extension.finalize(state, *this, material, h, D)) { return EigenStepperError::StepInvalid; } @@ -342,56 +339,53 @@ Acts::Result Acts::EigenStepper::step( // sub-matrices at all! assert((D.topLeftCorner<4, 4>().isIdentity())); assert((D.bottomLeftCorner<4, 4>().isZero())); - assert((state.stepping.jacTransport.template topLeftCorner<4, 4>() - .isIdentity())); - assert((state.stepping.jacTransport.template bottomLeftCorner<4, 4>() - .isZero())); + assert((state.jacTransport.template topLeftCorner<4, 4>().isIdentity())); + assert((state.jacTransport.template bottomLeftCorner<4, 4>().isZero())); - state.stepping.jacTransport.template topRightCorner<4, 4>() += + state.jacTransport.template topRightCorner<4, 4>() += D.topRightCorner<4, 4>() * - state.stepping.jacTransport.template bottomRightCorner<4, 4>(); - state.stepping.jacTransport.template bottomRightCorner<4, 4>() = + state.jacTransport.template bottomRightCorner<4, 4>(); + state.jacTransport.template bottomRightCorner<4, 4>() = (D.bottomRightCorner<4, 4>() * - state.stepping.jacTransport.template bottomRightCorner<4, 4>()) + state.jacTransport.template bottomRightCorner<4, 4>()) .eval(); } else { - if (!state.stepping.extension.finalize(state, *this, navigator, h)) { + if (!state.extension.finalize(state, *this, material, h)) { return EigenStepperError::StepInvalid; } } // Update the track parameters according to the equations of motion - state.stepping.pars.template segment<3>(eFreePos0) += + state.pars.template segment<3>(eFreePos0) += h * dir + h2 / 6. * (sd.k1 + sd.k2 + sd.k3); - state.stepping.pars.template segment<3>(eFreeDir0) += + state.pars.template segment<3>(eFreeDir0) += h / 6. * (sd.k1 + 2. * (sd.k2 + sd.k3) + sd.k4); - (state.stepping.pars.template segment<3>(eFreeDir0)).normalize(); + (state.pars.template segment<3>(eFreeDir0)).normalize(); - if (state.stepping.covTransport) { + if (state.covTransport) { // using the updated direction - state.stepping.derivative.template head<3>() = - state.stepping.pars.template segment<3>(eFreeDir0); - state.stepping.derivative.template segment<3>(4) = sd.k4; + state.derivative.template head<3>() = + state.pars.template segment<3>(eFreeDir0); + state.derivative.template segment<3>(4) = sd.k4; } - state.stepping.pathAccumulated += h; - ++state.stepping.nSteps; - state.stepping.nStepTrials += nStepTrials; + state.pathAccumulated += h; + ++state.nSteps; + state.nStepTrials += nStepTrials; - ++state.stepping.statistics.nSuccessfulSteps; - if (state.options.direction != - Direction::fromScalarZeroAsPositive(initialH)) { - ++state.stepping.statistics.nReverseSteps; + ++state.statistics.nSuccessfulSteps; + if (propDir != Direction::fromScalarZeroAsPositive(initialH)) { + ++state.statistics.nReverseSteps; } - state.stepping.statistics.pathLength += h; - state.stepping.statistics.absolutePathLength += std::abs(h); + state.statistics.pathLength += h; + state.statistics.absolutePathLength += std::abs(h); const double stepSizeScaling = calcStepSizeScaling(errorEstimate); const double nextAccuracy = std::abs(h * stepSizeScaling); - const double previousAccuracy = std::abs(state.stepping.stepSize.accuracy()); + const double previousAccuracy = std::abs(state.stepSize.accuracy()); const double initialStepLength = std::abs(initialH); if (nextAccuracy < initialStepLength || nextAccuracy > previousAccuracy) { - state.stepping.stepSize.setAccuracy(nextAccuracy); + state.stepSize.setAccuracy(nextAccuracy); } return h; diff --git a/Core/include/Acts/Propagator/EigenStepperDefaultExtension.hpp b/Core/include/Acts/Propagator/EigenStepperDefaultExtension.hpp index c48bd9bbcc2..cfc69def7f5 100644 --- a/Core/include/Acts/Propagator/EigenStepperDefaultExtension.hpp +++ b/Core/include/Acts/Propagator/EigenStepperDefaultExtension.hpp @@ -16,6 +16,8 @@ namespace Acts { +class IVolumeMaterial; + /// @brief Default evaluator of the k_i's and elements of the transport matrix /// D of the RKN4 stepping. This is a pure implementation by textbook. struct EigenStepperDefaultExtension { @@ -23,11 +25,9 @@ struct EigenStepperDefaultExtension { /// step sets up qop, too. /// /// @tparam i Index of the k_i, i = [0, 3] - /// @tparam propagator_state_t Type of the state of the propagator /// @tparam stepper_t Type of the stepper - /// @tparam navigator_t Type of the navigator /// - /// @param [in] state State of the propagator + /// @param [in] state State of the stepper /// @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 @@ -36,22 +36,22 @@ struct EigenStepperDefaultExtension { /// @param [in] kprev Evaluated k_{i - 1} /// /// @return Boolean flag if the calculation is valid - template - bool k(const propagator_state_t& state, const stepper_t& stepper, - const navigator_t& /*navigator*/, Vector3& knew, const Vector3& bField, - std::array& kQoP, const double h = 0., - const Vector3& kprev = Vector3::Zero()) + template + bool k(const stepper_t::State& state, const stepper_t& stepper, + const IVolumeMaterial* volumeMaterial, Vector3& knew, + const Vector3& bField, std::array& kQoP, + const double h = 0., const Vector3& kprev = Vector3::Zero()) requires(i >= 0 && i <= 3) { - auto qop = stepper.qOverP(state.stepping); + (void)volumeMaterial; + + auto qop = stepper.qOverP(state); // First step does not rely on previous data if constexpr (i == 0) { - knew = qop * stepper.direction(state.stepping).cross(bField); + knew = qop * stepper.direction(state).cross(bField); kQoP = {0., 0., 0., 0.}; } else { - knew = - qop * (stepper.direction(state.stepping) + h * kprev).cross(bField); + knew = qop * (stepper.direction(state) + h * kprev).cross(bField); } return true; } @@ -60,19 +60,18 @@ struct EigenStepperDefaultExtension { /// error of the step. Since the textbook does not deliver further vetos, /// this is a dummy function. /// - /// @tparam propagator_state_t Type of the state of the propagator /// @tparam stepper_t Type of the stepper - /// @tparam navigator_t Type of the navigator /// - /// @param [in] state State of the propagator + /// @param [in] state State of the stepper /// @param [in] stepper Stepper of the propagation /// @param [in] h Step size /// /// @return Boolean flag if the calculation is valid - template - bool finalize(propagator_state_t& state, const stepper_t& stepper, - const navigator_t& /*navigator*/, const double h) const { + template + bool finalize(stepper_t::State& state, const stepper_t& stepper, + const IVolumeMaterial* volumeMaterial, const double h) const { + (void)volumeMaterial; + propagateTime(state, stepper, h); return true; } @@ -81,21 +80,20 @@ struct EigenStepperDefaultExtension { /// error of the step. Since the textbook does not deliver further vetos, /// this is just for the evaluation of the transport matrix. /// - /// @tparam propagator_state_t Type of the state of the propagator /// @tparam stepper_t Type of the stepper - /// @tparam navigator_t Type of the navigator /// - /// @param [in] state State of the propagator + /// @param [in] state State of the stepper /// @param [in] stepper Stepper of the propagation /// @param [in] h Step size /// @param [out] D Transport matrix /// /// @return Boolean flag if the calculation is valid - template - bool finalize(propagator_state_t& state, const stepper_t& stepper, - const navigator_t& /*navigator*/, const double h, + template + bool finalize(stepper_t::State state, const stepper_t& stepper, + const IVolumeMaterial* volumeMaterial, const double h, FreeMatrix& D) const { + (void)volumeMaterial; + propagateTime(state, stepper, h); return transportMatrix(state, stepper, h, D); } @@ -103,40 +101,38 @@ struct EigenStepperDefaultExtension { private: /// @brief Propagation function for the time coordinate /// - /// @tparam propagator_state_t Type of the state of the propagator /// @tparam stepper_t Type of the stepper /// - /// @param [in, out] state State of the propagator + /// @param [in, out] state State of the stepper /// @param [in] stepper Stepper of the propagation /// @param [in] h Step size - template - void propagateTime(propagator_state_t& state, const stepper_t& stepper, + template + void propagateTime(stepper_t::State& state, const stepper_t& stepper, const double h) const { /// This evaluation is based on dt/ds = 1/v = 1/(beta * c) with the velocity /// v, the speed of light c and beta = v/c. This can be re-written as dt/ds /// = sqrt(m^2/p^2 + c^{-2}) with the mass m and the momentum p. - auto m = stepper.particleHypothesis(state.stepping).mass(); - auto p = stepper.absoluteMomentum(state.stepping); + auto m = stepper.particleHypothesis(state).mass(); + auto p = stepper.absoluteMomentum(state); auto dtds = std::sqrt(1 + m * m / (p * p)); - state.stepping.pars[eFreeTime] += h * dtds; - if (state.stepping.covTransport) { - state.stepping.derivative(3) = dtds; + state.pars[eFreeTime] += h * dtds; + if (state.covTransport) { + state.derivative(3) = dtds; } } /// @brief Calculates the transport matrix D for the jacobian /// - /// @tparam propagator_state_t Type of the state of the propagator /// @tparam stepper_t Type of the stepper /// - /// @param [in] state State of the propagator + /// @param [in] state State of the stepper /// @param [in] stepper Stepper of the propagation /// @param [in] h Step size /// @param [out] D Transport matrix /// /// @return Boolean flag if evaluation is valid - template - bool transportMatrix(propagator_state_t& state, const stepper_t& stepper, + template + bool transportMatrix(stepper_t::State& state, const stepper_t& stepper, const double h, FreeMatrix& D) const { /// The calculations are based on ATL-SOFT-PUB-2009-002. The update of the /// Jacobian matrix is requires only the calculation of eq. 17 and 18. @@ -157,11 +153,11 @@ struct EigenStepperDefaultExtension { /// constant offset does not exist for rectangular matrix dGdu' (due to the /// missing Lambda part) and only exists for dFdu' in dlambda/dlambda. - auto m = state.stepping.particleHypothesis.mass(); - auto& sd = state.stepping.stepData; - auto dir = stepper.direction(state.stepping); - auto qop = stepper.qOverP(state.stepping); - auto p = stepper.absoluteMomentum(state.stepping); + auto m = state.particleHypothesis.mass(); + auto& sd = state.stepData; + auto dir = stepper.direction(state); + auto qop = stepper.qOverP(state); + auto p = stepper.absoluteMomentum(state); auto dtds = std::sqrt(1 + m * m / (p * p)); D = FreeMatrix::Identity(); diff --git a/Core/include/Acts/Propagator/EigenStepperDenseExtension.hpp b/Core/include/Acts/Propagator/EigenStepperDenseExtension.hpp index 935e747f3aa..1aafd3967f9 100644 --- a/Core/include/Acts/Propagator/EigenStepperDenseExtension.hpp +++ b/Core/include/Acts/Propagator/EigenStepperDenseExtension.hpp @@ -12,11 +12,11 @@ #include "Acts/Utilities/detail/ReferenceWrapperAnyCompat.hpp" #include "Acts/Definitions/Algebra.hpp" +#include "Acts/Material/IVolumeMaterial.hpp" #include "Acts/Material/Interactions.hpp" #include "Acts/Material/Material.hpp" #include "Acts/Material/MaterialSlab.hpp" #include "Acts/Propagator/EigenStepperDefaultExtension.hpp" -#include "Acts/Propagator/Propagator.hpp" #include "Acts/Utilities/MathHelpers.hpp" namespace Acts { @@ -58,13 +58,10 @@ struct EigenStepperDenseExtension { /// step sets up member parameters, too. /// /// @tparam i Index of the k_i, i = [0, 3] - /// @tparam propagator_state_t Type of the state of the propagator /// @tparam stepper_t Type of the stepper - /// @tparam navigator_t Type of the navigator /// - /// @param [in] state State of the propagator + /// @param [in] state State of the stepper /// @param [in] stepper Stepper of the propagator - /// @param [in] navigator Navigator 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 @@ -72,36 +69,33 @@ struct EigenStepperDenseExtension { /// @param [in] kprev Evaluated k_{i - 1} /// /// @return Boolean flag if the calculation is valid - template - bool k(const propagator_state_t& state, const stepper_t& stepper, - const navigator_t& navigator, Vector3& knew, const Vector3& bField, - std::array& kQoP, const double h = 0., - const Vector3& kprev = Vector3::Zero()) + template + bool k(const stepper_t::State& state, const stepper_t& stepper, + const IVolumeMaterial* volumeMaterial, Vector3& knew, + const Vector3& bField, std::array& kQoP, + const double h = 0., const Vector3& kprev = Vector3::Zero()) requires(i >= 0 && i <= 3) { - const auto* volumeMaterial = - navigator.currentVolumeMaterial(state.navigation); if (volumeMaterial == nullptr) { - return defaultExtension.template k(state, stepper, navigator, knew, - bField, kQoP, h, kprev); + return defaultExtension.template k(state, stepper, volumeMaterial, + knew, bField, kQoP, h, kprev); } - double q = stepper.charge(state.stepping); - const auto& particleHypothesis = stepper.particleHypothesis(state.stepping); + double q = stepper.charge(state); + const auto& particleHypothesis = stepper.particleHypothesis(state); float mass = particleHypothesis.mass(); // i = 0 is used for setup and evaluation of k if constexpr (i == 0) { // Set up for energy loss - Vector3 position = stepper.position(state.stepping); + Vector3 position = stepper.position(state); material = volumeMaterial->material(position.template cast()); - initialMomentum = stepper.absoluteMomentum(state.stepping); + initialMomentum = stepper.absoluteMomentum(state); currentMomentum = initialMomentum; - qop[0] = stepper.qOverP(state.stepping); + qop[0] = stepper.qOverP(state); initializeEnergyLoss(state, stepper); // Evaluate k - knew = qop[0] * stepper.direction(state.stepping).cross(bField); + knew = qop[0] * stepper.direction(state).cross(bField); // Evaluate k for the time propagation Lambdappi[0] = -qop[0] * qop[0] * qop[0] * g * energy[0] / (q * q); //~ tKi[0] = std::hypot(1, mass / initialMomentum); @@ -110,12 +104,11 @@ struct EigenStepperDenseExtension { } else { // Update parameters and check for momentum condition updateEnergyLoss(mass, h, state, stepper, i); - if (currentMomentum < state.options.stepping.dense.momentumCutOff) { + if (currentMomentum < state.options.dense.momentumCutOff) { return false; } // Evaluate k - knew = qop[i] * - (stepper.direction(state.stepping) + h * kprev).cross(bField); + knew = qop[i] * (stepper.direction(state) + h * kprev).cross(bField); // Evaluate k_i for the time propagation auto qopNew = qop[0] + h * Lambdappi[i - 1]; Lambdappi[i] = -qopNew * qopNew * qopNew * g * energy[i] / (q * q); @@ -130,50 +123,43 @@ struct EigenStepperDenseExtension { /// of the energy loss and the therewith constrained to keep the momentum /// after the step in reasonable values. /// - /// @tparam propagator_state_t Type of the state of the propagator /// @tparam stepper_t Type of the stepper - /// @tparam navigator_t Type of the navigator /// - /// @param [in] state State of the propagator + /// @param [in] state State of the stepper /// @param [in] stepper Stepper of the propagator - /// @param [in] navigator Navigator of the propagator /// @param [in] h Step size /// /// @return Boolean flag if the calculation is valid - template - bool finalize(propagator_state_t& state, const stepper_t& stepper, - const navigator_t& navigator, const double h) const { - const auto* volumeMaterial = - navigator.currentVolumeMaterial(state.navigation); + template + bool finalize(stepper_t::State& state, const stepper_t& stepper, + const IVolumeMaterial* volumeMaterial, const double h) const { if (volumeMaterial == nullptr) { - return defaultExtension.finalize(state, stepper, navigator, h); + return defaultExtension.finalize(state, stepper, volumeMaterial, h); } - const auto& particleHypothesis = stepper.particleHypothesis(state.stepping); + const auto& particleHypothesis = stepper.particleHypothesis(state); float mass = particleHypothesis.mass(); // Evaluate the new momentum auto newMomentum = - stepper.absoluteMomentum(state.stepping) + + stepper.absoluteMomentum(state) + (h / 6.) * (dPds[0] + 2. * (dPds[1] + dPds[2]) + dPds[3]); // Break propagation if momentum becomes below cut-off - if (newMomentum < state.options.stepping.dense.momentumCutOff) { + if (newMomentum < state.options.dense.momentumCutOff) { return false; } // Add derivative dlambda/ds = Lambda'' - state.stepping.derivative(7) = -fastHypot(mass, newMomentum) * g / - (newMomentum * newMomentum * newMomentum); + state.derivative(7) = -fastHypot(mass, newMomentum) * g / + (newMomentum * newMomentum * newMomentum); // Update momentum - state.stepping.pars[eFreeQOverP] = - stepper.charge(state.stepping) / newMomentum; + state.pars[eFreeQOverP] = stepper.charge(state) / newMomentum; // Add derivative dt/ds = 1/(beta * c) = sqrt(m^2 * p^{-2} + c^{-2}) - state.stepping.derivative(3) = fastHypot(1, mass / newMomentum); + state.derivative(3) = fastHypot(1, mass / newMomentum); // Update time - state.stepping.pars[eFreeTime] += + state.pars[eFreeTime] += (h / 6.) * (tKi[0] + 2. * (tKi[1] + tKi[2]) + tKi[3]); return true; @@ -186,29 +172,23 @@ struct EigenStepperDenseExtension { /// after the step in reasonable values and the evaluation of the transport /// matrix. /// - /// @tparam propagator_state_t Type of the state of the propagator /// @tparam stepper_t Type of the stepper - /// @tparam navigator_t Type of the navigator /// - /// @param [in] state State of the propagator + /// @param [in] state State of the stepper /// @param [in] stepper Stepper of the propagator - /// @param [in] navigator Navigator of the propagator /// @param [in] h Step size /// @param [out] D Transport matrix /// /// @return Boolean flag if the calculation is valid - template - bool finalize(propagator_state_t& state, const stepper_t& stepper, - const navigator_t& navigator, const double h, + template + bool finalize(stepper_t::State& state, const stepper_t& stepper, + const IVolumeMaterial* volumeMaterial, const double h, FreeMatrix& D) const { - const auto* volumeMaterial = - navigator.currentVolumeMaterial(state.navigation); if (volumeMaterial == nullptr) { - return defaultExtension.finalize(state, stepper, navigator, h, D); + return defaultExtension.finalize(state, stepper, volumeMaterial, h, D); } - return finalize(state, stepper, navigator, h) && + return finalize(state, stepper, volumeMaterial, h) && transportMatrix(state, stepper, h, D); } @@ -224,8 +204,8 @@ struct EigenStepperDenseExtension { /// @param [out] D Transport matrix /// /// @return Boolean flag if evaluation is valid - template - bool transportMatrix(propagator_state_t& state, const stepper_t& stepper, + template + bool transportMatrix(stepper_t::State& state, const stepper_t& stepper, const double h, FreeMatrix& D) const { /// The calculations are based on ATL-SOFT-PUB-2009-002. The update of the /// Jacobian matrix is requires only the calculation of eq. 17 and 18. @@ -247,9 +227,9 @@ struct EigenStepperDenseExtension { /// constant offset does not exist for rectangular matrix dFdu' (due to the /// missing Lambda part) and only exists for dGdu' in dlambda/dlambda. - auto& sd = state.stepping.stepData; - auto dir = stepper.direction(state.stepping); - const auto& particleHypothesis = stepper.particleHypothesis(state.stepping); + auto& sd = state.stepData; + auto dir = stepper.direction(state); + const auto& particleHypothesis = stepper.particleHypothesis(state); float mass = particleHypothesis.mass(); D = FreeMatrix::Identity(); @@ -360,15 +340,14 @@ struct EigenStepperDenseExtension { /// @brief Initializer of all parameters related to a RKN4 step with energy /// loss of a particle in material /// - /// @tparam propagator_state_t Type of the state of the propagator /// @tparam stepper_t Type of the stepper /// /// @param [in] state Deliverer of configurations /// @param [in] stepper Stepper of the propagator - template - void initializeEnergyLoss(const propagator_state_t& state, + template + void initializeEnergyLoss(const stepper_t::State& state, const stepper_t& stepper) { - const auto& particleHypothesis = stepper.particleHypothesis(state.stepping); + const auto& particleHypothesis = stepper.particleHypothesis(state); float mass = particleHypothesis.mass(); PdgParticle absPdg = particleHypothesis.absolutePdg(); float absQ = particleHypothesis.absoluteCharge(); @@ -377,7 +356,7 @@ struct EigenStepperDenseExtension { // use unit length as thickness to compute the energy loss per unit length MaterialSlab slab(material, 1); // Use the same energy loss throughout the step. - if (state.options.stepping.dense.meanEnergyLoss) { + if (state.options.dense.meanEnergyLoss) { g = -computeEnergyLossMean(slab, absPdg, mass, static_cast(qop[0]), absQ); } else { @@ -389,11 +368,11 @@ struct EigenStepperDenseExtension { // Change of the momentum per path length // dPds = dPdE * dEds dPds[0] = g * energy[0] / initialMomentum; - if (state.stepping.covTransport) { + if (state.covTransport) { // Calculate the change of the energy loss per path length and // inverse momentum - if (state.options.stepping.dense.includeGradient) { - if (state.options.stepping.dense.meanEnergyLoss) { + if (state.options.dense.includeGradient) { + if (state.options.dense.meanEnergyLoss) { dgdqopValue = deriveEnergyLossMeanQOverP( slab, absPdg, mass, static_cast(qop[0]), absQ); } else { @@ -413,7 +392,6 @@ struct EigenStepperDenseExtension { /// @brief Update of the kinematic parameters of the RKN4 sub-steps after /// initialization with energy loss of a particle in material /// - /// @tparam propagator_state_t Type of the state of the propagator /// @tparam stepper_t Type of the stepper /// /// @param [in] h Stepped distance of the sub-step (1-3) @@ -421,17 +399,17 @@ struct EigenStepperDenseExtension { /// @param [in] state State of the stepper /// @param [in] stepper Stepper of the propagator /// @param [in] i Index of the sub-step (1-3) - template + template void updateEnergyLoss(const double mass, const double h, - const propagator_state_t& state, - const stepper_t& stepper, const int i) { + const stepper_t::State& state, const stepper_t& stepper, + const int i) { // Update parameters related to a changed momentum currentMomentum = initialMomentum + h * dPds[i - 1]; energy[i] = fastHypot(currentMomentum, mass); dPds[i] = g * energy[i] / currentMomentum; - qop[i] = stepper.charge(state.stepping) / currentMomentum; + qop[i] = stepper.charge(state) / currentMomentum; // Calculate term for later error propagation - if (state.stepping.covTransport) { + if (state.covTransport) { dLdl[i] = (-qop[i] * qop[i] * g * energy[i] * (3. - (currentMomentum * currentMomentum) / (energy[i] * energy[i])) - diff --git a/Core/include/Acts/Propagator/MultiEigenStepperLoop.hpp b/Core/include/Acts/Propagator/MultiEigenStepperLoop.hpp index 6e1364a16e4..a05e8281e33 100644 --- a/Core/include/Acts/Propagator/MultiEigenStepperLoop.hpp +++ b/Core/include/Acts/Propagator/MultiEigenStepperLoop.hpp @@ -675,13 +675,10 @@ class MultiEigenStepperLoop : public EigenStepper { /// 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 + /// @param [in, out] state The stepping state (thread-local cache) /// @return true if nothing is missing after this call, false otherwise. - template - bool prepareCurvilinearState( - [[maybe_unused]] propagator_state_t& prop_state, - [[maybe_unused]] const navigator_t& navigator) const { + bool prepareCurvilinearState(State& state) const { + (void)state; return true; } @@ -736,16 +733,13 @@ class MultiEigenStepperLoop : public EigenStepper { /// Perform a Runge-Kutta track parameter propagation step /// - /// @param [in,out] state is the propagation state associated with the track - /// parameters that are being propagated. - /// @param [in] navigator is the navigator of the propagation + /// @param [in,out] state The state of the stepper /// /// The state contains the desired step size. It can be negative during /// backwards track propagation, and since we're using an adaptive /// algorithm, it can be modified by the stepper class during propagation. - template - Result step(propagator_state_t& state, - const navigator_t& navigator) const; + Result step(State& state, Direction propDir, + const IVolumeMaterial* material) const; }; } // namespace Acts diff --git a/Core/include/Acts/Propagator/MultiEigenStepperLoop.ipp b/Core/include/Acts/Propagator/MultiEigenStepperLoop.ipp index 01b392dbf2a..c3153b2f278 100644 --- a/Core/include/Acts/Propagator/MultiEigenStepperLoop.ipp +++ b/Core/include/Acts/Propagator/MultiEigenStepperLoop.ipp @@ -64,8 +64,9 @@ auto MultiEigenStepperLoop::boundState( } template -auto MultiEigenStepperLoop::curvilinearState( - State& state, bool transportCov) const -> CurvilinearState { +auto MultiEigenStepperLoop::curvilinearState(State& state, + bool transportCov) const + -> CurvilinearState { assert(!state.components.empty()); std::vector> @@ -90,25 +91,23 @@ auto MultiEigenStepperLoop::curvilinearState( } template -template Result MultiEigenStepperLoop::step( - propagator_state_t& state, const navigator_t& navigator) const { + State& state, Direction propDir, const IVolumeMaterial* material) const { using Status = Acts::IntersectionStatus; - State& stepping = state.stepping; - auto& components = stepping.components; + auto& components = state.components; const Logger& logger = *m_logger; // Update step count - stepping.steps++; + state.steps++; // Check if we abort because of m_stepLimitAfterFirstComponentOnSurface - if (stepping.stepCounterAfterFirstComponentOnSurface) { - (*stepping.stepCounterAfterFirstComponentOnSurface)++; + if (state.stepCounterAfterFirstComponentOnSurface) { + (*state.stepCounterAfterFirstComponentOnSurface)++; // If the limit is reached, remove all components which are not on a // surface, reweight the components, perform no step and return 0 - if (*stepping.stepCounterAfterFirstComponentOnSurface >= + if (*state.stepCounterAfterFirstComponentOnSurface >= m_stepLimitAfterFirstComponentOnSurface) { for (auto& cmp : components) { if (cmp.status != Status::onSurface) { @@ -122,13 +121,13 @@ Result MultiEigenStepperLoop::step( ACTS_VERBOSE( "-> remove all components not on a surface, perform no step"); - removeMissedComponents(stepping); - reweightComponents(stepping); + removeMissedComponents(state); + reweightComponents(state); ACTS_VERBOSE(components.size() << " components left after removing missed components"); - stepping.stepCounterAfterFirstComponentOnSurface.reset(); + state.stepCounterAfterFirstComponentOnSurface.reset(); return 0.0; } @@ -145,7 +144,7 @@ Result MultiEigenStepperLoop::step( [&](auto& cmp) { return cmp.status == IntersectionStatus::onSurface; }); if (cmpsOnSurface > 0) { - removeMissedComponents(stepping); + removeMissedComponents(state); reweightNecessary = true; } @@ -155,12 +154,6 @@ Result MultiEigenStepperLoop::step( double accumulatedPathLength = 0.0; std::size_t errorSteps = 0; - // Type of the proxy single propagation2 state - using ThisSinglePropState = - detail::SinglePropState; - // Lambda that performs the step for a component and returns false if the step // went ok and true if there was an error auto errorInStep = [&](auto& component) { @@ -171,10 +164,8 @@ Result MultiEigenStepperLoop::step( return false; } - ThisSinglePropState single_state(component.state, state.navigation, - state.options, state.geoContext); - - results.emplace_back(SingleStepper::step(single_state, navigator)); + results.emplace_back( + SingleStepper::step(component.state, propDir, material)); if (results.back()->ok()) { accumulatedPathLength += component.weight * results.back()->value(); @@ -193,7 +184,7 @@ Result MultiEigenStepperLoop::step( // Reweight if necessary if (reweightNecessary) { - reweightComponents(stepping); + reweightComponents(state); } // Print the result vector to a string so we can log it @@ -231,7 +222,7 @@ Result MultiEigenStepperLoop::step( } // Return the weighted accumulated path length of all successful steps - stepping.pathAccumulated += accumulatedPathLength; + state.pathAccumulated += accumulatedPathLength; return accumulatedPathLength; } diff --git a/Core/include/Acts/Propagator/Propagator.hpp b/Core/include/Acts/Propagator/Propagator.hpp index d3664985fbc..37a53ec23a7 100644 --- a/Core/include/Acts/Propagator/Propagator.hpp +++ b/Core/include/Acts/Propagator/Propagator.hpp @@ -92,7 +92,7 @@ class BasePropagatorHelper : public BasePropagator { template class Propagator final : public std::conditional_t< - SupportsBoundParameters_v, + SupportsBoundParameters_v, detail::BasePropagatorHelper>, detail::PropagatorStub> { /// Re-define bound track parameters dependent on the stepper diff --git a/Core/include/Acts/Propagator/Propagator.ipp b/Core/include/Acts/Propagator/Propagator.ipp index 0790a5d906b..30d6442da44 100644 --- a/Core/include/Acts/Propagator/Propagator.ipp +++ b/Core/include/Acts/Propagator/Propagator.ipp @@ -18,14 +18,6 @@ #include -namespace Acts::detail { -template -concept propagator_stepper_compatible_with = - requires(const Stepper& s, StateType& st, const N& n) { - { s.step(st, n) } -> std::same_as>; - }; -} // namespace Acts::detail - template template auto Acts::Propagator::propagate(propagator_state_t& state) const @@ -89,7 +81,9 @@ auto Acts::Propagator::propagate(propagator_state_t& state) const // Stepping loop for (; state.steps < state.options.maxSteps; ++state.steps) { // Perform a step - Result res = m_stepper.step(state, m_navigator); + Result res = + m_stepper.step(state.stepping, state.options.direction, + m_navigator.currentVolumeMaterial(state.navigation)); if (!res.ok()) { ACTS_ERROR("Step failed with " << res.error() << ": " << res.error().message()); @@ -258,11 +252,6 @@ auto Acts::Propagator::makeState( StateType state{eOptions, m_stepper.makeState(eOptions.stepping, start), m_navigator.makeState(eOptions.navigation)}; - static_assert( - detail::propagator_stepper_compatible_with, - "Step method of the Stepper is not compatible with the propagator " - "state"); - initialize(state); return state; @@ -297,11 +286,6 @@ auto Acts::Propagator::makeState( StateType state{eOptions, m_stepper.makeState(eOptions.stepping, start), m_navigator.makeState(eOptions.navigation)}; - static_assert( - detail::propagator_stepper_compatible_with, - "Step method of the Stepper is not compatible with the propagator " - "state"); - initialize(state); return state; @@ -335,7 +319,7 @@ auto Acts::Propagator::makeResult(propagator_state_t state, moveStateToResult(state, result); if (makeCurvilinear) { - if (!m_stepper.prepareCurvilinearState(state, m_navigator)) { + if (!m_stepper.prepareCurvilinearState(state.stepping)) { // information to compute curvilinearState is incomplete. return propagationResult.error(); } diff --git a/Core/include/Acts/Propagator/PropagatorState.hpp b/Core/include/Acts/Propagator/PropagatorState.hpp index f8b639d1faa..e7abb7fc6ac 100644 --- a/Core/include/Acts/Propagator/PropagatorState.hpp +++ b/Core/include/Acts/Propagator/PropagatorState.hpp @@ -36,6 +36,8 @@ template struct PropagatorState : private detail::Extendable { using options_type = propagator_options_t; + using stepper_state_type = stepper_state_t; + using navigator_state_type = navigator_state_t; /// Create the propagator state from the options /// diff --git a/Core/include/Acts/Propagator/PropagatorTraits.hpp b/Core/include/Acts/Propagator/PropagatorTraits.hpp index 2bd1567fdb2..a19cd9c6617 100644 --- a/Core/include/Acts/Propagator/PropagatorTraits.hpp +++ b/Core/include/Acts/Propagator/PropagatorTraits.hpp @@ -9,11 +9,14 @@ #pragma once #include + namespace Acts { -template + +template struct SupportsBoundParameters : public std::false_type {}; -template +template constexpr bool SupportsBoundParameters_v = - SupportsBoundParameters::value; + SupportsBoundParameters::value; + } // namespace Acts diff --git a/Core/include/Acts/Propagator/StraightLineStepper.hpp b/Core/include/Acts/Propagator/StraightLineStepper.hpp index de472497090..d1e96f30e6b 100644 --- a/Core/include/Acts/Propagator/StraightLineStepper.hpp +++ b/Core/include/Acts/Propagator/StraightLineStepper.hpp @@ -36,6 +36,8 @@ namespace Acts { +class IVolumeMaterial; + /// @brief straight line stepper based on Surface intersection /// /// The straight line stepper is a simple navigation stepper @@ -315,28 +317,24 @@ class StraightLineStepper { /// 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 + /// @param [in, out] state The stepping state (thread-local cache) /// @return true if nothing is missing after this call, false otherwise. - template - bool prepareCurvilinearState( - [[maybe_unused]] propagator_state_t& prop_state, - [[maybe_unused]] const navigator_t& navigator) const { + bool prepareCurvilinearState(State& state) 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) = - fastHypot(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.; + if (state.pathAccumulated != 0) { + return true; } + + // dr/ds : + state.derivative.template head<3>() = direction(state); + // dt / ds + state.derivative(eFreeTime) = fastHypot( + 1., state.particleHypothesis.mass() / absoluteMomentum(state)); + // d (dr/ds) / ds : == 0 + state.derivative.template segment<3>(4) = Acts::Vector3::Zero().transpose(); + // d qop / ds == 0 + state.derivative(eFreeQOverP) = 0.; + return true; } @@ -401,7 +399,7 @@ class StraightLineStepper { /// Perform a straight line propagation step /// - /// @param [in,out] state is the propagation state associated with the track + /// @param [in,out] state The state of the stepper /// parameters that are being propagated. /// The state contains the desired step size, /// it can be negative during backwards track propagation, @@ -409,54 +407,54 @@ class StraightLineStepper { /// be modified by the stepper class during propagation. /// /// @return the step size taken - template - Result step(propagator_state_t& state, - const navigator_t& /*navigator*/) const { + Result step(State& state, Direction propDir, + const IVolumeMaterial* material) const { + (void)material; + // use the adjusted step size - const auto h = state.stepping.stepSize.value() * state.options.direction; - const auto m = state.stepping.particleHypothesis.mass(); - const auto p = absoluteMomentum(state.stepping); + const auto h = state.stepSize.value() * propDir; + const auto m = state.particleHypothesis.mass(); + const auto p = absoluteMomentum(state); // time propagates along distance as 1/b = sqrt(1 + m²/p²) const auto dtds = fastHypot(1., m / p); // Update the track parameters according to the equations of motion - Vector3 dir = direction(state.stepping); - state.stepping.pars.template segment<3>(eFreePos0) += h * dir; - state.stepping.pars[eFreeTime] += h * dtds; + Vector3 dir = direction(state); + state.pars.template segment<3>(eFreePos0) += h * dir; + state.pars[eFreeTime] += h * dtds; // Propagate the jacobian - if (state.stepping.covTransport) { + if (state.covTransport) { // The step transport matrix in global coordinates FreeMatrix D = FreeMatrix::Identity(); D.block<3, 3>(0, 4) = ActsSquareMatrix<3>::Identity() * h; // Extend the calculation by the time propagation // Evaluate dt/dlambda - D(3, 7) = h * m * m * state.stepping.pars[eFreeQOverP] / dtds; + D(3, 7) = h * m * m * state.pars[eFreeQOverP] / dtds; // Set the derivative factor the time - state.stepping.derivative(3) = dtds; + state.derivative(3) = dtds; // Update jacobian and derivative - state.stepping.jacTransport = D * state.stepping.jacTransport; - state.stepping.derivative.template head<3>() = dir; + state.jacTransport = D * state.jacTransport; + state.derivative.template head<3>() = dir; } // state the path length - state.stepping.pathAccumulated += h; - ++state.stepping.nSteps; - ++state.stepping.nStepTrials; - - ++state.stepping.statistics.nAttemptedSteps; - ++state.stepping.statistics.nSuccessfulSteps; - if (state.options.direction != Direction::fromScalarZeroAsPositive(h)) { - ++state.stepping.statistics.nReverseSteps; + state.pathAccumulated += h; + ++state.nSteps; + ++state.nStepTrials; + + ++state.statistics.nAttemptedSteps; + ++state.statistics.nSuccessfulSteps; + if (propDir != Direction::fromScalarZeroAsPositive(h)) { + ++state.statistics.nReverseSteps; } - state.stepping.statistics.pathLength += h; - state.stepping.statistics.absolutePathLength += std::abs(h); + state.statistics.pathLength += h; + state.statistics.absolutePathLength += std::abs(h); return h; } }; -template -struct SupportsBoundParameters - : public std::true_type {}; +template <> +struct SupportsBoundParameters : public std::true_type {}; } // namespace Acts diff --git a/Core/include/Acts/Propagator/SympyStepper.hpp b/Core/include/Acts/Propagator/SympyStepper.hpp index 53b3474b851..8fa694d6176 100644 --- a/Core/include/Acts/Propagator/SympyStepper.hpp +++ b/Core/include/Acts/Propagator/SympyStepper.hpp @@ -24,6 +24,8 @@ namespace Acts { +class IVolumeMaterial; + class SympyStepper { public: /// Jacobian, Covariance and State definitions @@ -297,15 +299,9 @@ class SympyStepper { /// 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 + /// @param [in, out] state State of the stepper /// @return true if nothing is missing after this call, false otherwise. - template - bool prepareCurvilinearState( - [[maybe_unused]] propagator_state_t& prop_state, - [[maybe_unused]] const navigator_t& navigator) const { - return true; - } + bool prepareCurvilinearState(State& state) const; /// Create and return a curvilinear state at the current position /// @@ -367,15 +363,13 @@ class SympyStepper { /// Perform a Runge-Kutta track parameter propagation step /// - /// @param [in,out] state the propagation state - /// @param [in] navigator the navigator of the propagation + /// @param [in,out] state State of the stepper /// @note The state contains the desired step size. It can be negative during /// backwards track propagation, and since we're using an adaptive /// algorithm, it can be modified by the stepper class during /// propagation. - template - Result step(propagator_state_t& state, - const navigator_t& navigator) const; + Result step(State& state, Direction propDir, + const IVolumeMaterial* material) const; /// Method that reset the Jacobian to the Identity for when no bound state are /// available @@ -388,15 +382,12 @@ class SympyStepper { std::shared_ptr m_bField; private: - Result stepImpl(State& state, Direction stepDirection, - double stepTolerance, double stepSizeCutOff, + Result stepImpl(State& state, Direction propDir, double stepTolerance, + double stepSizeCutOff, std::size_t maxRungeKuttaStepTrials) const; }; -template -struct SupportsBoundParameters - : public std::true_type {}; +template <> +struct SupportsBoundParameters : public std::true_type {}; } // namespace Acts - -#include "Acts/Propagator/SympyStepper.ipp" diff --git a/Core/include/Acts/Propagator/SympyStepper.ipp b/Core/include/Acts/Propagator/SympyStepper.ipp deleted file mode 100644 index 30e91427026..00000000000 --- a/Core/include/Acts/Propagator/SympyStepper.ipp +++ /dev/null @@ -1,22 +0,0 @@ -// This file is part of the ACTS project. -// -// Copyright (C) 2016 CERN for the benefit of the ACTS project -// -// This Source Code Form is subject to the terms of the Mozilla Public -// License, v. 2.0. If a copy of the MPL was not distributed with this -// file, You can obtain one at https://mozilla.org/MPL/2.0/. - -#include "Acts/Propagator/EigenStepperError.hpp" - -namespace Acts { - -template -Result SympyStepper::step(propagator_state_t& state, - const navigator_t& /*navigator*/) const { - return stepImpl(state.stepping, state.options.direction, - state.options.stepping.stepTolerance, - state.options.stepping.stepSizeCutOff, - state.options.stepping.maxRungeKuttaStepTrials); -} - -} // namespace Acts diff --git a/Core/include/Acts/Propagator/VoidNavigator.hpp b/Core/include/Acts/Propagator/VoidNavigator.hpp index d91b75450ec..1f1c2dfc56e 100644 --- a/Core/include/Acts/Propagator/VoidNavigator.hpp +++ b/Core/include/Acts/Propagator/VoidNavigator.hpp @@ -16,6 +16,8 @@ namespace Acts { +class TrackingVolume; +class IVolumeMaterial; class Surface; /// @brief A navigator that does nothing @@ -56,6 +58,14 @@ class VoidNavigator { return nullptr; } + const TrackingVolume* currentVolume(const State& /*state*/) const { + return nullptr; + } + + const IVolumeMaterial* currentVolumeMaterial(const State& /*state*/) const { + return nullptr; + } + const Surface* startSurface(const State& /*state*/) const { return nullptr; } const Surface* targetSurface(const State& /*state*/) const { return nullptr; } diff --git a/Core/include/Acts/Propagator/detail/LoopStepperUtils.hpp b/Core/include/Acts/Propagator/detail/LoopStepperUtils.hpp index 26ddb590a8f..02b9d2e5b41 100644 --- a/Core/include/Acts/Propagator/detail/LoopStepperUtils.hpp +++ b/Core/include/Acts/Propagator/detail/LoopStepperUtils.hpp @@ -8,6 +8,11 @@ #pragma once +#include +#include +#include +#include + namespace Acts::detail { /// A helper type for providinig a propagation state which can be used with @@ -41,6 +46,7 @@ struct LoopComponentProxyBase { // These are the const accessors, which are shared between the mutable // ComponentProxy and the ConstComponentProxy + const auto& state() const { return cmp.state; } auto status() const { return cmp.status; } auto weight() const { return cmp.weight; } auto pathAccumulated() const { return cmp.state.pathAccumulated; } @@ -92,6 +98,7 @@ struct LoopComponentProxy using Base::pathAccumulated; using Base::singleState; using Base::singleStepper; + using Base::state; using Base::status; using Base::weight; @@ -103,6 +110,7 @@ struct LoopComponentProxy // These are the mutable accessors, the const ones are inherited from the // ComponentProxyBase + auto& state() { return cmp.state; } auto& status() { return cmp.status; } auto& weight() { return cmp.weight; } auto& pathAccumulated() { return cmp.state.pathAccumulated; } diff --git a/Core/include/Acts/TrackFitting/KalmanFitter.hpp b/Core/include/Acts/TrackFitting/KalmanFitter.hpp index 1bc0858b3df..4832494e64f 100644 --- a/Core/include/Acts/TrackFitting/KalmanFitter.hpp +++ b/Core/include/Acts/TrackFitting/KalmanFitter.hpp @@ -657,7 +657,7 @@ class KalmanFitter { // detected or if the surface has material (no holes before the first // measurement) auto trackStateProxyRes = detail::kalmanHandleNoMeasurement( - state, stepper, *surface, *result.fittedStates, + state.stepping, stepper, *surface, *result.fittedStates, result.lastTrackIndex, true, logger(), precedingMeasurementExists, freeToBoundCorrection); diff --git a/Core/include/Acts/TrackFitting/detail/GsfActor.hpp b/Core/include/Acts/TrackFitting/detail/GsfActor.hpp index 0b3a8053870..9791cb55ca6 100644 --- a/Core/include/Acts/TrackFitting/detail/GsfActor.hpp +++ b/Core/include/Acts/TrackFitting/detail/GsfActor.hpp @@ -235,9 +235,8 @@ struct GsfActor { } for (auto cmp : stepper.componentIterable(state.stepping)) { - auto singleState = cmp.singleState(state); - cmp.singleStepper(stepper).transportCovarianceToBound( - singleState.stepping, surface); + cmp.singleStepper(stepper).transportCovarianceToBound(cmp.state(), + surface); } if (haveMaterial) { @@ -640,9 +639,8 @@ struct GsfActor { // components should behave the same bool isHole = true; - auto cmps = stepper.componentIterable(state.stepping); - for (auto cmp : cmps) { - auto singleState = cmp.singleState(state); + for (auto cmp : stepper.componentIterable(state.stepping)) { + auto& singleState = cmp.state(); const auto& singleStepper = cmp.singleStepper(stepper); // There is some redundant checking inside this function, but do this for diff --git a/Core/include/Acts/TrackFitting/detail/KalmanUpdateHelpers.hpp b/Core/include/Acts/TrackFitting/detail/KalmanUpdateHelpers.hpp index f4230e0e18d..1098c342672 100644 --- a/Core/include/Acts/TrackFitting/detail/KalmanUpdateHelpers.hpp +++ b/Core/include/Acts/TrackFitting/detail/KalmanUpdateHelpers.hpp @@ -10,7 +10,6 @@ #include "Acts/EventData/MultiTrajectory.hpp" #include "Acts/EventData/SourceLink.hpp" -#include "Acts/EventData/TrackParameters.hpp" #include "Acts/EventData/detail/CorrectedTransformationFreeToBound.hpp" #include "Acts/Surfaces/Surface.hpp" #include "Acts/Utilities/CalibrationContext.hpp" @@ -121,7 +120,6 @@ auto kalmanHandleMeasurement( /// If there are no source links on surface, add either a hole or passive /// material TrackState entry multi trajectory. No storage allocation for /// uncalibrated/calibrated measurement and filtered parameter -/// @tparam propagator_state_t The propagator state type /// @tparam stepper_t The stepper type /// @param state The propagator state /// @param stepper The stepper @@ -131,11 +129,12 @@ auto kalmanHandleMeasurement( /// @param doCovTransport Whether to perform a covariance transport when /// computing the bound state or not /// @param freeToBoundCorrection Correction for non-linearity effect during transform from free to bound (only corrected when performing CovTransport) -template +template auto kalmanHandleNoMeasurement( - propagator_state_t &state, const stepper_t &stepper, const Surface &surface, - traj_t &fittedStates, const std::size_t lastTrackIndex, bool doCovTransport, - const Logger &logger, const bool precedingMeasurementExists, + typename stepper_t::State &state, const stepper_t &stepper, + const Surface &surface, traj_t &fittedStates, + const std::size_t lastTrackIndex, bool doCovTransport, const Logger &logger, + const bool precedingMeasurementExists, const FreeToBoundCorrection &freeToBoundCorrection = FreeToBoundCorrection( false)) -> Result { // Add a TrackState entry multi trajectory. This allocates storage for @@ -151,7 +150,7 @@ auto kalmanHandleNoMeasurement( { trackStateProxy.setReferenceSurface(surface.getSharedPtr()); // Bind the transported state to the current surface - auto res = stepper.boundState(state.stepping, surface, doCovTransport, + auto res = stepper.boundState(state, surface, doCovTransport, freeToBoundCorrection); if (!res.ok()) { return res.error(); @@ -160,7 +159,7 @@ auto kalmanHandleNoMeasurement( // Fill the track state trackStateProxy.predicted() = boundParams.parameters(); - trackStateProxy.predictedCovariance() = state.stepping.cov; + trackStateProxy.predictedCovariance() = state.cov; trackStateProxy.jacobian() = jacobian; trackStateProxy.pathLength() = pathLength; diff --git a/Core/src/Propagator/SympyStepper.cpp b/Core/src/Propagator/SympyStepper.cpp index 932c8782c7b..4fb155a39fb 100644 --- a/Core/src/Propagator/SympyStepper.cpp +++ b/Core/src/Propagator/SympyStepper.cpp @@ -8,6 +8,7 @@ #include "Acts/Propagator/SympyStepper.hpp" +#include "Acts/Propagator/EigenStepperError.hpp" #include "Acts/Propagator/detail/SympyCovarianceEngine.hpp" #include "Acts/Propagator/detail/SympyJacobianEngine.hpp" @@ -37,7 +38,6 @@ SympyStepper::State SympyStepper::makeState( // Init the jacobian matrix if needed if (par.covariance()) { - // Get the reference surface for navigation const auto& surface = par.referenceSurface(); // set the covariance transport flag to true and copy state.covTransport = true; @@ -84,6 +84,12 @@ SympyStepper::boundState( state.pathAccumulated, freeToBoundCorrection); } +bool SympyStepper::prepareCurvilinearState(State& state) const { + // TODO implement like in EigenStepper + (void)state; + return true; +} + std::tuple SympyStepper::curvilinearState(State& state, bool transportCov) const { return detail::sympy::curvilinearState( @@ -127,8 +133,16 @@ void SympyStepper::transportCovarianceToBound( freeToBoundCorrection); } +Result SympyStepper::step(State& state, Direction propDir, + const IVolumeMaterial* material) const { + (void)material; + return stepImpl(state, propDir, state.options.stepTolerance, + state.options.maxStepSize, + state.options.maxRungeKuttaStepTrials); +} + Result SympyStepper::stepImpl( - State& state, Direction stepDirection, double stepTolerance, + State& state, Direction propDir, double stepTolerance, double stepSizeCutOff, std::size_t maxRungeKuttaStepTrials) const { auto pos = position(state); auto dir = direction(state); @@ -160,7 +174,7 @@ Result SympyStepper::stepImpl( return std::clamp(x, lower, upper); }; - double h = state.stepSize.value() * stepDirection; + double h = state.stepSize.value() * propDir; double initialH = h; std::size_t nStepTrials = 0; double errorEstimate = 0.; @@ -212,7 +226,7 @@ Result SympyStepper::stepImpl( state.nStepTrials += nStepTrials; ++state.statistics.nSuccessfulSteps; - if (stepDirection != Direction::fromScalarZeroAsPositive(initialH)) { + if (propDir != Direction::fromScalarZeroAsPositive(initialH)) { ++state.statistics.nReverseSteps; } state.statistics.pathLength += h; diff --git a/Tests/UnitTests/Core/Propagator/AtlasStepperTests.cpp b/Tests/UnitTests/Core/Propagator/AtlasStepperTests.cpp index a5dc3ecc331..93af62cfa95 100644 --- a/Tests/UnitTests/Core/Propagator/AtlasStepperTests.cpp +++ b/Tests/UnitTests/Core/Propagator/AtlasStepperTests.cpp @@ -55,26 +55,6 @@ using Covariance = BoundSquareMatrix; using Jacobian = BoundMatrix; using Stepper = Acts::AtlasStepper; -/// Simplified propagator state. -struct MockPropagatorState { - MockPropagatorState(Stepper::State stepperState) - : stepping(std::move(stepperState)) {} - - /// Stepper state. - Stepper::State stepping; - /// Propagator options with only the relevant components. - struct { - Direction direction = Direction::Backward(); - struct { - double stepTolerance = 10_um; - } stepping; - } options; -}; - -struct MockNavigator {}; - -static constexpr MockNavigator navigator; - // epsilon for floating point comparisons static constexpr auto eps = 1024 * std::numeric_limits::epsilon(); @@ -310,31 +290,31 @@ BOOST_AUTO_TEST_CASE(Step) { Stepper::Options options(geoCtx, magCtx); options.maxStepSize = stepSize; - MockPropagatorState state(stepper.makeState(options, cp)); - state.stepping.covTransport = false; + auto state = stepper.makeState(options, cp); + state.covTransport = false; // ensure step does not result in an error - auto res = stepper.step(state, navigator); + auto res = stepper.step(state, Direction::Forward(), nullptr); BOOST_CHECK(res.ok()); // extract the actual step size auto h = res.value(); - BOOST_CHECK_EQUAL(state.stepping.stepSize.value(), stepSize); - BOOST_CHECK_EQUAL(state.stepping.stepSize.value(), h * navDir); + BOOST_CHECK_EQUAL(state.stepSize.value(), stepSize); + BOOST_CHECK_EQUAL(state.stepSize.value(), h * navDir); // check that the position has moved - auto deltaPos = (stepper.position(state.stepping) - pos).eval(); + auto deltaPos = (stepper.position(state) - pos).eval(); BOOST_CHECK_LT(0, deltaPos.norm()); // check that time has changed - auto deltaTime = stepper.time(state.stepping) - time; + auto deltaTime = stepper.time(state) - time; BOOST_CHECK_LT(0, std::abs(deltaTime)); // check that the direction was rotated - auto projDir = stepper.direction(state.stepping).dot(unitDir); + auto projDir = stepper.direction(state).dot(unitDir); BOOST_CHECK_LT(projDir, 1); // momentum and charge should not change - CHECK_CLOSE_ABS(stepper.absoluteMomentum(state.stepping), absMom, eps); - BOOST_CHECK_EQUAL(stepper.charge(state.stepping), charge); + CHECK_CLOSE_ABS(stepper.absoluteMomentum(state), absMom, eps); + BOOST_CHECK_EQUAL(stepper.charge(state), charge); } // test step method with covariance transport @@ -347,34 +327,34 @@ BOOST_AUTO_TEST_CASE(StepWithCovariance) { Stepper::Options options(geoCtx, magCtx); options.maxStepSize = stepSize; - MockPropagatorState state(stepper.makeState(options, cp)); - state.stepping.covTransport = true; + auto state = stepper.makeState(options, cp); + state.covTransport = true; // ensure step does not result in an error - auto res = stepper.step(state, navigator); + auto res = stepper.step(state, Direction::Forward(), nullptr); BOOST_CHECK(res.ok()); // extract the actual step size auto h = res.value(); - BOOST_CHECK_EQUAL(state.stepping.stepSize.value(), stepSize); - BOOST_CHECK_EQUAL(state.stepping.stepSize.value(), h * navDir); + BOOST_CHECK_EQUAL(state.stepSize.value(), stepSize); + BOOST_CHECK_EQUAL(state.stepSize.value(), h * navDir); // check that the position has moved - auto deltaPos = (stepper.position(state.stepping) - pos).eval(); + auto deltaPos = (stepper.position(state) - pos).eval(); BOOST_CHECK_LT(0, deltaPos.norm()); // check that time has changed - auto deltaTime = stepper.time(state.stepping) - time; + auto deltaTime = stepper.time(state) - time; BOOST_CHECK_LT(0, std::abs(deltaTime)); // check that the direction was rotated - auto projDir = stepper.direction(state.stepping).dot(unitDir); + auto projDir = stepper.direction(state).dot(unitDir); BOOST_CHECK_LT(projDir, 1); // momentum and charge should not change - CHECK_CLOSE_ABS(stepper.absoluteMomentum(state.stepping), absMom, eps); - BOOST_CHECK_EQUAL(stepper.charge(state.stepping), charge); + CHECK_CLOSE_ABS(stepper.absoluteMomentum(state), absMom, eps); + BOOST_CHECK_EQUAL(stepper.charge(state), charge); - stepper.transportCovarianceToCurvilinear(state.stepping); - BOOST_CHECK_NE(state.stepping.cov, cov); + stepper.transportCovarianceToCurvilinear(state); + BOOST_CHECK_NE(state.cov, cov); } // test state reset method @@ -387,11 +367,11 @@ BOOST_AUTO_TEST_CASE(Reset) { Stepper::Options options(geoCtx, magCtx); options.maxStepSize = stepSize; - MockPropagatorState state(stepper.makeState(options, cp)); - state.stepping.covTransport = true; + auto state = stepper.makeState(options, cp); + state.covTransport = true; // ensure step does not result in an error - stepper.step(state, navigator); + stepper.step(state, Direction::Forward(), nullptr); // Construct the parameters Vector3 newPos(1.5, -2.5, 3.5); @@ -443,7 +423,7 @@ BOOST_AUTO_TEST_CASE(Reset) { }; // Reset all possible parameters - Stepper::State stateCopy = copyState(*magneticField, state.stepping); + Stepper::State stateCopy = copyState(*magneticField, state); BOOST_CHECK(cp.covariance().has_value()); stepper.resetState(stateCopy, cp.parameters(), *cp.covariance(), cp.referenceSurface(), stepSize); @@ -456,15 +436,14 @@ BOOST_AUTO_TEST_CASE(Reset) { freeParams.template segment<3>(eFreeDir0).normalized()); BOOST_CHECK_EQUAL(stepper.absoluteMomentum(stateCopy), std::abs(1. / freeParams[eFreeQOverP])); - BOOST_CHECK_EQUAL(stepper.charge(stateCopy), stepper.charge(state.stepping)); + BOOST_CHECK_EQUAL(stepper.charge(stateCopy), stepper.charge(state)); BOOST_CHECK_EQUAL(stepper.time(stateCopy), freeParams[eFreeTime]); BOOST_CHECK_EQUAL(stateCopy.pathAccumulated, 0.); BOOST_CHECK_EQUAL(stateCopy.stepSize.value(), navDir * stepSize); - BOOST_CHECK_EQUAL(stateCopy.previousStepSize, - state.stepping.previousStepSize); + BOOST_CHECK_EQUAL(stateCopy.previousStepSize, state.previousStepSize); // Reset all possible parameters except the step size - stateCopy = copyState(*magneticField, state.stepping); + stateCopy = copyState(*magneticField, state); stepper.resetState(stateCopy, cp.parameters(), *cp.covariance(), cp.referenceSurface()); // Test all components @@ -476,16 +455,15 @@ BOOST_AUTO_TEST_CASE(Reset) { freeParams.template segment<3>(eFreeDir0).normalized()); BOOST_CHECK_EQUAL(stepper.absoluteMomentum(stateCopy), std::abs(1. / freeParams[eFreeQOverP])); - BOOST_CHECK_EQUAL(stepper.charge(stateCopy), stepper.charge(state.stepping)); + BOOST_CHECK_EQUAL(stepper.charge(stateCopy), stepper.charge(state)); BOOST_CHECK_EQUAL(stepper.time(stateCopy), freeParams[eFreeTime]); BOOST_CHECK_EQUAL(stateCopy.pathAccumulated, 0.); BOOST_CHECK_EQUAL(stateCopy.stepSize.value(), std::numeric_limits::max()); - BOOST_CHECK_EQUAL(stateCopy.previousStepSize, - state.stepping.previousStepSize); + BOOST_CHECK_EQUAL(stateCopy.previousStepSize, state.previousStepSize); // Reset the least amount of parameters - stateCopy = copyState(*magneticField, state.stepping); + stateCopy = copyState(*magneticField, state); stepper.resetState(stateCopy, cp.parameters(), *cp.covariance(), cp.referenceSurface()); // Test all components @@ -497,13 +475,12 @@ BOOST_AUTO_TEST_CASE(Reset) { freeParams.template segment<3>(eFreeDir0).normalized()); BOOST_CHECK_EQUAL(stepper.absoluteMomentum(stateCopy), std::abs(1. / freeParams[eFreeQOverP])); - BOOST_CHECK_EQUAL(stepper.charge(stateCopy), stepper.charge(state.stepping)); + BOOST_CHECK_EQUAL(stepper.charge(stateCopy), stepper.charge(state)); BOOST_CHECK_EQUAL(stepper.time(stateCopy), freeParams[eFreeTime]); BOOST_CHECK_EQUAL(stateCopy.pathAccumulated, 0.); BOOST_CHECK_EQUAL(stateCopy.stepSize.value(), std::numeric_limits::max()); - BOOST_CHECK_EQUAL(stateCopy.previousStepSize, - state.stepping.previousStepSize); + BOOST_CHECK_EQUAL(stateCopy.previousStepSize, state.previousStepSize); // Reset using different surface shapes // 1) Disc surface @@ -521,13 +498,13 @@ BOOST_AUTO_TEST_CASE(Reset) { .value(); // Reset the state and test - Stepper::State stateDisc = copyState(*magneticField, state.stepping); + Stepper::State stateDisc = copyState(*magneticField, state); BOOST_CHECK(boundDisc.covariance().has_value()); stepper.resetState(stateDisc, boundDisc.parameters(), *boundDisc.covariance(), boundDisc.referenceSurface()); CHECK_NE_COLLECTIONS(stateDisc.pVector, stateCopy.pVector); - CHECK_NE_COLLECTIONS(stateDisc.pVector, state.stepping.pVector); + CHECK_NE_COLLECTIONS(stateDisc.pVector, state.pVector); // 2) Perigee surface // Setting some parameters @@ -544,13 +521,13 @@ BOOST_AUTO_TEST_CASE(Reset) { .value(); // Reset the state and test - Stepper::State statePerigee = copyState(*magneticField, state.stepping); + Stepper::State statePerigee = copyState(*magneticField, state); BOOST_CHECK(boundPerigee.covariance().has_value()); stepper.resetState(statePerigee, boundPerigee.parameters(), *boundPerigee.covariance(), boundPerigee.referenceSurface()); CHECK_NE_COLLECTIONS(statePerigee.pVector, stateCopy.pVector); - CHECK_NE_COLLECTIONS(statePerigee.pVector, state.stepping.pVector); + CHECK_NE_COLLECTIONS(statePerigee.pVector, state.pVector); CHECK_NE_COLLECTIONS(statePerigee.pVector, stateDisc.pVector); // 3) Straw surface @@ -562,12 +539,12 @@ BOOST_AUTO_TEST_CASE(Reset) { .value(); // Reset the state and test - Stepper::State stateStraw = copyState(*magneticField, state.stepping); + Stepper::State stateStraw = copyState(*magneticField, state); BOOST_CHECK(boundStraw.covariance().has_value()); stepper.resetState(stateStraw, boundStraw.parameters(), *boundStraw.covariance(), boundStraw.referenceSurface()); CHECK_NE_COLLECTIONS(stateStraw.pVector, stateCopy.pVector); - CHECK_NE_COLLECTIONS(stateStraw.pVector, state.stepping.pVector); + CHECK_NE_COLLECTIONS(stateStraw.pVector, state.pVector); CHECK_NE_COLLECTIONS(stateStraw.pVector, stateDisc.pVector); BOOST_CHECK_EQUAL_COLLECTIONS( stateStraw.pVector.begin(), stateStraw.pVector.end(), diff --git a/Tests/UnitTests/Core/Propagator/EigenStepperTests.cpp b/Tests/UnitTests/Core/Propagator/EigenStepperTests.cpp index 33ff3a68f29..909dc36b4da 100644 --- a/Tests/UnitTests/Core/Propagator/EigenStepperTests.cpp +++ b/Tests/UnitTests/Core/Propagator/EigenStepperTests.cpp @@ -50,7 +50,6 @@ #include "Acts/Utilities/Result.hpp" #include "Acts/Utilities/UnitVectors.hpp" -#include #include #include #include @@ -74,31 +73,6 @@ static constexpr auto eps = 3 * std::numeric_limits::epsilon(); GeometryContext tgContext = GeometryContext(); MagneticFieldContext mfContext = MagneticFieldContext(); -/// @brief Simplified propagator state -template -struct PropState { - /// @brief Constructor - PropState(Direction direction, stepper_state_t sState) - : stepping(std::move(sState)) { - options.direction = direction; - } - /// State of the eigen stepper - stepper_state_t stepping; - /// Propagator options which only carry the relevant components - struct { - Direction direction = Direction::Forward(); - struct { - double stepTolerance = 1e-4; - double stepSizeCutOff = 0.; - unsigned int maxRungeKuttaStepTrials = 10000; - } stepping; - } options; -}; - -struct MockNavigator {}; - -static constexpr MockNavigator mockNavigator; - /// @brief Aborter for the case that a particle leaves the detector or reaches /// a custom made threshold. /// @@ -303,27 +277,26 @@ BOOST_AUTO_TEST_CASE(eigen_stepper_test) { // Perform a step without and with covariance transport esState.cov = cov; - PropState ps(navDir, std::move(esState)); - - ps.stepping.covTransport = false; - es.step(ps, mockNavigator).value(); - CHECK_CLOSE_COVARIANCE(ps.stepping.cov, cov, eps); - BOOST_CHECK_NE(es.position(ps.stepping).norm(), newPos.norm()); - BOOST_CHECK_NE(es.direction(ps.stepping), newMom.normalized()); - BOOST_CHECK_EQUAL(es.charge(ps.stepping), charge); - BOOST_CHECK_LT(es.time(ps.stepping), newTime); - BOOST_CHECK_EQUAL(ps.stepping.derivative, FreeVector::Zero()); - BOOST_CHECK_EQUAL(ps.stepping.jacTransport, FreeMatrix::Identity()); - - ps.stepping.covTransport = true; - es.step(ps, mockNavigator).value(); - CHECK_CLOSE_COVARIANCE(ps.stepping.cov, cov, eps); - BOOST_CHECK_NE(es.position(ps.stepping).norm(), newPos.norm()); - BOOST_CHECK_NE(es.direction(ps.stepping), newMom.normalized()); - BOOST_CHECK_EQUAL(es.charge(ps.stepping), charge); - BOOST_CHECK_LT(es.time(ps.stepping), newTime); - BOOST_CHECK_NE(ps.stepping.derivative, FreeVector::Zero()); - BOOST_CHECK_NE(ps.stepping.jacTransport, FreeMatrix::Identity()); + + esState.covTransport = false; + es.step(esState, Direction::Forward(), nullptr).value(); + CHECK_CLOSE_COVARIANCE(esState.cov, cov, eps); + BOOST_CHECK_NE(es.position(esState).norm(), newPos.norm()); + BOOST_CHECK_NE(es.direction(esState), newMom.normalized()); + BOOST_CHECK_EQUAL(es.charge(esState), charge); + BOOST_CHECK_LT(es.time(esState), newTime); + BOOST_CHECK_EQUAL(esState.derivative, FreeVector::Zero()); + BOOST_CHECK_EQUAL(esState.jacTransport, FreeMatrix::Identity()); + + esState.covTransport = true; + es.step(esState, Direction::Forward(), nullptr).value(); + CHECK_CLOSE_COVARIANCE(esState.cov, cov, eps); + BOOST_CHECK_NE(es.position(esState).norm(), newPos.norm()); + BOOST_CHECK_NE(es.direction(esState), newMom.normalized()); + BOOST_CHECK_EQUAL(es.charge(esState), charge); + BOOST_CHECK_LT(es.time(esState), newTime); + BOOST_CHECK_NE(esState.derivative, FreeVector::Zero()); + BOOST_CHECK_NE(esState.jacTransport, FreeMatrix::Identity()); /// Test the state reset // Construct the parameters @@ -366,13 +339,13 @@ BOOST_AUTO_TEST_CASE(eigen_stepper_test) { }; // Reset all possible parameters - EigenStepper<>::State esStateCopy = copyState(*bField, ps.stepping); + EigenStepper<>::State esStateCopy = copyState(*bField, esState); BOOST_CHECK(cp2.covariance().has_value()); es.resetState(esStateCopy, cp2.parameters(), *cp2.covariance(), cp2.referenceSurface(), stepSize2); // Test all components BOOST_CHECK_NE(esStateCopy.jacToGlobal, BoundToFreeMatrix::Zero()); - BOOST_CHECK_NE(esStateCopy.jacToGlobal, ps.stepping.jacToGlobal); + BOOST_CHECK_NE(esStateCopy.jacToGlobal, esState.jacToGlobal); BOOST_CHECK_EQUAL(esStateCopy.jacTransport, FreeMatrix::Identity()); BOOST_CHECK_EQUAL(esStateCopy.derivative, FreeVector::Zero()); BOOST_CHECK(esStateCopy.covTransport); @@ -383,19 +356,19 @@ BOOST_AUTO_TEST_CASE(eigen_stepper_test) { freeParams.template segment<3>(eFreeDir0).normalized()); BOOST_CHECK_EQUAL(es.absoluteMomentum(esStateCopy), std::abs(1. / freeParams[eFreeQOverP])); - BOOST_CHECK_EQUAL(es.charge(esStateCopy), -es.charge(ps.stepping)); + BOOST_CHECK_EQUAL(es.charge(esStateCopy), -es.charge(esState)); BOOST_CHECK_EQUAL(es.time(esStateCopy), freeParams[eFreeTime]); BOOST_CHECK_EQUAL(esStateCopy.pathAccumulated, 0.); BOOST_CHECK_EQUAL(esStateCopy.stepSize.value(), navDir * stepSize2); - BOOST_CHECK_EQUAL(esStateCopy.previousStepSize, ps.stepping.previousStepSize); + BOOST_CHECK_EQUAL(esStateCopy.previousStepSize, esState.previousStepSize); // Reset all possible parameters except the step size - esStateCopy = copyState(*bField, ps.stepping); + esStateCopy = copyState(*bField, esState); es.resetState(esStateCopy, cp2.parameters(), *cp2.covariance(), cp2.referenceSurface()); // Test all components BOOST_CHECK_NE(esStateCopy.jacToGlobal, BoundToFreeMatrix::Zero()); - BOOST_CHECK_NE(esStateCopy.jacToGlobal, ps.stepping.jacToGlobal); + BOOST_CHECK_NE(esStateCopy.jacToGlobal, esState.jacToGlobal); BOOST_CHECK_EQUAL(esStateCopy.jacTransport, FreeMatrix::Identity()); BOOST_CHECK_EQUAL(esStateCopy.derivative, FreeVector::Zero()); BOOST_CHECK(esStateCopy.covTransport); @@ -406,20 +379,20 @@ BOOST_AUTO_TEST_CASE(eigen_stepper_test) { freeParams.template segment<3>(eFreeDir0).normalized()); BOOST_CHECK_EQUAL(es.absoluteMomentum(esStateCopy), std::abs(1. / freeParams[eFreeQOverP])); - BOOST_CHECK_EQUAL(es.charge(esStateCopy), -es.charge(ps.stepping)); + BOOST_CHECK_EQUAL(es.charge(esStateCopy), -es.charge(esState)); BOOST_CHECK_EQUAL(es.time(esStateCopy), freeParams[eFreeTime]); BOOST_CHECK_EQUAL(esStateCopy.pathAccumulated, 0.); BOOST_CHECK_EQUAL(esStateCopy.stepSize.value(), std::numeric_limits::max()); - BOOST_CHECK_EQUAL(esStateCopy.previousStepSize, ps.stepping.previousStepSize); + BOOST_CHECK_EQUAL(esStateCopy.previousStepSize, esState.previousStepSize); // Reset the least amount of parameters - esStateCopy = copyState(*bField, ps.stepping); + esStateCopy = copyState(*bField, esState); es.resetState(esStateCopy, cp2.parameters(), *cp2.covariance(), cp2.referenceSurface()); // Test all components BOOST_CHECK_NE(esStateCopy.jacToGlobal, BoundToFreeMatrix::Zero()); - BOOST_CHECK_NE(esStateCopy.jacToGlobal, ps.stepping.jacToGlobal); + BOOST_CHECK_NE(esStateCopy.jacToGlobal, esState.jacToGlobal); BOOST_CHECK_EQUAL(esStateCopy.jacTransport, FreeMatrix::Identity()); BOOST_CHECK_EQUAL(esStateCopy.derivative, FreeVector::Zero()); BOOST_CHECK(esStateCopy.covTransport); @@ -430,12 +403,12 @@ BOOST_AUTO_TEST_CASE(eigen_stepper_test) { freeParams.template segment<3>(eFreeDir0).normalized()); BOOST_CHECK_EQUAL(es.absoluteMomentum(esStateCopy), std::abs(1. / freeParams[eFreeQOverP])); - BOOST_CHECK_EQUAL(es.charge(esStateCopy), -es.charge(ps.stepping)); + BOOST_CHECK_EQUAL(es.charge(esStateCopy), -es.charge(esState)); BOOST_CHECK_EQUAL(es.time(esStateCopy), freeParams[eFreeTime]); BOOST_CHECK_EQUAL(esStateCopy.pathAccumulated, 0.); BOOST_CHECK_EQUAL(esStateCopy.stepSize.value(), std::numeric_limits::max()); - BOOST_CHECK_EQUAL(esStateCopy.previousStepSize, ps.stepping.previousStepSize); + BOOST_CHECK_EQUAL(esStateCopy.previousStepSize, esState.previousStepSize); /// Repeat with surface related methods auto plane = CurvilinearSurface(pos, dir.normalized()).planeSurface(); @@ -508,9 +481,9 @@ BOOST_AUTO_TEST_CASE(eigen_stepper_test) { CHECK_CLOSE_COVARIANCE(esState.cov, Covariance(2. * cov), eps); // Test a case where no step size adjustment is required - ps.options.stepping.stepTolerance = 2. * 4.4258e+09; + esState.options.stepTolerance = 2. * 4.4258e+09; double h0 = esState.stepSize.value(); - es.step(ps, mockNavigator); + es.step(esState, navDir, nullptr); CHECK_CLOSE_ABS(h0, esState.stepSize.value(), eps); // Produce some errors @@ -518,18 +491,18 @@ BOOST_AUTO_TEST_CASE(eigen_stepper_test) { EigenStepper<> nes(nBfield); EigenStepper<>::Options nesOptions(tgContext, mfContext); EigenStepper<>::State nesState = nes.makeState(nesOptions, cp); - PropState nps(navDir, copyState(*nBfield, nesState)); + auto nEsState = copyState(*nBfield, nesState); // Test that we can reach the minimum step size - nps.options.stepping.stepTolerance = 1e-21; - nps.options.stepping.stepSizeCutOff = 1e20; - auto res = nes.step(nps, mockNavigator); + nEsState.options.stepTolerance = 1e-21; + nEsState.options.stepSizeCutOff = 1e20; + auto res = nes.step(nEsState, navDir, nullptr); BOOST_CHECK(!res.ok()); BOOST_CHECK_EQUAL(res.error(), EigenStepperError::StepSizeStalled); // Test that the number of trials exceeds - nps.options.stepping.stepSizeCutOff = 0.; - nps.options.stepping.maxRungeKuttaStepTrials = 0.; - res = nes.step(nps, mockNavigator); + nEsState.options.stepSizeCutOff = 0.; + nEsState.options.maxRungeKuttaStepTrials = 0.; + res = nes.step(nEsState, navDir, nullptr); BOOST_CHECK(!res.ok()); BOOST_CHECK_EQUAL(res.error(), EigenStepperError::StepSizeAdjustmentFailed); } diff --git a/Tests/UnitTests/Core/Propagator/MultiStepperTests.cpp b/Tests/UnitTests/Core/Propagator/MultiStepperTests.cpp index 992e39c8d2e..14b064670a6 100644 --- a/Tests/UnitTests/Core/Propagator/MultiStepperTests.cpp +++ b/Tests/UnitTests/Core/Propagator/MultiStepperTests.cpp @@ -30,7 +30,6 @@ #include "Acts/Surfaces/PlaneSurface.hpp" #include "Acts/Surfaces/Surface.hpp" #include "Acts/Utilities/Intersection.hpp" -#include "Acts/Utilities/Logger.hpp" #include #include @@ -75,40 +74,6 @@ const auto defaultNullBField = std::make_shared(); const auto particleHypothesis = ParticleHypothesis::pion(); -struct Options { - Direction direction = defaultNDir; - - const Acts::Logger &logger = Acts::getDummyLogger(); - - struct { - double stepTolerance = 1e-4; - double stepSizeCutOff = 0.0; - std::size_t maxRungeKuttaStepTrials = 10; - } stepping; -}; - -struct MockNavigator {}; - -static constexpr MockNavigator mockNavigator; - -struct Navigation {}; - -template -struct DummyPropState { - stepper_state_t &stepping; - Options options; - Navigation navigation; - GeometryContext geoContext; - - DummyPropState(Direction direction, stepper_state_t &ss) - : stepping(ss), - options(Options{}), - navigation(Navigation{}), - geoContext(geoCtx) { - options.direction = direction; - } -}; - // Makes random bound parameters and covariance and a plane surface at {0,0,0} // with normal {1,0,0}. Optionally some external fixed bound parameters can be // supplied @@ -336,13 +301,12 @@ void test_multi_stepper_vs_eigen_stepper() { // Do some steps and check that the results match for (int i = 0; i < 10; ++i) { // Single stepper - auto single_prop_state = DummyPropState(defaultNDir, single_state); - auto single_result = single_stepper.step(single_prop_state, mockNavigator); + auto single_result = + single_stepper.step(single_state, defaultNDir, nullptr); single_stepper.transportCovarianceToCurvilinear(single_state); // Multi stepper; - auto multi_prop_state = DummyPropState(defaultNDir, multi_state); - auto multi_result = multi_stepper.step(multi_prop_state, mockNavigator); + auto multi_result = multi_stepper.step(multi_state, defaultNDir, nullptr); multi_stepper.transportCovarianceToCurvilinear(multi_state); // Check equality @@ -526,12 +490,10 @@ void test_multi_stepper_surface_status_update() { // Step forward now { - auto multi_prop_state = DummyPropState(Direction::Forward(), multi_state); - multi_stepper.step(multi_prop_state, mockNavigator); + multi_stepper.step(multi_state, Direction::Forward(), nullptr); // Single stepper - auto single_prop_state = DummyPropState(Direction::Forward(), single_state); - single_stepper.step(single_prop_state, mockNavigator); + single_stepper.step(single_state, Direction::Forward(), nullptr); } // Update surface status and check again @@ -627,16 +589,14 @@ void test_component_bound_state() { multi_state, *right_surface, 0, Direction::Forward(), BoundaryTolerance::Infinite(), s_onSurfaceTolerance, ConstrainedStep::Type::Navigator); - auto multi_prop_state = DummyPropState(Direction::Forward(), multi_state); - multi_stepper.step(multi_prop_state, mockNavigator); + multi_stepper.step(multi_state, Direction::Forward(), nullptr); // Single stepper single_stepper.updateSurfaceStatus( single_state, *right_surface, 0, Direction::Forward(), BoundaryTolerance::Infinite(), s_onSurfaceTolerance, ConstrainedStep::Type::Navigator); - auto single_prop_state = DummyPropState(Direction::Forward(), single_state); - single_stepper.step(single_prop_state, mockNavigator); + single_stepper.step(single_state, Direction::Forward(), nullptr); } // Check component-wise bound-state @@ -779,12 +739,10 @@ void test_single_component_interface_function() { MultiState multi_state = multi_stepper.makeState(options, multi_pars); - DummyPropState multi_prop_state(defaultNDir, multi_state); - // Check at least some properties at the moment auto check = [&](auto cmp) { auto sstepper = cmp.singleStepper(multi_stepper); - auto &sstepping = cmp.singleState(multi_prop_state).stepping; + auto &sstepping = cmp.state(); BOOST_CHECK_EQUAL(sstepper.position(sstepping), cmp.pars().template segment<3>(eFreePos0)); diff --git a/Tests/UnitTests/Core/Propagator/StraightLineStepperTests.cpp b/Tests/UnitTests/Core/Propagator/StraightLineStepperTests.cpp index 2a9a81e9d5c..8ace9416308 100644 --- a/Tests/UnitTests/Core/Propagator/StraightLineStepperTests.cpp +++ b/Tests/UnitTests/Core/Propagator/StraightLineStepperTests.cpp @@ -31,7 +31,6 @@ #include #include #include -#include using Acts::VectorHelpers::makeVector4; @@ -39,25 +38,6 @@ namespace Acts::Test { using Covariance = BoundSquareMatrix; -/// @brief Simplified propagator state -struct PropState { - /// @brief Constructor - PropState(Direction direction, StraightLineStepper::State sState) - : stepping(std::move(sState)) { - options.direction = direction; - } - /// State of the straight line stepper - StraightLineStepper::State stepping; - /// Propagator options which only carry the particle's mass - struct { - Direction direction = Direction::Forward(); - } options; -}; - -struct MockNavigator {}; - -static constexpr MockNavigator mockNavigator; - static constexpr auto eps = 2 * std::numeric_limits::epsilon(); /// These tests are aiming to test whether the state setup is working properly @@ -192,33 +172,32 @@ BOOST_AUTO_TEST_CASE(straight_line_stepper_test) { // Perform a step without and with covariance transport slsState.cov = cov; - PropState ps(navDir, slsState); - - ps.stepping.covTransport = false; - double h = sls.step(ps, mockNavigator).value(); - BOOST_CHECK_EQUAL(ps.stepping.stepSize.value(), stepSize); - BOOST_CHECK_EQUAL(ps.stepping.stepSize.value(), h * navDir); - CHECK_CLOSE_COVARIANCE(ps.stepping.cov, cov, 1e-6); - BOOST_CHECK_GT(sls.position(ps.stepping).norm(), newPos.norm()); - CHECK_CLOSE_ABS(sls.direction(ps.stepping), newMom.normalized(), 1e-6); - CHECK_CLOSE_ABS(sls.absoluteMomentum(ps.stepping), newMom.norm(), 1e-6); - CHECK_CLOSE_ABS(sls.charge(ps.stepping), charge, 1e-6); - BOOST_CHECK_LT(sls.time(ps.stepping), newTime); - BOOST_CHECK_EQUAL(ps.stepping.derivative, FreeVector::Zero()); - BOOST_CHECK_EQUAL(ps.stepping.jacTransport, FreeMatrix::Identity()); - - ps.stepping.covTransport = true; - double h2 = sls.step(ps, mockNavigator).value(); - BOOST_CHECK_EQUAL(ps.stepping.stepSize.value(), stepSize); + + slsState.covTransport = false; + double h = sls.step(slsState, Direction::Forward(), nullptr).value(); + BOOST_CHECK_EQUAL(slsState.stepSize.value(), stepSize); + BOOST_CHECK_EQUAL(slsState.stepSize.value(), h * navDir); + CHECK_CLOSE_COVARIANCE(slsState.cov, cov, 1e-6); + BOOST_CHECK_GT(sls.position(slsState).norm(), newPos.norm()); + CHECK_CLOSE_ABS(sls.direction(slsState), newMom.normalized(), 1e-6); + CHECK_CLOSE_ABS(sls.absoluteMomentum(slsState), newMom.norm(), 1e-6); + CHECK_CLOSE_ABS(sls.charge(slsState), charge, 1e-6); + BOOST_CHECK_LT(sls.time(slsState), newTime); + BOOST_CHECK_EQUAL(slsState.derivative, FreeVector::Zero()); + BOOST_CHECK_EQUAL(slsState.jacTransport, FreeMatrix::Identity()); + + slsState.covTransport = true; + double h2 = sls.step(slsState, Direction::Forward(), nullptr).value(); + BOOST_CHECK_EQUAL(slsState.stepSize.value(), stepSize); BOOST_CHECK_EQUAL(h2, h); - CHECK_CLOSE_COVARIANCE(ps.stepping.cov, cov, 1e-6); - BOOST_CHECK_GT(sls.position(ps.stepping).norm(), newPos.norm()); - CHECK_CLOSE_ABS(sls.direction(ps.stepping), newMom.normalized(), 1e-6); - CHECK_CLOSE_ABS(sls.absoluteMomentum(ps.stepping), newMom.norm(), 1e-6); - CHECK_CLOSE_ABS(sls.charge(ps.stepping), charge, 1e-6); - BOOST_CHECK_LT(sls.time(ps.stepping), newTime); - BOOST_CHECK_NE(ps.stepping.derivative, FreeVector::Zero()); - BOOST_CHECK_NE(ps.stepping.jacTransport, FreeMatrix::Identity()); + CHECK_CLOSE_COVARIANCE(slsState.cov, cov, 1e-6); + BOOST_CHECK_GT(sls.position(slsState).norm(), newPos.norm()); + CHECK_CLOSE_ABS(sls.direction(slsState), newMom.normalized(), 1e-6); + CHECK_CLOSE_ABS(sls.absoluteMomentum(slsState), newMom.norm(), 1e-6); + CHECK_CLOSE_ABS(sls.charge(slsState), charge, 1e-6); + BOOST_CHECK_LT(sls.time(slsState), newTime); + BOOST_CHECK_NE(slsState.derivative, FreeVector::Zero()); + BOOST_CHECK_NE(slsState.jacTransport, FreeMatrix::Identity()); /// Test the state reset // Construct the parameters @@ -238,12 +217,12 @@ BOOST_AUTO_TEST_CASE(straight_line_stepper_test) { double stepSize2 = -2. * stepSize; // Reset all possible parameters - StraightLineStepper::State slsStateCopy = ps.stepping; + StraightLineStepper::State slsStateCopy = slsState; sls.resetState(slsStateCopy, cp2.parameters(), *cp2.covariance(), cp2.referenceSurface(), stepSize2); // Test all components BOOST_CHECK_NE(slsStateCopy.jacToGlobal, BoundToFreeMatrix::Zero()); - BOOST_CHECK_NE(slsStateCopy.jacToGlobal, ps.stepping.jacToGlobal); + BOOST_CHECK_NE(slsStateCopy.jacToGlobal, slsState.jacToGlobal); BOOST_CHECK_EQUAL(slsStateCopy.jacTransport, FreeMatrix::Identity()); BOOST_CHECK_EQUAL(slsStateCopy.derivative, FreeVector::Zero()); BOOST_CHECK(slsStateCopy.covTransport); @@ -254,20 +233,19 @@ BOOST_AUTO_TEST_CASE(straight_line_stepper_test) { freeParams.template segment<3>(eFreeDir0).normalized(), 1e-6); CHECK_CLOSE_ABS(sls.absoluteMomentum(slsStateCopy), std::abs(1. / freeParams[eFreeQOverP]), 1e-6); - CHECK_CLOSE_ABS(sls.charge(slsStateCopy), -sls.charge(ps.stepping), 1e-6); + CHECK_CLOSE_ABS(sls.charge(slsStateCopy), -sls.charge(slsState), 1e-6); CHECK_CLOSE_ABS(sls.time(slsStateCopy), freeParams[eFreeTime], 1e-6); BOOST_CHECK_EQUAL(slsStateCopy.pathAccumulated, 0.); BOOST_CHECK_EQUAL(slsStateCopy.stepSize.value(), navDir * stepSize2); - BOOST_CHECK_EQUAL(slsStateCopy.previousStepSize, - ps.stepping.previousStepSize); + BOOST_CHECK_EQUAL(slsStateCopy.previousStepSize, slsState.previousStepSize); // Reset all possible parameters except the step size - slsStateCopy = ps.stepping; + slsStateCopy = slsState; sls.resetState(slsStateCopy, cp2.parameters(), *cp2.covariance(), cp2.referenceSurface()); // Test all components BOOST_CHECK_NE(slsStateCopy.jacToGlobal, BoundToFreeMatrix::Zero()); - BOOST_CHECK_NE(slsStateCopy.jacToGlobal, ps.stepping.jacToGlobal); + BOOST_CHECK_NE(slsStateCopy.jacToGlobal, slsState.jacToGlobal); BOOST_CHECK_EQUAL(slsStateCopy.jacTransport, FreeMatrix::Identity()); BOOST_CHECK_EQUAL(slsStateCopy.derivative, FreeVector::Zero()); BOOST_CHECK(slsStateCopy.covTransport); @@ -278,21 +256,20 @@ BOOST_AUTO_TEST_CASE(straight_line_stepper_test) { freeParams.template segment<3>(eFreeDir0), 1e-6); CHECK_CLOSE_ABS(sls.absoluteMomentum(slsStateCopy), std::abs(1. / freeParams[eFreeQOverP]), 1e-6); - CHECK_CLOSE_ABS(sls.charge(slsStateCopy), -sls.charge(ps.stepping), 1e-6); + CHECK_CLOSE_ABS(sls.charge(slsStateCopy), -sls.charge(slsState), 1e-6); CHECK_CLOSE_ABS(sls.time(slsStateCopy), freeParams[eFreeTime], 1e-6); BOOST_CHECK_EQUAL(slsStateCopy.pathAccumulated, 0.); BOOST_CHECK_EQUAL(slsStateCopy.stepSize.value(), std::numeric_limits::max()); - BOOST_CHECK_EQUAL(slsStateCopy.previousStepSize, - ps.stepping.previousStepSize); + BOOST_CHECK_EQUAL(slsStateCopy.previousStepSize, slsState.previousStepSize); // Reset the least amount of parameters - slsStateCopy = ps.stepping; + slsStateCopy = slsState; sls.resetState(slsStateCopy, cp2.parameters(), *cp2.covariance(), cp2.referenceSurface()); // Test all components BOOST_CHECK_NE(slsStateCopy.jacToGlobal, BoundToFreeMatrix::Zero()); - BOOST_CHECK_NE(slsStateCopy.jacToGlobal, ps.stepping.jacToGlobal); + BOOST_CHECK_NE(slsStateCopy.jacToGlobal, slsState.jacToGlobal); BOOST_CHECK_EQUAL(slsStateCopy.jacTransport, FreeMatrix::Identity()); BOOST_CHECK_EQUAL(slsStateCopy.derivative, FreeVector::Zero()); BOOST_CHECK(slsStateCopy.covTransport); @@ -303,13 +280,12 @@ BOOST_AUTO_TEST_CASE(straight_line_stepper_test) { freeParams.template segment<3>(eFreeDir0).normalized(), 1e-6); CHECK_CLOSE_ABS(sls.absoluteMomentum(slsStateCopy), std::abs(1. / freeParams[eFreeQOverP]), 1e-6); - CHECK_CLOSE_ABS(sls.charge(slsStateCopy), -sls.charge(ps.stepping), 1e-6); + CHECK_CLOSE_ABS(sls.charge(slsStateCopy), -sls.charge(slsState), 1e-6); CHECK_CLOSE_ABS(sls.time(slsStateCopy), freeParams[eFreeTime], 1e-6); BOOST_CHECK_EQUAL(slsStateCopy.pathAccumulated, 0.); BOOST_CHECK_EQUAL(slsStateCopy.stepSize.value(), std::numeric_limits::max()); - BOOST_CHECK_EQUAL(slsStateCopy.previousStepSize, - ps.stepping.previousStepSize); + BOOST_CHECK_EQUAL(slsStateCopy.previousStepSize, slsState.previousStepSize); /// Repeat with surface related methods auto plane = CurvilinearSurface(pos, dir).planeSurface(); diff --git a/Tests/UnitTests/Core/Propagator/SympyStepperTests.cpp b/Tests/UnitTests/Core/Propagator/SympyStepperTests.cpp index 1efc0dfa0ba..358c86e7556 100644 --- a/Tests/UnitTests/Core/Propagator/SympyStepperTests.cpp +++ b/Tests/UnitTests/Core/Propagator/SympyStepperTests.cpp @@ -30,7 +30,6 @@ #include "Acts/Utilities/Logger.hpp" #include "Acts/Utilities/Result.hpp" -#include #include #include #include @@ -53,31 +52,6 @@ static constexpr auto eps = 3 * std::numeric_limits::epsilon(); GeometryContext tgContext = GeometryContext(); MagneticFieldContext mfContext = MagneticFieldContext(); -/// @brief Simplified propagator state -template -struct PropState { - /// @brief Constructor - PropState(Direction direction, stepper_state_t sState) - : stepping(std::move(sState)) { - options.direction = direction; - } - /// State of the sympy stepper - stepper_state_t stepping; - /// Propagator options which only carry the relevant components - struct { - Direction direction = Direction::Forward(); - struct { - double stepTolerance = 1e-4; - double stepSizeCutOff = 0.; - unsigned int maxRungeKuttaStepTrials = 10000; - } stepping; - } options; -}; - -struct MockNavigator {}; - -static constexpr MockNavigator mockNavigator; - /// @brief Aborter for the case that a particle leaves the detector or reaches /// a custom made threshold. /// @@ -277,27 +251,26 @@ BOOST_AUTO_TEST_CASE(sympy_stepper_test) { // Perform a step without and with covariance transport esState.cov = cov; - PropState ps(navDir, std::move(esState)); - - ps.stepping.covTransport = false; - es.step(ps, mockNavigator).value(); - CHECK_CLOSE_COVARIANCE(ps.stepping.cov, cov, eps); - BOOST_CHECK_NE(es.position(ps.stepping).norm(), newPos.norm()); - BOOST_CHECK_NE(es.direction(ps.stepping), newMom.normalized()); - BOOST_CHECK_EQUAL(es.charge(ps.stepping), charge); - BOOST_CHECK_LT(es.time(ps.stepping), newTime); - BOOST_CHECK_EQUAL(ps.stepping.derivative, FreeVector::Zero()); - BOOST_CHECK_EQUAL(ps.stepping.jacTransport, FreeMatrix::Identity()); - - ps.stepping.covTransport = true; - es.step(ps, mockNavigator).value(); - CHECK_CLOSE_COVARIANCE(ps.stepping.cov, cov, eps); - BOOST_CHECK_NE(es.position(ps.stepping).norm(), newPos.norm()); - BOOST_CHECK_NE(es.direction(ps.stepping), newMom.normalized()); - BOOST_CHECK_EQUAL(es.charge(ps.stepping), charge); - BOOST_CHECK_LT(es.time(ps.stepping), newTime); - BOOST_CHECK_NE(ps.stepping.derivative, FreeVector::Zero()); - BOOST_CHECK_NE(ps.stepping.jacTransport, FreeMatrix::Identity()); + + esState.covTransport = false; + es.step(esState, Direction::Forward(), nullptr).value(); + CHECK_CLOSE_COVARIANCE(esState.cov, cov, eps); + BOOST_CHECK_NE(es.position(esState).norm(), newPos.norm()); + BOOST_CHECK_NE(es.direction(esState), newMom.normalized()); + BOOST_CHECK_EQUAL(es.charge(esState), charge); + BOOST_CHECK_LT(es.time(esState), newTime); + BOOST_CHECK_EQUAL(esState.derivative, FreeVector::Zero()); + BOOST_CHECK_EQUAL(esState.jacTransport, FreeMatrix::Identity()); + + esState.covTransport = true; + es.step(esState, Direction::Forward(), nullptr).value(); + CHECK_CLOSE_COVARIANCE(esState.cov, cov, eps); + BOOST_CHECK_NE(es.position(esState).norm(), newPos.norm()); + BOOST_CHECK_NE(es.direction(esState), newMom.normalized()); + BOOST_CHECK_EQUAL(es.charge(esState), charge); + BOOST_CHECK_LT(es.time(esState), newTime); + BOOST_CHECK_NE(esState.derivative, FreeVector::Zero()); + BOOST_CHECK_NE(esState.jacTransport, FreeMatrix::Identity()); /// Test the state reset // Construct the parameters @@ -337,13 +310,13 @@ BOOST_AUTO_TEST_CASE(sympy_stepper_test) { }; // Reset all possible parameters - SympyStepper::State esStateCopy = copyState(*bField, ps.stepping); + SympyStepper::State esStateCopy = copyState(*bField, esState); BOOST_CHECK(cp2.covariance().has_value()); es.resetState(esStateCopy, cp2.parameters(), *cp2.covariance(), cp2.referenceSurface(), stepSize2); // Test all components BOOST_CHECK_NE(esStateCopy.jacToGlobal, BoundToFreeMatrix::Zero()); - BOOST_CHECK_NE(esStateCopy.jacToGlobal, ps.stepping.jacToGlobal); + BOOST_CHECK_NE(esStateCopy.jacToGlobal, esState.jacToGlobal); BOOST_CHECK_EQUAL(esStateCopy.jacTransport, FreeMatrix::Identity()); BOOST_CHECK_EQUAL(esStateCopy.derivative, FreeVector::Zero()); BOOST_CHECK(esStateCopy.covTransport); @@ -354,19 +327,19 @@ BOOST_AUTO_TEST_CASE(sympy_stepper_test) { freeParams.template segment<3>(eFreeDir0).normalized()); BOOST_CHECK_EQUAL(es.absoluteMomentum(esStateCopy), std::abs(1. / freeParams[eFreeQOverP])); - BOOST_CHECK_EQUAL(es.charge(esStateCopy), -es.charge(ps.stepping)); + BOOST_CHECK_EQUAL(es.charge(esStateCopy), -es.charge(esState)); BOOST_CHECK_EQUAL(es.time(esStateCopy), freeParams[eFreeTime]); BOOST_CHECK_EQUAL(esStateCopy.pathAccumulated, 0.); BOOST_CHECK_EQUAL(esStateCopy.stepSize.value(), navDir * stepSize2); - BOOST_CHECK_EQUAL(esStateCopy.previousStepSize, ps.stepping.previousStepSize); + BOOST_CHECK_EQUAL(esStateCopy.previousStepSize, esState.previousStepSize); // Reset all possible parameters except the step size - esStateCopy = copyState(*bField, ps.stepping); + esStateCopy = copyState(*bField, esState); es.resetState(esStateCopy, cp2.parameters(), *cp2.covariance(), cp2.referenceSurface()); // Test all components BOOST_CHECK_NE(esStateCopy.jacToGlobal, BoundToFreeMatrix::Zero()); - BOOST_CHECK_NE(esStateCopy.jacToGlobal, ps.stepping.jacToGlobal); + BOOST_CHECK_NE(esStateCopy.jacToGlobal, esState.jacToGlobal); BOOST_CHECK_EQUAL(esStateCopy.jacTransport, FreeMatrix::Identity()); BOOST_CHECK_EQUAL(esStateCopy.derivative, FreeVector::Zero()); BOOST_CHECK(esStateCopy.covTransport); @@ -377,20 +350,20 @@ BOOST_AUTO_TEST_CASE(sympy_stepper_test) { freeParams.template segment<3>(eFreeDir0).normalized()); BOOST_CHECK_EQUAL(es.absoluteMomentum(esStateCopy), std::abs(1. / freeParams[eFreeQOverP])); - BOOST_CHECK_EQUAL(es.charge(esStateCopy), -es.charge(ps.stepping)); + BOOST_CHECK_EQUAL(es.charge(esStateCopy), -es.charge(esState)); BOOST_CHECK_EQUAL(es.time(esStateCopy), freeParams[eFreeTime]); BOOST_CHECK_EQUAL(esStateCopy.pathAccumulated, 0.); BOOST_CHECK_EQUAL(esStateCopy.stepSize.value(), std::numeric_limits::max()); - BOOST_CHECK_EQUAL(esStateCopy.previousStepSize, ps.stepping.previousStepSize); + BOOST_CHECK_EQUAL(esStateCopy.previousStepSize, esState.previousStepSize); // Reset the least amount of parameters - esStateCopy = copyState(*bField, ps.stepping); + esStateCopy = copyState(*bField, esState); es.resetState(esStateCopy, cp2.parameters(), *cp2.covariance(), cp2.referenceSurface()); // Test all components BOOST_CHECK_NE(esStateCopy.jacToGlobal, BoundToFreeMatrix::Zero()); - BOOST_CHECK_NE(esStateCopy.jacToGlobal, ps.stepping.jacToGlobal); + BOOST_CHECK_NE(esStateCopy.jacToGlobal, esState.jacToGlobal); BOOST_CHECK_EQUAL(esStateCopy.jacTransport, FreeMatrix::Identity()); BOOST_CHECK_EQUAL(esStateCopy.derivative, FreeVector::Zero()); BOOST_CHECK(esStateCopy.covTransport); @@ -401,12 +374,12 @@ BOOST_AUTO_TEST_CASE(sympy_stepper_test) { freeParams.template segment<3>(eFreeDir0).normalized()); BOOST_CHECK_EQUAL(es.absoluteMomentum(esStateCopy), std::abs(1. / freeParams[eFreeQOverP])); - BOOST_CHECK_EQUAL(es.charge(esStateCopy), -es.charge(ps.stepping)); + BOOST_CHECK_EQUAL(es.charge(esStateCopy), -es.charge(esState)); BOOST_CHECK_EQUAL(es.time(esStateCopy), freeParams[eFreeTime]); BOOST_CHECK_EQUAL(esStateCopy.pathAccumulated, 0.); BOOST_CHECK_EQUAL(esStateCopy.stepSize.value(), std::numeric_limits::max()); - BOOST_CHECK_EQUAL(esStateCopy.previousStepSize, ps.stepping.previousStepSize); + BOOST_CHECK_EQUAL(esStateCopy.previousStepSize, esState.previousStepSize); /// Repeat with surface related methods auto plane = CurvilinearSurface(pos, dir.normalized()).planeSurface(); @@ -480,9 +453,9 @@ BOOST_AUTO_TEST_CASE(sympy_stepper_test) { CHECK_CLOSE_COVARIANCE(esState.cov, Covariance(2. * cov), eps); // Test a case where no step size adjustment is required - ps.options.stepping.stepTolerance = 2. * 4.4258e+09; + esState.options.stepTolerance = 2. * 4.4258e+09; double h0 = esState.stepSize.value(); - es.step(ps, mockNavigator); + es.step(esState, Direction::Forward(), nullptr); CHECK_CLOSE_ABS(h0, esState.stepSize.value(), eps); } From aee976ba35e43ab7540d42da47ff202149b092eb Mon Sep 17 00:00:00 2001 From: Andreas Stefl Date: Mon, 20 Jan 2025 19:09:22 +0100 Subject: [PATCH 02/11] fix formatting --- Core/include/Acts/Propagator/EigenStepper.ipp | 10 ++++------ Core/include/Acts/Propagator/MultiEigenStepperLoop.ipp | 5 ++--- 2 files changed, 6 insertions(+), 9 deletions(-) diff --git a/Core/include/Acts/Propagator/EigenStepper.ipp b/Core/include/Acts/Propagator/EigenStepper.ipp index 6f66ea66302..08b853c1223 100644 --- a/Core/include/Acts/Propagator/EigenStepper.ipp +++ b/Core/include/Acts/Propagator/EigenStepper.ipp @@ -20,9 +20,8 @@ Acts::EigenStepper::EigenStepper( : m_bField(std::move(bField)) {} template -auto Acts::EigenStepper::makeState(const Options& options, - const BoundTrackParameters& par) const - -> State { +auto Acts::EigenStepper::makeState( + const Options& options, const BoundTrackParameters& par) const -> State { State state{options, m_bField->makeCache(options.magFieldContext)}; state.particleHypothesis = par.particleHypothesis(); @@ -120,9 +119,8 @@ bool Acts::EigenStepper::prepareCurvilinearState(State& state) const { } template -auto Acts::EigenStepper::curvilinearState(State& state, - bool transportCov) const - -> CurvilinearState { +auto Acts::EigenStepper::curvilinearState( + State& state, bool transportCov) const -> CurvilinearState { return detail::curvilinearState( state.cov, state.jacobian, state.jacTransport, state.derivative, state.jacToGlobal, state.pars, state.particleHypothesis, diff --git a/Core/include/Acts/Propagator/MultiEigenStepperLoop.ipp b/Core/include/Acts/Propagator/MultiEigenStepperLoop.ipp index c3153b2f278..44b534a69ec 100644 --- a/Core/include/Acts/Propagator/MultiEigenStepperLoop.ipp +++ b/Core/include/Acts/Propagator/MultiEigenStepperLoop.ipp @@ -64,9 +64,8 @@ auto MultiEigenStepperLoop::boundState( } template -auto MultiEigenStepperLoop::curvilinearState(State& state, - bool transportCov) const - -> CurvilinearState { +auto MultiEigenStepperLoop::curvilinearState( + State& state, bool transportCov) const -> CurvilinearState { assert(!state.components.empty()); std::vector> From e496bf09e5f32c5d7e2a652e358aad1230eeba80 Mon Sep 17 00:00:00 2001 From: Andreas Stefl Date: Mon, 20 Jan 2025 19:47:16 +0100 Subject: [PATCH 03/11] fix stepper tests --- Core/src/Propagator/SympyStepper.cpp | 2 +- Tests/UnitTests/Core/Propagator/AtlasStepperTests.cpp | 6 +++--- Tests/UnitTests/Core/Propagator/EigenStepperTests.cpp | 4 ++-- .../UnitTests/Core/Propagator/StraightLineStepperTests.cpp | 4 ++-- Tests/UnitTests/Core/Propagator/SympyStepperTests.cpp | 4 ++-- 5 files changed, 10 insertions(+), 10 deletions(-) diff --git a/Core/src/Propagator/SympyStepper.cpp b/Core/src/Propagator/SympyStepper.cpp index 4fb155a39fb..0d2253bec20 100644 --- a/Core/src/Propagator/SympyStepper.cpp +++ b/Core/src/Propagator/SympyStepper.cpp @@ -137,7 +137,7 @@ Result SympyStepper::step(State& state, Direction propDir, const IVolumeMaterial* material) const { (void)material; return stepImpl(state, propDir, state.options.stepTolerance, - state.options.maxStepSize, + state.options.stepSizeCutOff, state.options.maxRungeKuttaStepTrials); } diff --git a/Tests/UnitTests/Core/Propagator/AtlasStepperTests.cpp b/Tests/UnitTests/Core/Propagator/AtlasStepperTests.cpp index 93af62cfa95..d2476ba2f99 100644 --- a/Tests/UnitTests/Core/Propagator/AtlasStepperTests.cpp +++ b/Tests/UnitTests/Core/Propagator/AtlasStepperTests.cpp @@ -294,7 +294,7 @@ BOOST_AUTO_TEST_CASE(Step) { state.covTransport = false; // ensure step does not result in an error - auto res = stepper.step(state, Direction::Forward(), nullptr); + auto res = stepper.step(state, Direction::Backward(), nullptr); BOOST_CHECK(res.ok()); // extract the actual step size @@ -331,7 +331,7 @@ BOOST_AUTO_TEST_CASE(StepWithCovariance) { state.covTransport = true; // ensure step does not result in an error - auto res = stepper.step(state, Direction::Forward(), nullptr); + auto res = stepper.step(state, Direction::Backward(), nullptr); BOOST_CHECK(res.ok()); // extract the actual step size @@ -371,7 +371,7 @@ BOOST_AUTO_TEST_CASE(Reset) { state.covTransport = true; // ensure step does not result in an error - stepper.step(state, Direction::Forward(), nullptr); + stepper.step(state, Direction::Backward(), nullptr); // Construct the parameters Vector3 newPos(1.5, -2.5, 3.5); diff --git a/Tests/UnitTests/Core/Propagator/EigenStepperTests.cpp b/Tests/UnitTests/Core/Propagator/EigenStepperTests.cpp index 909dc36b4da..eb64438a4c7 100644 --- a/Tests/UnitTests/Core/Propagator/EigenStepperTests.cpp +++ b/Tests/UnitTests/Core/Propagator/EigenStepperTests.cpp @@ -279,7 +279,7 @@ BOOST_AUTO_TEST_CASE(eigen_stepper_test) { esState.cov = cov; esState.covTransport = false; - es.step(esState, Direction::Forward(), nullptr).value(); + es.step(esState, navDir, nullptr).value(); CHECK_CLOSE_COVARIANCE(esState.cov, cov, eps); BOOST_CHECK_NE(es.position(esState).norm(), newPos.norm()); BOOST_CHECK_NE(es.direction(esState), newMom.normalized()); @@ -289,7 +289,7 @@ BOOST_AUTO_TEST_CASE(eigen_stepper_test) { BOOST_CHECK_EQUAL(esState.jacTransport, FreeMatrix::Identity()); esState.covTransport = true; - es.step(esState, Direction::Forward(), nullptr).value(); + es.step(esState, navDir, nullptr).value(); CHECK_CLOSE_COVARIANCE(esState.cov, cov, eps); BOOST_CHECK_NE(es.position(esState).norm(), newPos.norm()); BOOST_CHECK_NE(es.direction(esState), newMom.normalized()); diff --git a/Tests/UnitTests/Core/Propagator/StraightLineStepperTests.cpp b/Tests/UnitTests/Core/Propagator/StraightLineStepperTests.cpp index 8ace9416308..ebda0e399a9 100644 --- a/Tests/UnitTests/Core/Propagator/StraightLineStepperTests.cpp +++ b/Tests/UnitTests/Core/Propagator/StraightLineStepperTests.cpp @@ -174,7 +174,7 @@ BOOST_AUTO_TEST_CASE(straight_line_stepper_test) { slsState.cov = cov; slsState.covTransport = false; - double h = sls.step(slsState, Direction::Forward(), nullptr).value(); + double h = sls.step(slsState, navDir, nullptr).value(); BOOST_CHECK_EQUAL(slsState.stepSize.value(), stepSize); BOOST_CHECK_EQUAL(slsState.stepSize.value(), h * navDir); CHECK_CLOSE_COVARIANCE(slsState.cov, cov, 1e-6); @@ -187,7 +187,7 @@ BOOST_AUTO_TEST_CASE(straight_line_stepper_test) { BOOST_CHECK_EQUAL(slsState.jacTransport, FreeMatrix::Identity()); slsState.covTransport = true; - double h2 = sls.step(slsState, Direction::Forward(), nullptr).value(); + double h2 = sls.step(slsState, navDir, nullptr).value(); BOOST_CHECK_EQUAL(slsState.stepSize.value(), stepSize); BOOST_CHECK_EQUAL(h2, h); CHECK_CLOSE_COVARIANCE(slsState.cov, cov, 1e-6); diff --git a/Tests/UnitTests/Core/Propagator/SympyStepperTests.cpp b/Tests/UnitTests/Core/Propagator/SympyStepperTests.cpp index 358c86e7556..9c468188685 100644 --- a/Tests/UnitTests/Core/Propagator/SympyStepperTests.cpp +++ b/Tests/UnitTests/Core/Propagator/SympyStepperTests.cpp @@ -253,7 +253,7 @@ BOOST_AUTO_TEST_CASE(sympy_stepper_test) { esState.cov = cov; esState.covTransport = false; - es.step(esState, Direction::Forward(), nullptr).value(); + es.step(esState, navDir, nullptr).value(); CHECK_CLOSE_COVARIANCE(esState.cov, cov, eps); BOOST_CHECK_NE(es.position(esState).norm(), newPos.norm()); BOOST_CHECK_NE(es.direction(esState), newMom.normalized()); @@ -263,7 +263,7 @@ BOOST_AUTO_TEST_CASE(sympy_stepper_test) { BOOST_CHECK_EQUAL(esState.jacTransport, FreeMatrix::Identity()); esState.covTransport = true; - es.step(esState, Direction::Forward(), nullptr).value(); + es.step(esState, navDir, nullptr).value(); CHECK_CLOSE_COVARIANCE(esState.cov, cov, eps); BOOST_CHECK_NE(es.position(esState).norm(), newPos.norm()); BOOST_CHECK_NE(es.direction(esState), newMom.normalized()); From 1727f9ce7109ba2629c8c2f508b632847a2cb716 Mon Sep 17 00:00:00 2001 From: Andreas Stefl Date: Mon, 20 Jan 2025 19:57:19 +0100 Subject: [PATCH 04/11] fix doc --- Core/include/Acts/Propagator/AtlasStepper.hpp | 3 +++ Core/include/Acts/Propagator/EigenStepper.hpp | 5 ++++- .../Propagator/EigenStepperDefaultExtension.hpp | 3 +++ .../Propagator/EigenStepperDenseExtension.hpp | 3 +++ .../Acts/Propagator/MultiEigenStepperLoop.hpp | 3 +++ .../Acts/Propagator/StraightLineStepper.hpp | 15 ++++++++------- Core/include/Acts/Propagator/SympyStepper.hpp | 5 ++++- 7 files changed, 28 insertions(+), 9 deletions(-) diff --git a/Core/include/Acts/Propagator/AtlasStepper.hpp b/Core/include/Acts/Propagator/AtlasStepper.hpp index c698e17dd52..4dfe45da30c 100644 --- a/Core/include/Acts/Propagator/AtlasStepper.hpp +++ b/Core/include/Acts/Propagator/AtlasStepper.hpp @@ -1153,6 +1153,9 @@ class AtlasStepper { /// Perform the actual step on the state /// /// @param state is the provided stepper state (caller keeps thread locality) + /// @param propDir is the direction of propagation + /// @param material is the material properties + /// @return the result of the step Result step(State& state, Direction propDir, const IVolumeMaterial* material) const { (void)material; diff --git a/Core/include/Acts/Propagator/EigenStepper.hpp b/Core/include/Acts/Propagator/EigenStepper.hpp index e28f45d3063..aaccdf87f4a 100644 --- a/Core/include/Acts/Propagator/EigenStepper.hpp +++ b/Core/include/Acts/Propagator/EigenStepper.hpp @@ -397,7 +397,10 @@ class EigenStepper { /// Perform a Runge-Kutta track parameter propagation step /// /// @param [in,out] state The state of the stepper - /// @note The state contains the desired step size. It can be negative during + /// @param propDir is the direction of propagation + /// @param material is the material properties + /// @return the result of the step + /// @note The state contains the desired step size. It can be negative during /// backwards track propagation, and since we're using an adaptive /// algorithm, it can be modified by the stepper class during /// propagation. diff --git a/Core/include/Acts/Propagator/EigenStepperDefaultExtension.hpp b/Core/include/Acts/Propagator/EigenStepperDefaultExtension.hpp index cfc69def7f5..549ef05f4a4 100644 --- a/Core/include/Acts/Propagator/EigenStepperDefaultExtension.hpp +++ b/Core/include/Acts/Propagator/EigenStepperDefaultExtension.hpp @@ -29,6 +29,7 @@ struct EigenStepperDefaultExtension { /// /// @param [in] state State of the stepper /// @param [in] stepper Stepper of the propagation + /// @param [in] volumeMaterial Material of the volume /// @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 @@ -64,6 +65,7 @@ struct EigenStepperDefaultExtension { /// /// @param [in] state State of the stepper /// @param [in] stepper Stepper of the propagation + /// @param [in] volumeMaterial Material of the volume /// @param [in] h Step size /// /// @return Boolean flag if the calculation is valid @@ -84,6 +86,7 @@ struct EigenStepperDefaultExtension { /// /// @param [in] state State of the stepper /// @param [in] stepper Stepper of the propagation + /// @param [in] volumeMaterial Material of the volume /// @param [in] h Step size /// @param [out] D Transport matrix /// diff --git a/Core/include/Acts/Propagator/EigenStepperDenseExtension.hpp b/Core/include/Acts/Propagator/EigenStepperDenseExtension.hpp index 1aafd3967f9..4ca85fea628 100644 --- a/Core/include/Acts/Propagator/EigenStepperDenseExtension.hpp +++ b/Core/include/Acts/Propagator/EigenStepperDenseExtension.hpp @@ -62,6 +62,7 @@ struct EigenStepperDenseExtension { /// /// @param [in] state State of the stepper /// @param [in] stepper Stepper of the propagator + /// @param [in] volumeMaterial Material of the volume /// @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 @@ -127,6 +128,7 @@ struct EigenStepperDenseExtension { /// /// @param [in] state State of the stepper /// @param [in] stepper Stepper of the propagator + /// @param [in] volumeMaterial Material of the volume /// @param [in] h Step size /// /// @return Boolean flag if the calculation is valid @@ -176,6 +178,7 @@ struct EigenStepperDenseExtension { /// /// @param [in] state State of the stepper /// @param [in] stepper Stepper of the propagator + /// @param [in] volumeMaterial Material of the volume /// @param [in] h Step size /// @param [out] D Transport matrix /// diff --git a/Core/include/Acts/Propagator/MultiEigenStepperLoop.hpp b/Core/include/Acts/Propagator/MultiEigenStepperLoop.hpp index a05e8281e33..b0f58612785 100644 --- a/Core/include/Acts/Propagator/MultiEigenStepperLoop.hpp +++ b/Core/include/Acts/Propagator/MultiEigenStepperLoop.hpp @@ -734,6 +734,9 @@ class MultiEigenStepperLoop : public EigenStepper { /// Perform a Runge-Kutta track parameter propagation step /// /// @param [in,out] state The state of the stepper + /// @param propDir is the direction of propagation + /// @param material is the material properties + /// @return the result of the step /// /// The state contains the desired step size. It can be negative during /// backwards track propagation, and since we're using an adaptive diff --git a/Core/include/Acts/Propagator/StraightLineStepper.hpp b/Core/include/Acts/Propagator/StraightLineStepper.hpp index d1e96f30e6b..686e5a828bd 100644 --- a/Core/include/Acts/Propagator/StraightLineStepper.hpp +++ b/Core/include/Acts/Propagator/StraightLineStepper.hpp @@ -400,13 +400,14 @@ class StraightLineStepper { /// Perform a straight line propagation step /// /// @param [in,out] state The state of the stepper - /// parameters that are being propagated. - /// The state contains the desired step size, - /// it can be negative during backwards track propagation, - /// and since we're using an adaptive algorithm, it can - /// be modified by the stepper class during propagation. - /// - /// @return the step size taken + /// @param propDir is the direction of propagation + /// @param material is the material properties + /// @return the result of the step + /// + /// The state contains the desired step size, + /// it can be negative during backwards track propagation, + /// and since we're using an adaptive algorithm, it can + /// be modified by the stepper class during propagation. Result step(State& state, Direction propDir, const IVolumeMaterial* material) const { (void)material; diff --git a/Core/include/Acts/Propagator/SympyStepper.hpp b/Core/include/Acts/Propagator/SympyStepper.hpp index 8fa694d6176..cf5ed8eae9c 100644 --- a/Core/include/Acts/Propagator/SympyStepper.hpp +++ b/Core/include/Acts/Propagator/SympyStepper.hpp @@ -364,7 +364,10 @@ class SympyStepper { /// Perform a Runge-Kutta track parameter propagation step /// /// @param [in,out] state State of the stepper - /// @note The state contains the desired step size. It can be negative during + /// @param propDir is the direction of propagation + /// @param material is the material properties + /// @return the result of the step + /// @note The state contains the desired step size. It can be negative during /// backwards track propagation, and since we're using an adaptive /// algorithm, it can be modified by the stepper class during /// propagation. From 6ba53fb03832532c4e9846dd2059db52ac4b0960 Mon Sep 17 00:00:00 2001 From: Andreas Stefl Date: Mon, 20 Jan 2025 20:51:29 +0100 Subject: [PATCH 05/11] temporary switch off tests --- Tests/UnitTests/Core/Propagator/CMakeLists.txt | 2 +- Tests/UnitTests/Core/Vertexing/CMakeLists.txt | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/Tests/UnitTests/Core/Propagator/CMakeLists.txt b/Tests/UnitTests/Core/Propagator/CMakeLists.txt index 34f3e88f23f..2c4a7ed1b2f 100644 --- a/Tests/UnitTests/Core/Propagator/CMakeLists.txt +++ b/Tests/UnitTests/Core/Propagator/CMakeLists.txt @@ -14,5 +14,5 @@ add_unittest(Propagator PropagatorTests.cpp) add_unittest(EigenStepper EigenStepperTests.cpp) add_unittest(StraightLineStepper StraightLineStepperTests.cpp) add_unittest(VolumeMaterialInteraction VolumeMaterialInteractionTests.cpp) -add_unittest(BoundToCurvilinearConversionTests BoundToCurvilinearConversionTests.cpp) +#add_unittest(BoundToCurvilinearConversionTests BoundToCurvilinearConversionTests.cpp) add_unittest(SympyStepper SympyStepperTests.cpp) diff --git a/Tests/UnitTests/Core/Vertexing/CMakeLists.txt b/Tests/UnitTests/Core/Vertexing/CMakeLists.txt index ab7a2873500..c12ca43deb7 100644 --- a/Tests/UnitTests/Core/Vertexing/CMakeLists.txt +++ b/Tests/UnitTests/Core/Vertexing/CMakeLists.txt @@ -1,5 +1,5 @@ add_unittest(FullBilloirVertexFitter FullBilloirVertexFitterTests.cpp) -add_unittest(ImpactPointEstimator ImpactPointEstimatorTests.cpp) +#add_unittest(ImpactPointEstimator ImpactPointEstimatorTests.cpp) add_unittest(IterativeVertexFinder IterativeVertexFinderTests.cpp) add_unittest(KalmanVertexUpdater KalmanVertexUpdaterTests.cpp) add_unittest(LinearizedTrackFactory LinearizedTrackFactoryTests.cpp) From afaadc8dd4eacebc70d8bace953c3bc1f42b47e9 Mon Sep 17 00:00:00 2001 From: Andreas Stefl Date: Mon, 20 Jan 2025 21:15:07 +0100 Subject: [PATCH 06/11] try workaround compiler bug --- Core/include/Acts/Propagator/MultiEigenStepperLoop.ipp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/Core/include/Acts/Propagator/MultiEigenStepperLoop.ipp b/Core/include/Acts/Propagator/MultiEigenStepperLoop.ipp index 44b534a69ec..0d9cb418d38 100644 --- a/Core/include/Acts/Propagator/MultiEigenStepperLoop.ipp +++ b/Core/include/Acts/Propagator/MultiEigenStepperLoop.ipp @@ -155,7 +155,8 @@ Result MultiEigenStepperLoop::step( // Lambda that performs the step for a component and returns false if the step // went ok and true if there was an error - auto errorInStep = [&](auto& component) { + auto errorInStep = [this, &results, propDir, material, &accumulatedPathLength, + &errorSteps, &reweightNecessary](auto& component) { if (component.status == Status::onSurface) { // We need to add these, so the propagation does not fail if we have only // components on surfaces and failing states From 3be0f78dd9941f500824a21134d0d481855aebf6 Mon Sep 17 00:00:00 2001 From: Andreas Stefl Date: Mon, 20 Jan 2025 23:01:06 +0100 Subject: [PATCH 07/11] fix eigen stepper --- Core/include/Acts/Propagator/EigenStepperDefaultExtension.hpp | 2 +- Tests/UnitTests/Core/Propagator/CMakeLists.txt | 2 +- Tests/UnitTests/Core/Vertexing/CMakeLists.txt | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/Core/include/Acts/Propagator/EigenStepperDefaultExtension.hpp b/Core/include/Acts/Propagator/EigenStepperDefaultExtension.hpp index 549ef05f4a4..35991286d17 100644 --- a/Core/include/Acts/Propagator/EigenStepperDefaultExtension.hpp +++ b/Core/include/Acts/Propagator/EigenStepperDefaultExtension.hpp @@ -92,7 +92,7 @@ struct EigenStepperDefaultExtension { /// /// @return Boolean flag if the calculation is valid template - bool finalize(stepper_t::State state, const stepper_t& stepper, + bool finalize(stepper_t::State& state, const stepper_t& stepper, const IVolumeMaterial* volumeMaterial, const double h, FreeMatrix& D) const { (void)volumeMaterial; diff --git a/Tests/UnitTests/Core/Propagator/CMakeLists.txt b/Tests/UnitTests/Core/Propagator/CMakeLists.txt index 2c4a7ed1b2f..34f3e88f23f 100644 --- a/Tests/UnitTests/Core/Propagator/CMakeLists.txt +++ b/Tests/UnitTests/Core/Propagator/CMakeLists.txt @@ -14,5 +14,5 @@ add_unittest(Propagator PropagatorTests.cpp) add_unittest(EigenStepper EigenStepperTests.cpp) add_unittest(StraightLineStepper StraightLineStepperTests.cpp) add_unittest(VolumeMaterialInteraction VolumeMaterialInteractionTests.cpp) -#add_unittest(BoundToCurvilinearConversionTests BoundToCurvilinearConversionTests.cpp) +add_unittest(BoundToCurvilinearConversionTests BoundToCurvilinearConversionTests.cpp) add_unittest(SympyStepper SympyStepperTests.cpp) diff --git a/Tests/UnitTests/Core/Vertexing/CMakeLists.txt b/Tests/UnitTests/Core/Vertexing/CMakeLists.txt index c12ca43deb7..ab7a2873500 100644 --- a/Tests/UnitTests/Core/Vertexing/CMakeLists.txt +++ b/Tests/UnitTests/Core/Vertexing/CMakeLists.txt @@ -1,5 +1,5 @@ add_unittest(FullBilloirVertexFitter FullBilloirVertexFitterTests.cpp) -#add_unittest(ImpactPointEstimator ImpactPointEstimatorTests.cpp) +add_unittest(ImpactPointEstimator ImpactPointEstimatorTests.cpp) add_unittest(IterativeVertexFinder IterativeVertexFinderTests.cpp) add_unittest(KalmanVertexUpdater KalmanVertexUpdaterTests.cpp) add_unittest(LinearizedTrackFactory LinearizedTrackFactoryTests.cpp) From 1d63f03f52b9f53cab1aae56315b6d1b1a162364 Mon Sep 17 00:00:00 2001 From: Andreas Stefl Date: Mon, 20 Jan 2025 23:06:59 +0100 Subject: [PATCH 08/11] fix typename --- .../Propagator/EigenStepperDefaultExtension.hpp | 13 +++++++------ .../Propagator/EigenStepperDenseExtension.hpp | 17 +++++++++-------- .../Acts/Propagator/detail/SteppingHelper.hpp | 2 +- 3 files changed, 17 insertions(+), 15 deletions(-) diff --git a/Core/include/Acts/Propagator/EigenStepperDefaultExtension.hpp b/Core/include/Acts/Propagator/EigenStepperDefaultExtension.hpp index 35991286d17..8d14627f950 100644 --- a/Core/include/Acts/Propagator/EigenStepperDefaultExtension.hpp +++ b/Core/include/Acts/Propagator/EigenStepperDefaultExtension.hpp @@ -38,7 +38,7 @@ struct EigenStepperDefaultExtension { /// /// @return Boolean flag if the calculation is valid template - bool k(const stepper_t::State& state, const stepper_t& stepper, + bool k(const typename stepper_t::State& state, const stepper_t& stepper, const IVolumeMaterial* volumeMaterial, Vector3& knew, const Vector3& bField, std::array& kQoP, const double h = 0., const Vector3& kprev = Vector3::Zero()) @@ -70,7 +70,7 @@ struct EigenStepperDefaultExtension { /// /// @return Boolean flag if the calculation is valid template - bool finalize(stepper_t::State& state, const stepper_t& stepper, + bool finalize(typename stepper_t::State& state, const stepper_t& stepper, const IVolumeMaterial* volumeMaterial, const double h) const { (void)volumeMaterial; @@ -92,7 +92,7 @@ struct EigenStepperDefaultExtension { /// /// @return Boolean flag if the calculation is valid template - bool finalize(stepper_t::State& state, const stepper_t& stepper, + bool finalize(typename stepper_t::State& state, const stepper_t& stepper, const IVolumeMaterial* volumeMaterial, const double h, FreeMatrix& D) const { (void)volumeMaterial; @@ -110,7 +110,7 @@ struct EigenStepperDefaultExtension { /// @param [in] stepper Stepper of the propagation /// @param [in] h Step size template - void propagateTime(stepper_t::State& state, const stepper_t& stepper, + void propagateTime(typename stepper_t::State& state, const stepper_t& stepper, const double h) const { /// This evaluation is based on dt/ds = 1/v = 1/(beta * c) with the velocity /// v, the speed of light c and beta = v/c. This can be re-written as dt/ds @@ -135,8 +135,9 @@ struct EigenStepperDefaultExtension { /// /// @return Boolean flag if evaluation is valid template - bool transportMatrix(stepper_t::State& state, const stepper_t& stepper, - const double h, FreeMatrix& D) const { + bool transportMatrix(typename stepper_t::State& state, + const stepper_t& stepper, const double h, + FreeMatrix& D) const { /// The calculations are based on ATL-SOFT-PUB-2009-002. The update of the /// Jacobian matrix is requires only the calculation of eq. 17 and 18. /// Since the terms of eq. 18 are currently 0, this matrix is not needed diff --git a/Core/include/Acts/Propagator/EigenStepperDenseExtension.hpp b/Core/include/Acts/Propagator/EigenStepperDenseExtension.hpp index 4ca85fea628..6ecae8856b2 100644 --- a/Core/include/Acts/Propagator/EigenStepperDenseExtension.hpp +++ b/Core/include/Acts/Propagator/EigenStepperDenseExtension.hpp @@ -71,7 +71,7 @@ struct EigenStepperDenseExtension { /// /// @return Boolean flag if the calculation is valid template - bool k(const stepper_t::State& state, const stepper_t& stepper, + bool k(const typename stepper_t::State& state, const stepper_t& stepper, const IVolumeMaterial* volumeMaterial, Vector3& knew, const Vector3& bField, std::array& kQoP, const double h = 0., const Vector3& kprev = Vector3::Zero()) @@ -133,7 +133,7 @@ struct EigenStepperDenseExtension { /// /// @return Boolean flag if the calculation is valid template - bool finalize(stepper_t::State& state, const stepper_t& stepper, + bool finalize(typename stepper_t::State& state, const stepper_t& stepper, const IVolumeMaterial* volumeMaterial, const double h) const { if (volumeMaterial == nullptr) { return defaultExtension.finalize(state, stepper, volumeMaterial, h); @@ -184,7 +184,7 @@ struct EigenStepperDenseExtension { /// /// @return Boolean flag if the calculation is valid template - bool finalize(stepper_t::State& state, const stepper_t& stepper, + bool finalize(typename stepper_t::State& state, const stepper_t& stepper, const IVolumeMaterial* volumeMaterial, const double h, FreeMatrix& D) const { if (volumeMaterial == nullptr) { @@ -208,8 +208,9 @@ struct EigenStepperDenseExtension { /// /// @return Boolean flag if evaluation is valid template - bool transportMatrix(stepper_t::State& state, const stepper_t& stepper, - const double h, FreeMatrix& D) const { + bool transportMatrix(typename stepper_t::State& state, + const stepper_t& stepper, const double h, + FreeMatrix& D) const { /// The calculations are based on ATL-SOFT-PUB-2009-002. The update of the /// Jacobian matrix is requires only the calculation of eq. 17 and 18. /// Since the terms of eq. 18 are currently 0, this matrix is not needed @@ -348,7 +349,7 @@ struct EigenStepperDenseExtension { /// @param [in] state Deliverer of configurations /// @param [in] stepper Stepper of the propagator template - void initializeEnergyLoss(const stepper_t::State& state, + void initializeEnergyLoss(const typename stepper_t::State& state, const stepper_t& stepper) { const auto& particleHypothesis = stepper.particleHypothesis(state); float mass = particleHypothesis.mass(); @@ -404,8 +405,8 @@ struct EigenStepperDenseExtension { /// @param [in] i Index of the sub-step (1-3) template void updateEnergyLoss(const double mass, const double h, - const stepper_t::State& state, const stepper_t& stepper, - const int i) { + const typename stepper_t::State& state, + const stepper_t& stepper, const int i) { // Update parameters related to a changed momentum currentMomentum = initialMomentum + h * dPds[i - 1]; energy[i] = fastHypot(currentMomentum, mass); diff --git a/Core/include/Acts/Propagator/detail/SteppingHelper.hpp b/Core/include/Acts/Propagator/detail/SteppingHelper.hpp index 9c7c77720c5..97111e71d45 100644 --- a/Core/include/Acts/Propagator/detail/SteppingHelper.hpp +++ b/Core/include/Acts/Propagator/detail/SteppingHelper.hpp @@ -39,7 +39,7 @@ namespace Acts::detail { /// @param logger [in] A @c Logger instance template IntersectionStatus updateSingleSurfaceStatus( - const stepper_t& stepper, typename stepper_t::State& state, + const stepper_t& stepper, typename typename stepper_t::State& state, const Surface& surface, std::uint8_t index, Direction direction, const BoundaryTolerance& boundaryTolerance, double surfaceTolerance, ConstrainedStep::Type stype, const Logger& logger) { From 7c2993d64222a0832ad42529b43ce89d6a317230 Mon Sep 17 00:00:00 2001 From: Andreas Stefl Date: Mon, 20 Jan 2025 23:49:22 +0100 Subject: [PATCH 09/11] fix silly --- Core/include/Acts/Propagator/detail/SteppingHelper.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Core/include/Acts/Propagator/detail/SteppingHelper.hpp b/Core/include/Acts/Propagator/detail/SteppingHelper.hpp index 97111e71d45..9c7c77720c5 100644 --- a/Core/include/Acts/Propagator/detail/SteppingHelper.hpp +++ b/Core/include/Acts/Propagator/detail/SteppingHelper.hpp @@ -39,7 +39,7 @@ namespace Acts::detail { /// @param logger [in] A @c Logger instance template IntersectionStatus updateSingleSurfaceStatus( - const stepper_t& stepper, typename typename stepper_t::State& state, + const stepper_t& stepper, typename stepper_t::State& state, const Surface& surface, std::uint8_t index, Direction direction, const BoundaryTolerance& boundaryTolerance, double surfaceTolerance, ConstrainedStep::Type stype, const Logger& logger) { From 01a57152fd756ab8a259cf0e4b71c67c574a72cc Mon Sep 17 00:00:00 2001 From: Andreas Stefl Date: Sat, 8 Feb 2025 13:08:57 +0100 Subject: [PATCH 10/11] fix --- Core/include/Acts/Propagator/Propagator.ipp | 10 ---------- 1 file changed, 10 deletions(-) diff --git a/Core/include/Acts/Propagator/Propagator.ipp b/Core/include/Acts/Propagator/Propagator.ipp index 8a373480ac9..af6089b6e85 100644 --- a/Core/include/Acts/Propagator/Propagator.ipp +++ b/Core/include/Acts/Propagator/Propagator.ipp @@ -252,11 +252,6 @@ auto Acts::Propagator::makeState( StateType state{eOptions, m_stepper.makeState(eOptions.stepping), m_navigator.makeState(eOptions.navigation)}; - static_assert( - detail::propagator_stepper_compatible_with, - "Step method of the Stepper is not compatible with the propagator " - "state"); - initialize(state, start); return state; @@ -291,11 +286,6 @@ auto Acts::Propagator::makeState( StateType state{eOptions, m_stepper.makeState(eOptions.stepping), m_navigator.makeState(eOptions.navigation)}; - static_assert( - detail::propagator_stepper_compatible_with, - "Step method of the Stepper is not compatible with the propagator " - "state"); - initialize(state, start); return state; From 63409512b996263773d8180d68d8be6cdd06bac1 Mon Sep 17 00:00:00 2001 From: Andreas Stefl Date: Mon, 10 Feb 2025 11:13:51 +0100 Subject: [PATCH 11/11] doc --- Core/include/Acts/Propagator/AtlasStepper.hpp | 10 ++++++++-- Core/include/Acts/Propagator/EigenStepper.hpp | 6 ++++-- Core/include/Acts/Propagator/StraightLineStepper.hpp | 11 +++++------ Core/include/Acts/Propagator/SympyStepper.hpp | 4 +++- 4 files changed, 20 insertions(+), 11 deletions(-) diff --git a/Core/include/Acts/Propagator/AtlasStepper.hpp b/Core/include/Acts/Propagator/AtlasStepper.hpp index 61752c09414..16cafc17639 100644 --- a/Core/include/Acts/Propagator/AtlasStepper.hpp +++ b/Core/include/Acts/Propagator/AtlasStepper.hpp @@ -1147,10 +1147,16 @@ class AtlasStepper { /// Perform the actual step on the state /// - /// @param state is the provided stepper state (caller keeps thread locality) + /// @param [in,out] state State of the stepper /// @param propDir is the direction of propagation - /// @param material is the material properties + /// @param material is the optional volume material we are stepping through. + // This is simply ignored if `nullptr`. /// @return the result of the step + /// + /// @note The state contains the desired step size. It can be negative during + /// backwards track propagation, and since we're using an adaptive + /// algorithm, it can be modified by the stepper class during + /// propagation. Result step(State& state, Direction propDir, const IVolumeMaterial* material) const { (void)material; diff --git a/Core/include/Acts/Propagator/EigenStepper.hpp b/Core/include/Acts/Propagator/EigenStepper.hpp index c918fc91159..b17b2160a19 100644 --- a/Core/include/Acts/Propagator/EigenStepper.hpp +++ b/Core/include/Acts/Propagator/EigenStepper.hpp @@ -389,10 +389,12 @@ class EigenStepper { /// Perform a Runge-Kutta track parameter propagation step /// - /// @param [in,out] state The state of the stepper + /// @param [in,out] state State of the stepper /// @param propDir is the direction of propagation - /// @param material is the material properties + /// @param material is the optional volume material we are stepping through. + // This is simply ignored if `nullptr`. /// @return the result of the step + /// /// @note The state contains the desired step size. It can be negative during /// backwards track propagation, and since we're using an adaptive /// algorithm, it can be modified by the stepper class during diff --git a/Core/include/Acts/Propagator/StraightLineStepper.hpp b/Core/include/Acts/Propagator/StraightLineStepper.hpp index f32ba9a5116..b527ea3a04c 100644 --- a/Core/include/Acts/Propagator/StraightLineStepper.hpp +++ b/Core/include/Acts/Propagator/StraightLineStepper.hpp @@ -364,15 +364,14 @@ class StraightLineStepper { /// Perform a straight line propagation step /// - /// @param [in,out] state The state of the stepper + /// @param [in,out] state State of the stepper /// @param propDir is the direction of propagation - /// @param material is the material properties + /// @param material is the optional volume material we are stepping through. + // This is simply ignored if `nullptr`. /// @return the result of the step /// - /// The state contains the desired step size, - /// it can be negative during backwards track propagation, - /// and since we're using an adaptive algorithm, it can - /// be modified by the stepper class during propagation. + /// @note The state contains the desired step size. It can be negative during + /// backwards track propagation. Result step(State& state, Direction propDir, const IVolumeMaterial* material) const { (void)material; diff --git a/Core/include/Acts/Propagator/SympyStepper.hpp b/Core/include/Acts/Propagator/SympyStepper.hpp index fab189e0014..f64be7cb8f0 100644 --- a/Core/include/Acts/Propagator/SympyStepper.hpp +++ b/Core/include/Acts/Propagator/SympyStepper.hpp @@ -359,8 +359,10 @@ class SympyStepper { /// /// @param [in,out] state State of the stepper /// @param propDir is the direction of propagation - /// @param material is the material properties + /// @param material is the optional volume material we are stepping through. + // This is simply ignored if `nullptr`. /// @return the result of the step + /// /// @note The state contains the desired step size. It can be negative during /// backwards track propagation, and since we're using an adaptive /// algorithm, it can be modified by the stepper class during