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/Material/IVolumeMaterial.hpp"
0012 #include "Acts/Propagator/EigenStepperError.hpp"
0013 #include "Acts/Propagator/detail/CovarianceEngine.hpp"
0014 #include "Acts/Propagator/detail/SympyBoundToFreeScaling.hpp"
0015 #include "Acts/Propagator/detail/SympyCovarianceEngine.hpp"
0016 #include "Acts/Propagator/detail/SympyJacobianEngine.hpp"
0017 #include "Acts/Surfaces/BoundaryTolerance.hpp"
0018 #include "Acts/Surfaces/SurfaceError.hpp"
0019
0020 #include <cmath>
0021 #include <span>
0022
0023 #include "detail/SympyStepperStep.hpp"
0024
0025 namespace Acts {
0026
0027 SympyStepper::SympyStepper(std::shared_ptr<const MagneticFieldProvider> bField)
0028     : m_bField(std::move(bField)) {}
0029
0030 SympyStepper::SympyStepper(const Config& config) : m_bField(config.bField) {}
0031
0032 SympyStepper::State SympyStepper::makeState(const Options& options) const {
0033   State state{options, m_bField->makeCache(options.magFieldContext)};
0034   return state;
0035 }
0036
0037 void SympyStepper::initialize(State& state, const BoundParameters& par) const {
0038   return initialize(state, par.parameters(), par.covariance(),
0039                     par.particleHypothesis(), par.referenceSurface());
0040 }
0041
0042 void SympyStepper::initialize(State& state, const BoundVector& boundParams,
0043                               const std::optional<BoundMatrix>& cov,
0044                               ParticleHypothesis particleHypothesis,
0045                               const Surface& surface) const {
0046   FreeVector freeParams = transformBoundToFreeParameters(
0047       surface, state.options.geoContext, boundParams);
0048
0049   state.particleHypothesis = particleHypothesis;
0050
0051   state.pathAccumulated = 0;
0052   state.nSteps = 0;
0053   state.stepSize = ConstrainedStep();
0054   state.stepSize.setAccuracy(state.options.initialStepSize);
0055   state.stepSize.setUser(state.options.maxStepSize);
0056   state.previousStepSize = 0;
0057   state.statistics = StepperStatistics();
0058
0059   state.pars = freeParams;
0060   state.field.reset();
0061   state.dtds = detail::sympyDtds(state);
0062
0063   // Init the jacobian matrix if needed
0064   state.cov = cov;
0065   if (state.cov.has_value()) {
0066     state.jacToGlobal = surface.boundToFreeJacobian(
0067         state.options.geoContext, freeParams.segment<3>(eFreePos0),
0068         freeParams.segment<3>(eFreeDir0));
0069     detail::sympy::toScaledBoundToFree(state.jacToGlobal,
0070                                        freeParams[eFreeQOverP]);
0071     state.derivative = FreeVector::Zero();
0072   }
0073 }
0074
0075 Result<SympyStepper::BoundParameters> SympyStepper::boundParameters(
0076     const State& state, const Surface& surface) const {
0077   return detail::boundParameters(state.options.geoContext, surface, state.pars,
0078                                  state.cov, state.particleHypothesis);
0079 }
0080
0081 bool SympyStepper::prepareCurvilinearState(State& state) const {
0082   // TODO implement like in EigenStepper
0083   static_cast<void>(state);
0084   return true;
0085 }
0086
0087 SympyStepper::BoundParameters SympyStepper::curvilinearParameters(
0088     const State& state) const {
0089   return detail::curvilinearParameters(state.pars, state.cov,
0090                                        state.particleHypothesis);
0091 }
0092
0093 void SympyStepper::update(State& state, const FreeVector& freeParams,
0094                           const BoundVector& /*boundParams*/,
0095                           const Covariance& covariance,
0096                           const Surface& surface) const {
0097   state.pars = freeParams;
0098   state.field.reset();
0099   state.dtds = detail::sympyDtds(state);
0100   if (state.cov.has_value()) {
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     detail::sympy::toScaledBoundToFree(state.jacToGlobal,
0106                                        freeParams[eFreeQOverP]);
0107     state.derivative = FreeVector::Zero();
0108   }
0109   state.materialEffectsAccumulator.reset();
0110 }
0111
0112 void SympyStepper::update(State& state, const Vector3& uposition,
0113                           const Vector3& udirection, double qOverP,
0114                           double time) const {
0115   if (state.cov.has_value()) {
0116     detail::sympy::rescaleBoundToFree(state.jacToGlobal,
0117                                       state.pars[eFreeQOverP], qOverP);
0118   }
0119   state.pars.template segment<3>(eFreePos0) = uposition;
0120   state.pars.template segment<3>(eFreeDir0) = udirection;
0121   state.pars[eFreeTime] = time;
0122   state.pars[eFreeQOverP] = qOverP;
0123   state.dtds = detail::sympyDtds(state);
0124   state.field.reset();
0125 }
0126
0127 SympyStepper::Jacobian SympyStepper::transportToCurvilinear(
0128     State& state) const {
0129   Jacobian jacobian = Jacobian::Identity();
0130   if (!state.cov.has_value()) {
0131     state.materialEffectsAccumulator.reset();
0132     return jacobian;
0133   }
0134   const std::optional<FreeMatrix> additionalFreeCovariance =
0135       state.materialEffectsAccumulator.computeAdditionalFreeCovariance(
0136           direction(state));
0137   state.materialEffectsAccumulator.reset();
0138   detail::sympy::fromScaledBoundToFree(state.jacToGlobal, qOverP(state));
0139   detail::sympy::transportCovarianceToCurvilinear(
0140       *state.cov, jacobian, state.derivative, state.jacToGlobal,
0141       additionalFreeCovariance, state.pars.template segment<3>(eFreeDir0));
0142   detail::sympy::toScaledBoundToFree(state.jacToGlobal, qOverP(state));
0143   return jacobian;
0144 }
0145
0146 Result<SympyStepper::Jacobian> SympyStepper::transportToBound(
0147     State& state, const Surface& surface,
0148     const FreeToBoundCorrection& freeToBoundCorrection) const {
0149   if (!surface.isOnSurface(state.options.geoContext, position(state),
0150                            direction(state), BoundaryTolerance::Infinite())) {
0151     return Result<Jacobian>::failure(SurfaceError::GlobalPositionNotOnSurface);
0152   }
0153
0154   Jacobian jacobian = Jacobian::Identity();
0155   if (!state.cov.has_value()) {
0156     state.materialEffectsAccumulator.reset();
0157     return Result<Jacobian>::success(jacobian);
0158   }
0159   const std::optional<FreeMatrix> additionalFreeCovariance =
0160       state.materialEffectsAccumulator.computeAdditionalFreeCovariance(
0161           direction(state));
0162   state.materialEffectsAccumulator.reset();
0163   detail::sympy::fromScaledBoundToFree(state.jacToGlobal, qOverP(state));
0164   detail::sympy::transportCovarianceToBound(
0165       state.options.geoContext, surface, *state.cov, jacobian, state.derivative,
0166       state.jacToGlobal, additionalFreeCovariance, state.pars,
0167       freeToBoundCorrection);
0168   detail::sympy::toScaledBoundToFree(state.jacToGlobal, qOverP(state));
0169   return Result<Jacobian>::success(jacobian);
0170 }
0171
0172 Result<double> SympyStepper::step(State& state, Direction propDir,
0173                                   const IVolumeMaterial* material) const {
0174   if (state.options.doDense &&
0175       (material != nullptr || !state.materialEffectsAccumulator.isVacuum())) {
0176     if (state.cov.has_value()) {
0177       return detail::sympyStep<detail::SympyStepMode::Dense, true>(
0178           *this, state, propDir, material);
0179     }
0180     return detail::sympyStep<detail::SympyStepMode::Dense, false>(
0181         *this, state, propDir, material);
0182   }
0183   if (state.cov.has_value()) {
0184     return detail::sympyStep<detail::SympyStepMode::Vacuum, true>(
0185         *this, state, propDir, nullptr);
0186   }
0187   return detail::sympyStep<detail::SympyStepMode::Vacuum, false>(
0188       *this, state, propDir, nullptr);
0189 }
0190
0191 }  // namespace Acts