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