Propagator/detail/MaterialEffectsAccumulator.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/detail/MaterialEffectsAccumulator.hpp"
0010
0011 #include "Acts/Definitions/Direction.hpp"
0012 #include "Acts/Material/Interactions.hpp"
0013 #include "Acts/Utilities/MathHelpers.hpp"
0014
0015 namespace Acts::detail {
0016
0017 void MaterialEffectsAccumulator::initialize(
0018     double maxXOverX0Step, const ParticleHypothesis& particleHypothesis,
0019     double initialMomentum) {
0020   reset();
0021   m_maxXOverX0Step = maxXOverX0Step;
0022   m_particleHypothesis = particleHypothesis;
0023   m_initialMomentum = initialMomentum;
0024 }
0025
0026 void MaterialEffectsAccumulator::accumulate(const Material& material,
0027                                             double pathLength, double qOverPin,
0028                                             double qOverPout) {
0029   const Direction direction = Direction::fromScalarZeroAsPositive(pathLength);
0030   const MaterialSlab slab(material, std::abs(pathLength));
0031
0032   const float mass = m_particleHypothesis.mass();
0033   const float absQ = m_particleHypothesis.absoluteCharge();
0034   const PdgParticle absPdg = m_particleHypothesis.absolutePdg();
0035
0036   const double momentumIn = m_particleHypothesis.extractMomentum(qOverPin);
0037   const double momentumOut = m_particleHypothesis.extractMomentum(qOverPout);
0038
0039   const std::size_t substepCount =
0040       material.isVacuum() ? 1
0041                           : static_cast<std::size_t>(std::ceil(
0042                                 slab.thicknessInX0() / m_maxXOverX0Step));
0043   const double substep = pathLength / substepCount;
0044
0045   for (std::size_t i = 0; i < substepCount; ++i) {
0046     const double momentumMean =
0047         momentumIn + (momentumOut - momentumIn) * (i + 0.5) / substepCount;
0048     const double qOverPmean = m_particleHypothesis.qOverP(momentumMean, absQ);
0049
0050     const double theta0in = computeMultipleScatteringTheta0(
0051         m_accumulatedMaterial, absPdg, mass, qOverPmean, absQ);
0052
0053     m_accumulatedMaterial =
0054         MaterialSlab::combine(m_accumulatedMaterial, material, substep);
0055
0056     const double theta0out = computeMultipleScatteringTheta0(
0057         m_accumulatedMaterial, absPdg, mass, qOverPmean, absQ);
0058
0059     const double deltaVarTheta = square(theta0out) - square(theta0in);
0060     const double deltaVarPos =
0061         direction * m_varAngle * square(substep) +
0062         2 * m_covAnglePosition * substep +
0063         direction * deltaVarTheta * (square(substep) / 3);
0064     const double deltaCovAnglePosition =
0065         m_varAngle * substep + deltaVarTheta * substep / 2;
0066     m_varAngle += deltaVarTheta;
0067     m_varPosition += deltaVarPos;
0068     m_covAnglePosition += deltaCovAnglePosition;
0069   }
0070 }
0071
0072 std::optional<FreeMatrix>
0073 MaterialEffectsAccumulator::computeAdditionalFreeCovariance(
0074     const Vector3& direction) {
0075   if (isVacuum()) {
0076     return std::nullopt;
0077   }
0078
0079   FreeMatrix additionalFreeCovariance = FreeMatrix::Zero();
0080
0081   // handle multiple scattering
0082   {
0083     // for derivation see
0084     // https://github.com/andiwand/cern-scripts/blob/5f0ebf1bef35db65322f28c2e840c1db1aaaf9a7/notebooks/2023-12-07_qp-dense-nav.ipynb
0085     //
0086     SquareMatrix3 directionProjection =
0087         (ActsSquareMatrix<3>::Identity() - direction * direction.transpose());
0088
0089     additionalFreeCovariance.template block<3, 3>(eFreeDir0, eFreeDir0) =
0090         m_varAngle * directionProjection;
0091     additionalFreeCovariance.template block<3, 3>(eFreePos0, eFreePos0) =
0092         m_varPosition * directionProjection;
0093     additionalFreeCovariance.template block<3, 3>(eFreePos0, eFreeDir0) =
0094         m_covAnglePosition * directionProjection;
0095     additionalFreeCovariance.template block<3, 3>(eFreeDir0, eFreePos0) =
0096         additionalFreeCovariance.template block<3, 3>(eFreePos0, eFreeDir0);
0097   }
0098
0099   // handle energy loss covariance
0100   {
0101     const double mass = m_particleHypothesis.mass();
0102     const double absQ = m_particleHypothesis.absoluteCharge();
0103     const double qOverP = m_particleHypothesis.qOverP(
0104         m_initialMomentum, m_particleHypothesis.absoluteCharge());
0105
0106     const double qOverPSigma = computeEnergyLossLandauSigmaQOverP(
0107         m_accumulatedMaterial, mass, qOverP, absQ);
0108
0109     additionalFreeCovariance(eFreeQOverP, eFreeQOverP) =
0110         qOverPSigma * qOverPSigma;
0111
0112     // in principle the energy loss uncertainty also affects the time
0113     // uncertainty continuously. these terms are not included here.
0114   }
0115
0116   return additionalFreeCovariance;
0117 }
0118
0119 }  // namespace Acts::detail