Propagator/detail/PointwiseMaterialInteraction.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/Algebra.hpp"
0012 #include "Acts/Definitions/Common.hpp"
0013 #include "Acts/Definitions/Direction.hpp"
0014 #include "Acts/Definitions/TrackParametrization.hpp"
0015 #include "Acts/Definitions/Units.hpp"
0016 #include "Acts/EventData/ParticleHypothesis.hpp"
0017 #include "Acts/Geometry/GeometryContext.hpp"
0018 #include "Acts/Material/MaterialSlab.hpp"
0019 #include "Acts/Surfaces/Surface.hpp"
0020 #include "Acts/Utilities/MathHelpers.hpp"
0021
0022 #include <algorithm>
0023 #include <cmath>
0024
0025 namespace Acts::detail {
0026
0027 /// Determine the material update mode to be used for a given surface. Depending
0028 /// on whether the surface is the start or target surface, the update mode is
0029 /// restricted to pre- or post-update only. This is necessary to avoid double
0030 /// counting of material effects at the start and target surfaces.
0031 /// @param surface The current surface
0032 /// @param startSurface The starting surface of the propagation
0033 /// @param targetSurface The target surface of the propagation
0034 /// @param requestedMode The requested material update mode
0035 /// @return The determined material update mode
0036 MaterialUpdateMode determineMaterialUpdateMode(
0037     const Surface& surface, const Surface* startSurface,
0038     const Surface* targetSurface, MaterialUpdateMode requestedMode);
0039
0040 /// Determine the material update mode to be used for the current surface given
0041 /// the propagation state. Depending on whether the surface is the start or
0042 /// target surface, the update mode is restricted to pre- or post-update only.
0043 /// This is necessary to avoid double counting of material effects at the start
0044 /// and target surfaces.
0045 /// @tparam propagator_state_t The type of the propagator state
0046 /// @tparam navigator_t The type of the navigator
0047 /// @param state The current propagator state
0048 /// @param navigator The navigator used for the propagation
0049 /// @param requestedMode The requested material update mode
0050 /// @return The determined material update mode
0051 template <typename propagator_state_t, typename navigator_t>
0052 MaterialUpdateMode determineMaterialUpdateMode(
0053     const propagator_state_t& state, const navigator_t& navigator,
0054     MaterialUpdateMode requestedMode) {
0055   const Surface* surface = navigator.currentSurface(state.navigation);
0056   const Surface* startSurface = navigator.startSurface(state.navigation);
0057   const Surface* targetSurface = navigator.targetSurface(state.navigation);
0058
0059   if (surface == nullptr) {
0060     return MaterialUpdateMode::NoUpdate;
0061   }
0062
0063   return determineMaterialUpdateMode(*surface, startSurface, targetSurface,
0064                                      requestedMode);
0065 }
0066
0067 /// Evaluate the material slab at a given surface, position, and direction
0068 /// @param geoContext The geometry context
0069 /// @param surface The surface at which to evaluate the material slab
0070 /// @param propagationDirection The propagation direction
0071 /// @param position The position at which to evaluate the material slab
0072 /// @param direction The direction at which to evaluate the material slab
0073 /// @param updateMode The material update mode
0074 /// @return The evaluated material slab
0075 Result<MaterialSlab> evaluateMaterialSlab(const GeometryContext& geoContext,
0076                                           const Surface& surface,
0077                                           Direction propagationDirection,
0078                                           const Vector3& position,
0079                                           const Vector3& direction,
0080                                           MaterialUpdateMode updateMode);
0081
0082 /// Evaluate the material slab at the propagation state and surface
0083 /// @tparam propagator_state_t The type of the propagator state
0084 /// @tparam stepper_t The type of the stepper
0085 /// @param state The current propagator state
0086 /// @param stepper The stepper used for the propagation
0087 /// @param surface The surface at which to evaluate the material slab
0088 /// @param updateMode The material update mode
0089 /// @return The evaluated material slab
0090 template <typename propagator_state_t, typename stepper_t>
0091 Result<MaterialSlab> evaluateMaterialSlab(const propagator_state_t& state,
0092                                           const stepper_t& stepper,
0093                                           const Surface& surface,
0094                                           MaterialUpdateMode updateMode) {
0095   const GeometryContext& geoContext = state.options.geoContext;
0096   const Direction propagationDirection = state.options.direction;
0097   const Vector3 position = stepper.position(state.stepping);
0098   const Vector3 direction = stepper.direction(state.stepping);
0099
0100   return evaluateMaterialSlab(geoContext, surface, propagationDirection,
0101                               position, direction, updateMode);
0102 }
0103
0104 /// Struct to hold the material effects computed at a pointwise interaction
0105 struct PointwiseMaterialEffects {
0106   double eLoss = 0;
0107   double variancePhi = 0;
0108   double varianceTheta = 0;
0109   double varianceQoverP = 0;
0110 };
0111
0112 /// Compute the material effects given a material slab and particle properties
0113 /// @param slab The material slab
0114 /// @param particleHypothesis The particle hypothesis
0115 /// @param direction The direction of the particle
0116 /// @param qOverP The charge over momentum of the particle
0117 /// @param multipleScattering Whether to compute multiple scattering effects
0118 /// @param energyLoss Whether to compute energy loss effects
0119 /// @param covTransport Whether to compute covariance transport effects
0120 /// @return The computed material effects
0121 PointwiseMaterialEffects computeMaterialEffects(
0122     const MaterialSlab& slab, const ParticleHypothesis& particleHypothesis,
0123     const Vector3& direction, float qOverP, bool multipleScattering,
0124     bool energyLoss, bool covTransport);
0125
0126 /// Compute the material effects given the propagation state and material slab
0127 /// @tparam propagator_state_t The type of the propagator state
0128 /// @tparam stepper_t The type of the stepper
0129 /// @param state The current propagator state
0130 /// @param stepper The stepper used for the propagation
0131 /// @param slab The material slab
0132 /// @param multipleScattering Whether to compute multiple scattering effects
0133 /// @param energyLoss Whether to compute energy loss effects
0134 /// @return The computed material effects
0135 template <typename propagator_state_t, typename stepper_t>
0136 PointwiseMaterialEffects computeMaterialEffects(const propagator_state_t& state,
0137                                                 const stepper_t& stepper,
0138                                                 const MaterialSlab& slab,
0139                                                 bool multipleScattering,
0140                                                 bool energyLoss) {
0141   const bool covTransport = state.stepping.covTransport;
0142   const Vector3 direction = stepper.direction(state.stepping);
0143   const float qOverP = stepper.qOverP(state.stepping);
0144   const ParticleHypothesis& particleHypothesis =
0145       stepper.particleHypothesis(state.stepping);
0146
0147   return computeMaterialEffects(slab, particleHypothesis, direction, qOverP,
0148                                 multipleScattering, energyLoss, covTransport);
0149 }
0150
0151 /// Perform the material interaction given the propagation state and material
0152 /// slab
0153 /// @tparam propagator_state_t The type of the propagator state
0154 /// @tparam stepper_t The type of the stepper
0155 /// @param state The current propagator state
0156 /// @param stepper The stepper used for the propagation
0157 /// @param slab The material slab
0158 /// @param noiseUpdateMode The noise update mode
0159 /// @param multipleScattering Whether to compute multiple scattering effects
0160 /// @param energyLoss Whether to compute energy loss effects
0161 /// @return The computed material effects
0162 template <typename propagator_state_t, typename stepper_t>
0163 PointwiseMaterialEffects performMaterialInteraction(
0164     propagator_state_t& state, const stepper_t& stepper,
0165     const MaterialSlab& slab, NoiseUpdateMode noiseUpdateMode,
0166     bool multipleScattering, bool energyLoss) {
0167   if (slab.isVacuum()) {
0168     return {};
0169   }
0170
0171   const PointwiseMaterialEffects effects = computeMaterialEffects(
0172       state, stepper, slab, multipleScattering, energyLoss);
0173
0174   const Direction propDir = state.options.direction;
0175   const ParticleHypothesis& particleHypothesis =
0176       stepper.particleHypothesis(state.stepping);
0177   const double mass = particleHypothesis.mass();
0178   const double absQ = particleHypothesis.absoluteCharge();
0179   const Vector3 position = stepper.position(state.stepping);
0180   const double time = stepper.time(state.stepping);
0181   const Vector3 direction = stepper.direction(state.stepping);
0182   const float qOverP = stepper.qOverP(state.stepping);
0183   const double momentum = stepper.absoluteMomentum(state.stepping);
0184
0185   //! [energy loss update]
0186   // in forward(backward) propagation, energy decreases(increases) and
0187   // variances increase(decrease)
0188   const double nextE = fastHypot(mass, momentum) - effects.eLoss * propDir;
0189   // put particle at rest if energy loss is too large
0190   double nextP = (mass < nextE) ? fastCathetus(nextE, mass) : 0;
0191
0192   // minimum momentum below which we will not push particles via material
0193   // update
0194   // TODO 10 MeV might be quite low and we should make this configurable
0195   static constexpr double minP = 10 * Acts::UnitConstants::MeV;
0196   nextP = std::max(minP, nextP);
0197   const double nextQOverP =
0198       particleHypothesis.qOverP(nextP, std::copysign(absQ, qOverP));
0199
0200   // update track parameters
0201   stepper.update(state.stepping, position, direction, nextQOverP, time);
0202   //! [energy loss update]
0203
0204   // Convenience method to update a variance given a change and noise update
0205   // mode
0206   const auto updateVariance = [](double variance, double change,
0207                                  NoiseUpdateMode updateMode) {
0208     // Add/Subtract the change
0209     // Protect the variance against becoming negative
0210     return std::max(0.,
0211                     variance + std::copysign(change, toUnderlying(updateMode)));
0212   };
0213
0214   // update covariance matrix
0215   //! [covariance update]
0216   state.stepping.cov(eBoundPhi, eBoundPhi) =
0217       updateVariance(state.stepping.cov(eBoundPhi, eBoundPhi),
0218                      effects.variancePhi, noiseUpdateMode);
0219   state.stepping.cov(eBoundTheta, eBoundTheta) =
0220       updateVariance(state.stepping.cov(eBoundTheta, eBoundTheta),
0221                      effects.varianceTheta, noiseUpdateMode);
0222   state.stepping.cov(eBoundQOverP, eBoundQOverP) =
0223       updateVariance(state.stepping.cov(eBoundQOverP, eBoundQOverP),
0224                      effects.varianceQoverP, noiseUpdateMode);
0225   //! [covariance update]
0226
0227   return effects;
0228 }
0229
0230 /// Perform the material interaction at the current surface given the
0231 /// propagation state
0232 /// @tparam propagator_state_t The type of the propagator state
0233 /// @tparam stepper_t The type of the stepper
0234 /// @param state The current propagator state
0235 /// @param stepper The stepper used for the propagation
0236 /// @param surface The surface at which to perform the material interaction
0237 /// @param updateMode The material update mode
0238 /// @param noiseUpdateMode The noise update mode
0239 /// @param multipleScattering Whether to compute multiple scattering effects
0240 /// @param energyLoss Whether to compute energy loss effects
0241 /// @param logger The logger to use for verbose output
0242 /// @return The computed material effects
0243 template <typename propagator_state_t, typename stepper_t>
0244 Result<PointwiseMaterialEffects> performMaterialInteraction(
0245     propagator_state_t& state, const stepper_t& stepper, const Surface& surface,
0246     MaterialUpdateMode updateMode, NoiseUpdateMode noiseUpdateMode,
0247     bool multipleScattering, bool energyLoss, const Logger& logger) {
0248   const Result<MaterialSlab> slabResult =
0249       evaluateMaterialSlab(state, stepper, surface, updateMode);
0250   if (!slabResult.ok()) {
0251     return slabResult.error();
0252   }
0253   const MaterialSlab& slab = slabResult.value();
0254   if (slab.isVacuum()) {
0255     ACTS_VERBOSE("No material effects on surface: " << surface.geometryId()
0256                                                     << " with update mode: "
0257                                                     << updateMode);
0258   }
0259
0260   const PointwiseMaterialEffects effects = performMaterialInteraction(
0261       state, stepper, slab, noiseUpdateMode, multipleScattering, energyLoss);
0262
0263   const Direction propDir = state.options.direction;
0264
0265   ACTS_VERBOSE("Material effects on surface: "
0266                << surface.geometryId() << " with update mode: " << updateMode
0267                << " are :");
0268   ACTS_VERBOSE("eLoss = " << effects.eLoss * propDir << ", "
0269                           << "variancePhi = " << effects.variancePhi << ", "
0270                           << "varianceTheta = " << effects.varianceTheta << ", "
0271                           << "varianceQoverP = " << effects.varianceQoverP);
0272
0273   return effects;
0274 }
0275
0276 }  // namespace Acts::detail