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/Definitions/PdgParticle.hpp"
0012 #include "Acts/Material/IVolumeMaterial.hpp"
0013 #include "Acts/Material/Interactions.hpp"
0014 #include "Acts/Propagator/EigenStepperError.hpp"
0015 #include "Acts/Propagator/detail/SympyCovarianceEngine.hpp"
0016 #include "Acts/Propagator/detail/SympyJacobianEngine.hpp"
0017
0018 #include <cmath>
0019
0020 #include "codegen/sympy_stepper_math.hpp"
0021
0022 namespace Acts {
0023
0024 SympyStepper::SympyStepper(std::shared_ptr<const MagneticFieldProvider> bField)
0025     : m_bField(std::move(bField)) {}
0026
0027 SympyStepper::SympyStepper(const Config& config) : m_bField(config.bField) {}
0028
0029 SympyStepper::State SympyStepper::makeState(const Options& options) const {
0030   State state{options, m_bField->makeCache(options.magFieldContext)};
0031   return state;
0032 }
0033
0034 void SympyStepper::initialize(State& state,
0035                               const BoundTrackParameters& par) const {
0036   return initialize(state, par.parameters(), par.covariance(),
0037                     par.particleHypothesis(), par.referenceSurface());
0038 }
0039
0040 void SympyStepper::initialize(State& state, const BoundVector& boundParams,
0041                               const std::optional<BoundMatrix>& cov,
0042                               ParticleHypothesis particleHypothesis,
0043                               const Surface& surface) const {
0044   FreeVector freeParams = transformBoundToFreeParameters(
0045       surface, state.options.geoContext, boundParams);
0046
0047   state.particleHypothesis = particleHypothesis;
0048
0049   state.pathAccumulated = 0;
0050   state.nSteps = 0;
0051   state.nStepTrials = 0;
0052   state.stepSize = ConstrainedStep();
0053   state.stepSize.setAccuracy(state.options.initialStepSize);
0054   state.stepSize.setUser(state.options.maxStepSize);
0055   state.previousStepSize = 0;
0056   state.statistics = StepperStatistics();
0057
0058   state.pars = freeParams;
0059
0060   // Init the jacobian matrix if needed
0061   state.covTransport = cov.has_value();
0062   if (state.covTransport) {
0063     // set the covariance transport flag to true and copy
0064     state.cov = *cov;
0065     state.jacToGlobal = surface.boundToFreeJacobian(
0066         state.options.geoContext, freeParams.segment<3>(eFreePos0),
0067         freeParams.segment<3>(eFreeDir0));
0068     state.jacobian = BoundMatrix::Identity();
0069     state.jacTransport = FreeMatrix::Identity();
0070     state.derivative = FreeVector::Zero();
0071   }
0072 }
0073
0074 Result<std::tuple<BoundTrackParameters, BoundMatrix, double>>
0075 SympyStepper::boundState(
0076     State& state, const Surface& surface, bool transportCov,
0077     const FreeToBoundCorrection& freeToBoundCorrection) const {
0078   std::optional<FreeMatrix> additionalFreeCovariance =
0079       state.materialEffectsAccumulator.computeAdditionalFreeCovariance(
0080           direction(state));
0081   state.materialEffectsAccumulator.reset();
0082   return detail::sympy::boundState(
0083       state.options.geoContext, surface, state.cov, state.jacobian,
0084       state.jacTransport, state.derivative, state.jacToGlobal,
0085       additionalFreeCovariance, state.pars, state.particleHypothesis,
0086       state.covTransport && transportCov, state.pathAccumulated,
0087       freeToBoundCorrection);
0088 }
0089
0090 bool SympyStepper::prepareCurvilinearState(State& state) const {
0091   // TODO implement like in EigenStepper
0092   (void)state;
0093   return true;
0094 }
0095
0096 std::tuple<BoundTrackParameters, BoundMatrix, double>
0097 SympyStepper::curvilinearState(State& state, bool transportCov) const {
0098   std::optional<FreeMatrix> additionalFreeCovariance =
0099       state.materialEffectsAccumulator.computeAdditionalFreeCovariance(
0100           direction(state));
0101   state.materialEffectsAccumulator.reset();
0102   return detail::sympy::curvilinearState(
0103       state.cov, state.jacobian, state.jacTransport, state.derivative,
0104       state.jacToGlobal, additionalFreeCovariance, state.pars,
0105       state.particleHypothesis, state.covTransport && transportCov,
0106       state.pathAccumulated);
0107 }
0108
0109 void SympyStepper::update(State& state, const FreeVector& freeParams,
0110                           const BoundVector& /*boundParams*/,
0111                           const Covariance& covariance,
0112                           const Surface& surface) const {
0113   state.pars = freeParams;
0114   state.cov = covariance;
0115   state.jacToGlobal = surface.boundToFreeJacobian(
0116       state.options.geoContext, freeParams.template segment<3>(eFreePos0),
0117       freeParams.template segment<3>(eFreeDir0));
0118 }
0119
0120 void SympyStepper::update(State& state, const Vector3& uposition,
0121                           const Vector3& udirection, double qOverP,
0122                           double time) const {
0123   state.pars.template segment<3>(eFreePos0) = uposition;
0124   state.pars.template segment<3>(eFreeDir0) = udirection;
0125   state.pars[eFreeTime] = time;
0126   state.pars[eFreeQOverP] = qOverP;
0127 }
0128
0129 void SympyStepper::transportCovarianceToCurvilinear(State& state) const {
0130   detail::sympy::transportCovarianceToCurvilinear(
0131       state.cov, state.jacobian, state.jacTransport, state.derivative,
0132       state.jacToGlobal, std::nullopt,
0133       state.pars.template segment<3>(eFreeDir0));
0134 }
0135
0136 void SympyStepper::transportCovarianceToBound(
0137     State& state, const Surface& surface,
0138     const FreeToBoundCorrection& freeToBoundCorrection) const {
0139   detail::sympy::transportCovarianceToBound(
0140       state.options.geoContext, surface, state.cov, state.jacobian,
0141       state.jacTransport, state.derivative, state.jacToGlobal, std::nullopt,
0142       state.pars, freeToBoundCorrection);
0143 }
0144
0145 Result<double> SympyStepper::step(State& state, Direction propDir,
0146                                   const IVolumeMaterial* material) const {
0147   double h = state.stepSize.value() * propDir;
0148
0149   const double initialH = h;
0150   const Direction timeDirection = Direction::fromScalarZeroAsPositive(h);
0151
0152   const Vector3 pos = position(state);
0153   const Vector3 dir = direction(state);
0154   const double t = time(state);
0155   const double qop = qOverP(state);
0156   const double pabs = absoluteMomentum(state);
0157   const double m = particleHypothesis(state).mass();
0158   const PdgParticle absPdg = particleHypothesis(state).absolutePdg();
0159   const double q = charge(state);
0160   const double absQ = std::abs(q);
0161
0162   if (state.options.doDense && material != nullptr &&
0163       pabs < state.options.dense.momentumCutOff) {
0164     return EigenStepperError::StepInvalid;
0165   }
0166
0167   const auto getB = [&](const double* p) -> Result<Vector3> {
0168     return getField(state, {p[0], p[1], p[2]});
0169   };
0170
0171   const auto getG = [&](const double* p, double l) -> double {
0172     double newPabs = particleHypothesis(state).extractMomentum(l);
0173     if (newPabs < state.options.dense.momentumCutOff) {
0174       return 0.;
0175     }
0176
0177     if (state.options.dense.meanEnergyLoss) {
0178       return timeDirection *
0179              computeEnergyLossMean(
0180                  MaterialSlab(material->material({p[0], p[1], p[2]}),
0181                               1.0f * UnitConstants::mm),
0182                  absPdg, m, l, absQ);
0183     } else {
0184       return timeDirection *
0185              computeEnergyLossMode(
0186                  MaterialSlab(material->material({p[0], p[1], p[2]}),
0187                               1.0f * UnitConstants::mm),
0188                  absPdg, m, l, absQ);
0189     }
0190   };
0191
0192   const auto calcStepSizeScaling = [&](const double errorEstimate_) -> double {
0193     // For details about these values see ATL-SOFT-PUB-2009-001
0194     constexpr double lower = 0.25;
0195     constexpr double upper = 4.0;
0196     // This is given by the order of the Runge-Kutta method
0197     constexpr double exponent = 0.25;
0198
0199     double x = state.options.stepTolerance / errorEstimate_;
0200
0201     if constexpr (exponent == 0.25) {
0202       // This is 3x faster than std::pow
0203       x = std::sqrt(std::sqrt(x));
0204     } else {
0205       x = std::pow(x, exponent);
0206     }
0207
0208     return std::clamp(x, lower, upper);
0209   };
0210
0211   std::size_t nStepTrials = 0;
0212   double errorEstimate = 0.;
0213
0214   while (true) {
0215     ++nStepTrials;
0216     ++state.statistics.nAttemptedSteps;
0217
0218     // For details about the factor 4 see ATL-SOFT-PUB-2009-001
0219     Result<bool> res = Result<bool>::success(false);
0220     if (!state.options.doDense || material == nullptr) {
0221       res =
0222           rk4_vacuum(pos.data(), dir.data(), t, h, qop, m, pabs, getB,
0223                      &errorEstimate, 4 * state.options.stepTolerance,
0224                      state.pars.template segment<3>(eFreePos0).data(),
0225                      state.pars.template segment<1>(eFreeTime).data(),
0226                      state.pars.template segment<3>(eFreeDir0).data(),
0227                      state.derivative.data(),
0228                      state.covTransport ? state.jacTransport.data() : nullptr);
0229     } else {
0230       res = rk4_dense(pos.data(), dir.data(), t, h, qop, m, q, pabs, getB, getG,
0231                       &errorEstimate, 4 * state.options.stepTolerance,
0232                       state.pars.template segment<3>(eFreePos0).data(),
0233                       state.pars.template segment<1>(eFreeTime).data(),
0234                       state.pars.template segment<3>(eFreeDir0).data(),
0235                       state.pars.template segment<1>(eFreeQOverP).data(),
0236                       state.derivative.data(),
0237                       state.covTransport ? state.jacTransport.data() : nullptr);
0238     }
0239     if (!res.ok()) {
0240       return res.error();
0241     }
0242     // Protect against division by zero
0243     errorEstimate = std::max(1e-20, errorEstimate);
0244
0245     if (*res) {
0246       break;
0247     }
0248
0249     ++state.statistics.nRejectedSteps;
0250
0251     const double stepSizeScaling = calcStepSizeScaling(errorEstimate);
0252     h *= stepSizeScaling;
0253
0254     // If step size becomes too small the particle remains at the initial
0255     // place
0256     if (std::abs(h) < std::abs(state.options.stepSizeCutOff)) {
0257       // Not moving due to too low momentum needs an aborter
0258       return EigenStepperError::StepSizeStalled;
0259     }
0260
0261     // If the parameter is off track too much or given stepSize is not
0262     // appropriate
0263     if (nStepTrials > state.options.maxRungeKuttaStepTrials) {
0264       // Too many trials, have to abort
0265       return EigenStepperError::StepSizeAdjustmentFailed;
0266     }
0267   }
0268
0269   state.pathAccumulated += h;
0270   ++state.nSteps;
0271   state.nStepTrials += nStepTrials;
0272
0273   ++state.statistics.nSuccessfulSteps;
0274   if (propDir != Direction::fromScalarZeroAsPositive(initialH)) {
0275     ++state.statistics.nReverseSteps;
0276   }
0277   state.statistics.pathLength += h;
0278   state.statistics.absolutePathLength += std::abs(h);
0279
0280   const double stepSizeScaling = calcStepSizeScaling(errorEstimate);
0281   const double nextAccuracy = std::abs(h * stepSizeScaling);
0282   const double previousAccuracy = std::abs(state.stepSize.accuracy());
0283   const double initialStepLength = std::abs(initialH);
0284   if (nextAccuracy < initialStepLength || nextAccuracy > previousAccuracy) {
0285     state.stepSize.setAccuracy(nextAccuracy);
0286   }
0287
0288   if (state.options.doDense &&
0289       (material != nullptr || !state.materialEffectsAccumulator.isVacuum())) {
0290     if (state.materialEffectsAccumulator.isVacuum()) {
0291       state.materialEffectsAccumulator.initialize(
0292           state.options.maxXOverX0Step, particleHypothesis(state), pabs);
0293     }
0294
0295     Material mat =
0296         material != nullptr ? material->material(pos) : Material::Vacuum();
0297
0298     state.materialEffectsAccumulator.accumulate(mat, propDir * h, qop,
0299                                                 qOverP(state));
0300   }
0301
0302   return h;
0303 }
0304
0305 void SympyStepper::setIdentityJacobian(State& state) const {
0306   state.jacobian = BoundMatrix::Identity();
0307 }
0308
0309 }  // namespace Acts