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