Propagator/detail/SympyStepperDenseStepImpl.hpp

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 #pragma once
0010
0011 #include "Acts/Definitions/PdgParticle.hpp"
0012 #include "Acts/Definitions/Units.hpp"
0013 #include "Acts/Material/IVolumeMaterial.hpp"
0014 #include "Acts/Material/Interactions.hpp"
0015 #include "Acts/Material/MaterialSlab.hpp"
0016 #include "Acts/Propagator/detail/SympyStepperDenseStep.hpp"
0017
0018 #include <cmath>
0019
0020 #include "codegen/sympy_stepper_math.hpp"
0021
0022 namespace Acts {
0023
0024 template <bool WithJac>
0025 detail::Rk4Status detail::sympyDenseStep(
0026     const SympyStepper& stepper, SympyStepper::State& state,
0027     const IVolumeMaterial& material, double h, double errTol,
0028     double& errorEstimate, Vector3& lastField, std::error_code& fieldErr,
0029     std::span<double> jac) {
0030   const Vector3 pos = stepper.position(state);
0031   const Vector3 dir = stepper.direction(state);
0032   const double t = stepper.time(state);
0033   const double qop = stepper.qOverP(state);
0034   const double pabs = stepper.absoluteMomentum(state);
0035   const double m = stepper.particleHypothesis(state).mass();
0036   const double q = stepper.charge(state);
0037   const auto absQ = static_cast<float>(std::abs(q));
0038   const PdgParticle absPdg = stepper.particleHypothesis(state).absolutePdg();
0039
0040   const auto getB = [&](std::span<const double, 3> p) {
0041     return stepper.getField(state, {p[0], p[1], p[2]});
0042   };
0043
0044   const auto getG = [&](std::span<const double, 3> p, double l) -> double {
0045     if (const double newPabs =
0046             stepper.particleHypothesis(state).extractMomentum(l);
0047         newPabs < state.options.dense.momentumCutOff) {
0048       return 0.;
0049     }
0050
0051     const MaterialSlab slab(material.material({p[0], p[1], p[2]}),
0052                             1.0f * UnitConstants::mm);
0053     // Unsigned: the step's own sign gives the energy back on a backward step,
0054     // matching `PointwiseMaterialInteraction`.
0055     if (state.options.dense.meanEnergyLoss) {
0056       return computeEnergyLossMean(slab, absPdg, static_cast<float>(m),
0057                                    static_cast<float>(l), absQ);
0058     }
0059     return computeEnergyLossMode(slab, absPdg, static_cast<float>(m),
0060                                  static_cast<float>(l), absQ);
0061   };
0062
0063   if constexpr (WithJac) {
0064     return rk4_dense_jac(
0065         std::span<const double, 3>(pos.data(), 3),
0066         std::span<const double, 3>(dir.data(), 3), t, h, qop, m, q, pabs,
0067         std::span<const double, 3>(state.field->data(), 3), getB, getG,
0068         errorEstimate, errTol, fieldErr,
0069         std::span<double, 3>(state.pars.segment<3>(eFreePos0).data(), 3),
0070         state.pars[eFreeTime],
0071         std::span<double, 3>(state.pars.segment<3>(eFreeDir0).data(), 3),
0072         state.pars[eFreeQOverP], std::span<double, 3>(lastField.data(), 3),
0073         std::span<double, 8>(state.derivative.data(), 8), jac);
0074   } else {
0075     return rk4_dense_nojac(
0076         std::span<const double, 3>(pos.data(), 3),
0077         std::span<const double, 3>(dir.data(), 3), t, h, qop, m, q, pabs,
0078         std::span<const double, 3>(state.field->data(), 3), getB, getG,
0079         errorEstimate, errTol, fieldErr,
0080         std::span<double, 3>(state.pars.segment<3>(eFreePos0).data(), 3),
0081         state.pars[eFreeTime],
0082         std::span<double, 3>(state.pars.segment<3>(eFreeDir0).data(), 3),
0083         state.pars[eFreeQOverP], std::span<double, 3>(lastField.data(), 3));
0084   }
0085 }
0086
0087 }  // namespace Acts