src/Propagator/HelixStepper.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/HelixStepper.hpp"
0010
0011 #include "Acts/Definitions/TrackParametrization.hpp"
0012 #include "Acts/EventData/TransformationHelpers.hpp"
0013 #include "Acts/Propagator/EigenStepperError.hpp"
0014 #include "Acts/Propagator/detail/CovarianceEngine.hpp"
0015 #include "Acts/Surfaces/BoundaryTolerance.hpp"
0016 #include "Acts/Surfaces/SurfaceError.hpp"
0017 #include "Acts/Utilities/MathHelpers.hpp"
0018
0019 #include <algorithm>
0020 #include <cmath>
0021
0022 namespace Acts {
0023
0024 namespace {
0025
0026 /// Even functions of the turning angle theta = |omega| * |h| that the helix
0027 /// and its jacobian are built from. Each one is normalised so that it is
0028 /// finite at theta = 0.
0029 struct HelixCoefficients {
0030   /// sin(theta) / theta
0031   double f1 = 0;
0032   /// (1 - cos(theta)) / theta^2
0033   double f2 = 0;
0034   /// (theta - sin(theta)) / theta^3
0035   double g2 = 0;
0036   /// (sin(theta) - theta * cos(theta)) / theta^3
0037   double h1 = 0;
0038   /// (theta^2 / 2 + 1 - cos(theta) - theta * sin(theta)) / theta^4
0039   double h2 = 0;
0040 };
0041
0042 HelixCoefficients helixCoefficients(double theta2) {
0043   HelixCoefficients c;
0044   // The closed forms cancel for a small angle, so use the series there. The
0045   // first omitted term is below 1e-17 at the threshold.
0046   if (theta2 < 1e-3) {
0047     const double t2 = theta2;
0048     c.f1 = 1. + t2 * (-1. / 6. + t2 * (1. / 120. + t2 * (-1. / 5040.)));
0049     c.f2 = 1. / 2. + t2 * (-1. / 24. + t2 * (1. / 720. + t2 * (-1. / 40320.)));
0050     c.g2 =
0051         1. / 6. + t2 * (-1. / 120. + t2 * (1. / 5040. + t2 * (-1. / 362880.)));
0052     c.h1 = 1. / 3. + t2 * (-1. / 30. + t2 * (1. / 840. + t2 * (-1. / 45360.)));
0053     c.h2 =
0054         1. / 8. + t2 * (-1. / 144. + t2 * (1. / 5760. + t2 * (-1. / 403200.)));
0055     return c;
0056   }
0057   const double theta = std::sqrt(theta2);
0058   // Half-angle forms: two trigonometric calls, and 1 - cos(theta) without
0059   // the cancellation
0060   const double halfSin = std::sin(0.5 * theta);
0061   const double halfCos = std::cos(0.5 * theta);
0062   const double sinTheta = 2. * halfSin * halfCos;
0063   const double oneMinusCos = 2. * halfSin * halfSin;
0064   const double cosTheta = 1. - oneMinusCos;
0065   const double theta3 = theta2 * theta;
0066   c.f1 = sinTheta / theta;
0067   c.f2 = oneMinusCos / theta2;
0068   c.g2 = (theta - sinTheta) / theta3;
0069   c.h1 = (sinTheta - theta * cosTheta) / theta3;
0070   c.h2 = (0.5 * theta2 + oneMinusCos - theta * sinTheta) / (theta2 * theta2);
0071   return c;
0072 }
0073
0074 /// Cross product matrix W with W * v = omega x v
0075 SquareMatrix3 crossMatrix(const Vector3& omega) {
0076   SquareMatrix3 w;
0077   // clang-format off
0078   w <<         0., -omega.z(),  omega.y(),
0079         omega.z(),         0., -omega.x(),
0080        -omega.y(),  omega.x(),         0.;
0081   // clang-format on
0082   return w;
0083 }
0084
0085 }  // namespace
0086
0087 HelixStepper::HelixStepper(std::shared_ptr<const MagneticFieldProvider> bField)
0088     : m_bField(std::move(bField)) {}
0089
0090 HelixStepper::HelixStepper(const Config& config) : m_bField(config.bField) {}
0091
0092 HelixStepper::State HelixStepper::makeState(const Options& options) const {
0093   State state{options, m_bField->makeCache(options.magFieldContext)};
0094   return state;
0095 }
0096
0097 void HelixStepper::initialize(State& state, const BoundParameters& par) const {
0098   initialize(state, par.parameters(), par.covariance(),
0099              par.particleHypothesis(), par.referenceSurface());
0100 }
0101
0102 void HelixStepper::initialize(State& state, const BoundVector& boundParams,
0103                               const std::optional<BoundMatrix>& cov,
0104                               ParticleHypothesis particleHypothesis,
0105                               const Surface& surface) const {
0106   FreeVector freeParams = transformBoundToFreeParameters(
0107       surface, state.options.geoContext, boundParams);
0108
0109   state.particleHypothesis = particleHypothesis;
0110
0111   state.pathAccumulated = 0;
0112   state.nSteps = 0;
0113   state.stepSize = ConstrainedStep();
0114   state.stepSize.setAccuracy(state.options.initialStepSize);
0115   state.stepSize.setUser(state.options.maxStepSize);
0116   state.previousStepSize = 0;
0117   state.statistics = StepperStatistics();
0118
0119   state.pars = freeParams;
0120   state.field.reset();
0121
0122   state.cov = cov;
0123   if (state.cov.has_value()) {
0124     state.jacToGlobal = surface.boundToFreeJacobian(
0125         state.options.geoContext, freeParams.segment<3>(eFreePos0),
0126         freeParams.segment<3>(eFreeDir0));
0127     state.jacTransport = FreeMatrix::Identity();
0128     state.derivative = FreeVector::Zero();
0129   }
0130 }
0131
0132 Result<HelixStepper::BoundParameters> HelixStepper::boundParameters(
0133     const State& state, const Surface& surface) const {
0134   return detail::boundParameters(state.options.geoContext, surface, state.pars,
0135                                  state.cov, state.particleHypothesis);
0136 }
0137
0138 bool HelixStepper::prepareCurvilinearState(State& state) const {
0139   // The derivatives are only missing if no step was executed yet.
0140   if (state.pathAccumulated != 0) {
0141     return true;
0142   }
0143
0144   if (!state.field.has_value()) {
0145     auto fieldRes = getField(state, position(state));
0146     if (!fieldRes.ok()) {
0147       return false;
0148     }
0149     state.field = *fieldRes;
0150   }
0151
0152   const double qop = qOverP(state);
0153   const double mass = state.particleHypothesis.mass();
0154   state.derivative.segment<3>(eFreePos0) = direction(state);
0155   state.derivative[eFreeTime] = fastHypot(1., mass / absoluteMomentum(state));
0156   state.derivative.segment<3>(eFreeDir0) =
0157       qop * direction(state).cross(*state.field);
0158   state.derivative[eFreeQOverP] = 0.;
0159   return true;
0160 }
0161
0162 HelixStepper::BoundParameters HelixStepper::curvilinearParameters(
0163     const State& state) const {
0164   return detail::curvilinearParameters(state.pars, state.cov,
0165                                        state.particleHypothesis);
0166 }
0167
0168 void HelixStepper::update(State& state, const FreeVector& freeParams,
0169                           const BoundVector& /*boundParams*/,
0170                           const Covariance& covariance,
0171                           const Surface& surface) const {
0172   state.pars = freeParams;
0173   state.field.reset();
0174   if (state.cov.has_value()) {
0175     state.cov = covariance;
0176   }
0177   state.jacToGlobal = surface.boundToFreeJacobian(
0178       state.options.geoContext, freeParams.template segment<3>(eFreePos0),
0179       freeParams.template segment<3>(eFreeDir0));
0180   state.jacTransport = FreeMatrix::Identity();
0181   state.derivative = FreeVector::Zero();
0182 }
0183
0184 void HelixStepper::update(State& state, const Vector3& uposition,
0185                           const Vector3& udirection, double qop,
0186                           double time) const {
0187   // Material interactions keep the position, so the cached field stays valid.
0188   if (uposition != position(state)) {
0189     state.field.reset();
0190   }
0191   state.pars.template segment<3>(eFreePos0) = uposition;
0192   state.pars.template segment<3>(eFreeDir0) = udirection;
0193   state.pars[eFreeTime] = time;
0194   state.pars[eFreeQOverP] = qop;
0195 }
0196
0197 HelixStepper::Jacobian HelixStepper::transportToCurvilinear(
0198     State& state) const {
0199   Jacobian jacobian = Jacobian::Identity();
0200   if (!state.cov.has_value()) {
0201     return jacobian;
0202   }
0203   detail::transportCovarianceToCurvilinear(
0204       *state.cov, jacobian, state.jacTransport, state.derivative,
0205       state.jacToGlobal, std::nullopt, direction(state));
0206   return jacobian;
0207 }
0208
0209 Result<HelixStepper::Jacobian> HelixStepper::transportToBound(
0210     State& state, const Surface& surface,
0211     const FreeToBoundCorrection& freeToBoundCorrection) const {
0212   if (!surface.isOnSurface(state.options.geoContext, position(state),
0213                            direction(state), BoundaryTolerance::Infinite())) {
0214     return Result<Jacobian>::failure(SurfaceError::GlobalPositionNotOnSurface);
0215   }
0216
0217   Jacobian jacobian = Jacobian::Identity();
0218   if (!state.cov.has_value()) {
0219     return Result<Jacobian>::success(jacobian);
0220   }
0221   detail::transportCovarianceToBound(
0222       state.options.geoContext, surface, *state.cov, jacobian,
0223       state.jacTransport, state.derivative, state.jacToGlobal, std::nullopt,
0224       state.pars, freeToBoundCorrection);
0225   return Result<Jacobian>::success(jacobian);
0226 }
0227
0228 Result<double> HelixStepper::step(State& state, Direction propDir,
0229                                   const IVolumeMaterial* /*material*/) const {
0230   const Vector3 pos = position(state);
0231   const Vector3 dir = direction(state);
0232   const double qop = qOverP(state);
0233
0234   if (!state.field.has_value()) {
0235     auto fieldRes = getField(state, pos);
0236     if (!fieldRes.ok()) {
0237       return fieldRes.error();
0238     }
0239     state.field = *fieldRes;
0240   }
0241   const Vector3 startField = *state.field;
0242
0243   // The equation of motion is dT/ds = T x omega. Its solution is a rotation
0244   // of T about omega:
0245   //   T(h) = T - h f1 a + h^2 f2 b
0246   //   r(h) = r + h T - h^2 f2 a + h^3 g2 b
0247   // with a = omega x T and b = omega x a.
0248   const Vector3 omega = qop * startField;
0249   const Vector3 a = omega.cross(dir);
0250   const Vector3 b = omega.cross(a);
0251   const double omega2 = omega.squaredNorm();
0252
0253   const auto calcStepSizeScaling = [&](const double errorEstimate) -> double {
0254     constexpr double lower = 0.25;
0255     constexpr double upper = 4.0;
0256     // The error of a linear field change along the step grows with h^3
0257     const double x = std::cbrt(state.options.stepTolerance / errorEstimate);
0258     return std::clamp(x, lower, upper);
0259   };
0260
0261   const double initialH = state.stepSize.value() * propDir;
0262   double h = initialH;
0263   HelixCoefficients coeffs;
0264   Vector3 endPos;
0265   Vector3 endDir;
0266   Vector3 endField;
0267   double errorEstimate = 0.;
0268   std::size_t nStepTrials = 0;
0269
0270   while (true) {
0271     ++nStepTrials;
0272     ++state.statistics.nAttemptedSteps;
0273
0274     coeffs = helixCoefficients(omega2 * h * h);
0275     const double h2 = h * h;
0276     endDir = dir - h * coeffs.f1 * a + h2 * coeffs.f2 * b;
0277     endPos = pos + h * dir - h2 * coeffs.f2 * a + h2 * h * coeffs.g2 * b;
0278
0279     // The end field is the start field of the next step, so an accepted step
0280     // costs no extra lookup.
0281     auto fieldRes = getField(state, endPos);
0282     if (!fieldRes.ok()) {
0283       return fieldRes.error();
0284     }
0285     endField = *fieldRes;
0286
0287     // Position deviation from a field that changes linearly along the step
0288     errorEstimate = std::max(
0289         1e-20, h2 / 6. * (qop * dir.cross(endField - startField)).norm());
0290
0291     if (errorEstimate <= state.options.stepTolerance) {
0292       break;
0293     }
0294
0295     ++state.statistics.nRejectedSteps;
0296
0297     h *= calcStepSizeScaling(errorEstimate);
0298
0299     if (std::abs(h) < std::abs(state.options.stepSizeCutOff)) {
0300       return EigenStepperError::StepSizeStalled;
0301     }
0302
0303     if (nStepTrials > state.options.maxRungeKuttaStepTrials) {
0304       return EigenStepperError::StepSizeAdjustmentFailed;
0305     }
0306   }
0307
0308   const double mass = state.particleHypothesis.mass();
0309   const double absMom = absoluteMomentum(state);
0310   // time propagates along distance as 1/b = sqrt(1 + m²/p²)
0311   const double dtds = fastHypot(1., mass / absMom);
0312
0313   if (state.cov.has_value()) {
0314     const double h2 = h * h;
0315     const SquareMatrix3 w = crossMatrix(omega);
0316     const SquareMatrix3 w2 = w * w;
0317
0318     // The helix does not depend on the start position and time, so the step
0319     // jacobian D has D11 = I and D21 = 0. Only the right half is non-trivial.
0320     Matrix<4, 4> dTopRight = Matrix<4, 4>::Zero();
0321     Matrix<4, 4> dBottomRight = Matrix<4, 4>::Identity();
0322
0323     // dr/dT and dT/dT
0324     dTopRight.topLeftCorner<3, 3>() = h * SquareMatrix3::Identity() -
0325                                       h2 * coeffs.f2 * w +
0326                                       h2 * h * coeffs.g2 * w2;
0327     dBottomRight.topLeftCorner<3, 3>() =
0328         SquareMatrix3::Identity() - h * coeffs.f1 * w + h2 * coeffs.f2 * w2;
0329
0330     // omega = qop * B, so d/d(qop) = B x d/d(omega) applied to the
0331     // integrals of T:
0332     //   dT(h)/d(qop) = h T(h) x B
0333     //   dr(h)/d(qop) = (h^2/2 T - h^3 h1 a + h^4 h2 b) x B
0334     dTopRight.block<3, 1>(0, 3) =
0335         (0.5 * h2 * dir - h2 * h * coeffs.h1 * a + h2 * h2 * coeffs.h2 * b)
0336             .cross(startField);
0337     dBottomRight.block<3, 1>(0, 3) = h * endDir.cross(startField);
0338
0339     // dt/d(qop) = h m^2 (q/p) / (q^2 dt/ds). A neutral particle has a fixed
0340     // momentum hypothesis, so its time does not depend on qop.
0341     const double q = charge(state);
0342     dTopRight(3, 3) = q != 0. ? h * mass * mass * qop / (q * q * dtds) : 0.;
0343
0344     // See EigenStepper::step for the blocked multiplication.
0345     state.jacTransport.template topRightCorner<4, 4>() +=
0346         dTopRight * state.jacTransport.template bottomRightCorner<4, 4>();
0347     state.jacTransport.template bottomRightCorner<4, 4>() =
0348         (dBottomRight * state.jacTransport.template bottomRightCorner<4, 4>())
0349             .eval();
0350
0351     state.derivative.template segment<3>(eFreePos0) = endDir;
0352     state.derivative[eFreeTime] = dtds;
0353     state.derivative.template segment<3>(eFreeDir0) = endDir.cross(omega);
0354     state.derivative[eFreeQOverP] = 0.;
0355   }
0356
0357   state.pars.template segment<3>(eFreePos0) = endPos;
0358   state.pars.template segment<3>(eFreeDir0) = endDir;
0359   state.pars[eFreeTime] += h * dtds;
0360   state.field = endField;
0361
0362   state.pathAccumulated += h;
0363   ++state.nSteps;
0364
0365   ++state.statistics.nSuccessfulSteps;
0366   if (propDir != Direction::fromScalarZeroAsPositive(initialH)) {
0367     ++state.statistics.nReverseSteps;
0368   }
0369   state.statistics.pathLength += h;
0370   state.statistics.absolutePathLength += std::abs(h);
0371
0372   const double nextAccuracy = std::abs(h * calcStepSizeScaling(errorEstimate));
0373   const double previousAccuracy = std::abs(state.stepSize.accuracy());
0374   const double initialStepLength = std::abs(initialH);
0375   if (nextAccuracy < initialStepLength || nextAccuracy > previousAccuracy) {
0376     state.stepSize.setAccuracy(nextAccuracy);
0377   }
0378
0379   return h;
0380 }
0381
0382 }  // namespace Acts