src/Propagator/SympyStepper.cpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 #include "Acts/Propagator/SympyStepper.hpp"
0010
0011 #include "Acts/Propagator/EigenStepperError.hpp"
0012 #include "Acts/Propagator/detail/SympyCovarianceEngine.hpp"
0013 #include "Acts/Propagator/detail/SympyJacobianEngine.hpp"
0014
0015 #include <cmath>
0016
0017 #include "codegen/sympy_stepper_math.hpp"
0018
0019 namespace Acts {
0020
0021 SympyStepper::SympyStepper(std::shared_ptr<const MagneticFieldProvider> bField)
0022     : m_bField(std::move(bField)) {}
0023
0024 SympyStepper::SympyStepper(const Config& config) : m_bField(config.bField) {}
0025
0026 SympyStepper::State SympyStepper::makeState(const Options& options) const {
0027   State state{options, m_bField->makeCache(options.magFieldContext)};
0028   return state;
0029 }
0030
0031 void SympyStepper::initialize(State& state,
0032                               const BoundTrackParameters& par) const {
0033   return initialize(state, par.parameters(), par.covariance(),
0034                     par.particleHypothesis(), par.referenceSurface());
0035 }
0036
0037 void SympyStepper::initialize(State& state, const BoundVector& boundParams,
0038                               const std::optional<BoundMatrix>& cov,
0039                               ParticleHypothesis particleHypothesis,
0040                               const Surface& surface) const {
0041   FreeVector freeParams = transformBoundToFreeParameters(
0042       surface, state.options.geoContext, boundParams);
0043
0044   state.particleHypothesis = particleHypothesis;
0045
0046   state.pathAccumulated = 0;
0047   state.nSteps = 0;
0048   state.nStepTrials = 0;
0049   state.stepSize = ConstrainedStep();
0050   state.stepSize.setAccuracy(state.options.initialStepSize);
0051   state.stepSize.setUser(state.options.maxStepSize);
0052   state.previousStepSize = 0;
0053   state.statistics = StepperStatistics();
0054
0055   state.pars = freeParams;
0056
0057   // Init the jacobian matrix if needed
0058   state.covTransport = cov.has_value();
0059   if (state.covTransport) {
0060     // set the covariance transport flag to true and copy
0061     state.cov = *cov;
0062     state.jacToGlobal = surface.boundToFreeJacobian(
0063         state.options.geoContext, freeParams.segment<3>(eFreePos0),
0064         freeParams.segment<3>(eFreeDir0));
0065     state.jacobian = BoundMatrix::Identity();
0066     state.jacTransport = FreeMatrix::Identity();
0067     state.derivative = FreeVector::Zero();
0068   }
0069 }
0070
0071 Result<std::tuple<BoundTrackParameters, BoundMatrix, double>>
0072 SympyStepper::boundState(
0073     State& state, const Surface& surface, bool transportCov,
0074     const FreeToBoundCorrection& freeToBoundCorrection) const {
0075   return detail::sympy::boundState(
0076       state.options.geoContext, surface, state.cov, state.jacobian,
0077       state.jacTransport, state.derivative, state.jacToGlobal, state.pars,
0078       state.particleHypothesis, state.covTransport && transportCov,
0079       state.pathAccumulated, freeToBoundCorrection);
0080 }
0081
0082 bool SympyStepper::prepareCurvilinearState(State& state) const {
0083   // TODO implement like in EigenStepper
0084   (void)state;
0085   return true;
0086 }
0087
0088 std::tuple<BoundTrackParameters, BoundMatrix, double>
0089 SympyStepper::curvilinearState(State& state, bool transportCov) const {
0090   return detail::sympy::curvilinearState(
0091       state.cov, state.jacobian, state.jacTransport, state.derivative,
0092       state.jacToGlobal, state.pars, state.particleHypothesis,
0093       state.covTransport && transportCov, state.pathAccumulated);
0094 }
0095
0096 void SympyStepper::update(State& state, const FreeVector& freeParams,
0097                           const BoundVector& /*boundParams*/,
0098                           const Covariance& covariance,
0099                           const Surface& surface) const {
0100   state.pars = freeParams;
0101   state.cov = covariance;
0102   state.jacToGlobal = surface.boundToFreeJacobian(
0103       state.options.geoContext, freeParams.template segment<3>(eFreePos0),
0104       freeParams.template segment<3>(eFreeDir0));
0105 }
0106
0107 void SympyStepper::update(State& state, const Vector3& uposition,
0108                           const Vector3& udirection, double qOverP,
0109                           double time) const {
0110   state.pars.template segment<3>(eFreePos0) = uposition;
0111   state.pars.template segment<3>(eFreeDir0) = udirection;
0112   state.pars[eFreeTime] = time;
0113   state.pars[eFreeQOverP] = qOverP;
0114 }
0115
0116 void SympyStepper::transportCovarianceToCurvilinear(State& state) const {
0117   detail::sympy::transportCovarianceToCurvilinear(
0118       state.cov, state.jacobian, state.jacTransport, state.derivative,
0119       state.jacToGlobal, state.pars.template segment<3>(eFreeDir0));
0120 }
0121
0122 void SympyStepper::transportCovarianceToBound(
0123     State& state, const Surface& surface,
0124     const FreeToBoundCorrection& freeToBoundCorrection) const {
0125   detail::sympy::transportCovarianceToBound(
0126       state.options.geoContext, surface, state.cov, state.jacobian,
0127       state.jacTransport, state.derivative, state.jacToGlobal, state.pars,
0128       freeToBoundCorrection);
0129 }
0130
0131 Result<double> SympyStepper::step(State& state, Direction propDir,
0132                                   const IVolumeMaterial* material) const {
0133   (void)material;
0134
0135   auto pos = position(state);
0136   auto dir = direction(state);
0137   double t = time(state);
0138   double qop = qOverP(state);
0139   double m = particleHypothesis(state).mass();
0140   double p_abs = absoluteMomentum(state);
0141
0142   auto getB = [&](const double* p) -> Result<Vector3> {
0143     return getField(state, {p[0], p[1], p[2]});
0144   };
0145
0146   const auto calcStepSizeScaling = [&](const double errorEstimate_) -> double {
0147     // For details about these values see ATL-SOFT-PUB-2009-001
0148     constexpr double lower = 0.25;
0149     constexpr double upper = 4.0;
0150     // This is given by the order of the Runge-Kutta method
0151     constexpr double exponent = 0.25;
0152
0153     double x = state.options.stepTolerance / errorEstimate_;
0154
0155     if constexpr (exponent == 0.25) {
0156       // This is 3x faster than std::pow
0157       x = std::sqrt(std::sqrt(x));
0158     } else {
0159       x = std::pow(x, exponent);
0160     }
0161
0162     return std::clamp(x, lower, upper);
0163   };
0164
0165   double h = state.stepSize.value() * propDir;
0166   double initialH = h;
0167   std::size_t nStepTrials = 0;
0168   double errorEstimate = 0.;
0169
0170   while (true) {
0171     ++nStepTrials;
0172     ++state.statistics.nAttemptedSteps;
0173
0174     // For details about the factor 4 see ATL-SOFT-PUB-2009-001
0175     Result<bool> res =
0176         rk4(pos.data(), dir.data(), t, h, qop, m, p_abs, getB, &errorEstimate,
0177             4 * state.options.stepTolerance,
0178             state.pars.template segment<3>(eFreePos0).data(),
0179             state.pars.template segment<3>(eFreeDir0).data(),
0180             state.pars.template segment<1>(eFreeTime).data(),
0181             state.derivative.data(),
0182             state.covTransport ? state.jacTransport.data() : nullptr);
0183     if (!res.ok()) {
0184       return res.error();
0185     }
0186     // Protect against division by zero
0187     errorEstimate = std::max(1e-20, errorEstimate);
0188
0189     if (*res) {
0190       break;
0191     }
0192
0193     ++state.statistics.nRejectedSteps;
0194
0195     const double stepSizeScaling = calcStepSizeScaling(errorEstimate);
0196     h *= stepSizeScaling;
0197
0198     // If step size becomes too small the particle remains at the initial
0199     // place
0200     if (std::abs(h) < std::abs(state.options.stepSizeCutOff)) {
0201       // Not moving due to too low momentum needs an aborter
0202       return EigenStepperError::StepSizeStalled;
0203     }
0204
0205     // If the parameter is off track too much or given stepSize is not
0206     // appropriate
0207     if (nStepTrials > state.options.maxRungeKuttaStepTrials) {
0208       // Too many trials, have to abort
0209       return EigenStepperError::StepSizeAdjustmentFailed;
0210     }
0211   }
0212
0213   state.pathAccumulated += h;
0214   ++state.nSteps;
0215   state.nStepTrials += nStepTrials;
0216
0217   ++state.statistics.nSuccessfulSteps;
0218   if (propDir != Direction::fromScalarZeroAsPositive(initialH)) {
0219     ++state.statistics.nReverseSteps;
0220   }
0221   state.statistics.pathLength += h;
0222   state.statistics.absolutePathLength += std::abs(h);
0223
0224   const double stepSizeScaling = calcStepSizeScaling(errorEstimate);
0225   const double nextAccuracy = std::abs(h * stepSizeScaling);
0226   const double previousAccuracy = std::abs(state.stepSize.accuracy());
0227   const double initialStepLength = std::abs(initialH);
0228   if (nextAccuracy < initialStepLength || nextAccuracy > previousAccuracy) {
0229     state.stepSize.setAccuracy(nextAccuracy);
0230   }
0231
0232   return h;
0233 }
0234
0235 void SympyStepper::setIdentityJacobian(State& state) const {
0236   state.jacobian = BoundMatrix::Identity();
0237 }
0238
0239 }  // namespace Acts