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/PdgParticle.hpp"
0015 #include "Acts/Definitions/TrackParametrization.hpp"
0016 #include "Acts/Definitions/Units.hpp"
0017 #include "Acts/Material/ISurfaceMaterial.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 /// @brief Struct to handle pointwise material interaction
0028 struct PointwiseMaterialInteraction {
0029   /// Data from the propagation state
0030   const Surface* surface = nullptr;
0031
0032   /// The particle position at the interaction.
0033   const Vector3 pos = Vector3(0., 0., 0);
0034   /// The particle time at the interaction.
0035   const double time = 0.0;
0036   /// The particle direction at the interaction.
0037   const Vector3 dir = Vector3(0., 0., 0);
0038   /// The particle q/p at the interaction
0039   const float qOverP = 0.0;
0040   /// The absolute particle charge
0041   const float absQ = 0.0;
0042   /// The particle momentum at the interaction
0043   const float momentum = 0.0;
0044   /// The particle mass
0045   const float mass = 0.0;
0046   /// The particle absolute pdg
0047   const PdgParticle absPdg = PdgParticle::eInvalid;
0048   /// The covariance transport decision at the interaction
0049   const bool performCovarianceTransport = false;
0050   /// The navigation direction
0051   const Direction navDir;
0052
0053   /// The effective, passed material properties including the path correction.
0054   MaterialSlab slab = MaterialSlab(0.);
0055   /// The path correction factor due to non-zero incidence on the surface.
0056   double pathCorrection = 0.;
0057   /// Expected phi variance due to the interactions.
0058   double variancePhi = 0.;
0059   /// Expected theta variance due to the interactions.
0060   double varianceTheta = 0.;
0061   /// Expected q/p variance due to the interactions.
0062   double varianceQoverP = 0.;
0063   /// The energy change due to the interaction.
0064   double Eloss = 0.;
0065   /// The momentum after the interaction
0066   double nextP = 0.;
0067
0068   /// @brief Constructor
0069   ///
0070   /// @tparam propagator_state_t Type of the propagator state
0071   /// @tparam stepper_t Type of the stepper
0072   ///
0073   /// @param [in] sSurface The current surface
0074   /// @param [in] state State of the propagation
0075   /// @param [in] stepper Stepper in use
0076   template <typename propagator_state_t, typename stepper_t>
0077   PointwiseMaterialInteraction(const Surface* sSurface,
0078                                const propagator_state_t& state,
0079                                const stepper_t& stepper)
0080       : surface(sSurface),
0081         pos(stepper.position(state.stepping)),
0082         time(stepper.time(state.stepping)),
0083         dir(stepper.direction(state.stepping)),
0084         qOverP(stepper.qOverP(state.stepping)),
0085         absQ(stepper.particleHypothesis(state.stepping).absoluteCharge()),
0086         momentum(stepper.absoluteMomentum(state.stepping)),
0087         mass(stepper.particleHypothesis(state.stepping).mass()),
0088         absPdg(stepper.particleHypothesis(state.stepping).absolutePdg()),
0089         performCovarianceTransport(state.stepping.covTransport),
0090         navDir(state.options.direction) {}
0091
0092   /// @brief This function evaluates the material properties to interact with
0093   /// This updates the slab and then returns, if the resulting slab is valid
0094   ///
0095   /// @tparam propagator_state_t Type of the propagator state
0096   /// @tparam navigator_t Type of the navigator
0097   ///
0098   /// @param [in] state State of the propagation
0099   /// @param [in] navigator Navigator of the propagation
0100   /// @param [in] updateStage The stage of the material update
0101   ///
0102   /// @return Boolean statement whether the material is valid
0103   template <typename propagator_state_t, typename navigator_t>
0104   bool evaluateMaterialSlab(
0105       const propagator_state_t& state, const navigator_t& navigator,
0106       MaterialUpdateStage updateStage = MaterialUpdateStage::FullUpdate) {
0107     // We are at the start surface
0108     if (surface == navigator.startSurface(state.navigation)) {
0109       updateStage = MaterialUpdateStage::PostUpdate;
0110       // Or is it the target surface ?
0111     } else if (surface == navigator.targetSurface(state.navigation)) {
0112       updateStage = MaterialUpdateStage::PreUpdate;
0113     }
0114
0115     // Retrieve the material properties
0116     slab = navigator.currentSurface(state.navigation)
0117                ->surfaceMaterial()
0118                ->materialSlab(pos, navDir, updateStage);
0119
0120     // Correct the material properties for non-zero incidence
0121     pathCorrection = surface->pathCorrection(state.geoContext, pos, dir);
0122     slab.scaleThickness(pathCorrection);
0123
0124     // Check if the evaluated material is valid
0125     return slab.isValid();
0126   }
0127
0128   /// @brief This function evaluate the material effects
0129   ///
0130   /// @param [in] multipleScattering Boolean to indicate the application of
0131   /// multiple scattering
0132   /// @param [in] energyLoss Boolean to indicate the application of energy loss
0133   void evaluatePointwiseMaterialInteraction(bool multipleScattering,
0134                                             bool energyLoss);
0135
0136   /// @brief Update the state
0137   ///
0138   /// @tparam propagator_state_t Type of the propagator state
0139   /// @tparam stepper_t Type of the stepper
0140   ///
0141   /// @param [in] state State of the propagation
0142   /// @param [in] stepper Stepper in use
0143   /// @param [in] updateMode The noise update mode (in default: add noise)
0144   template <typename propagator_state_t, typename stepper_t>
0145   void updateState(propagator_state_t& state, const stepper_t& stepper,
0146                    NoiseUpdateMode updateMode = addNoise) {
0147     const auto& particleHypothesis = stepper.particleHypothesis(state.stepping);
0148     // in forward(backward) propagation, energy decreases(increases) and
0149     // variances increase(decrease)
0150     const auto nextE = fastHypot(mass, momentum) - Eloss * navDir;
0151     // put particle at rest if energy loss is too large
0152     nextP = (mass < nextE) ? std::sqrt(nextE * nextE - mass * mass) : 0;
0153     // minimum momentum below which we will not push particles via material
0154     // update
0155     // TODO 10 MeV might be quite low and we should make this configurable
0156     static constexpr double minP = 10 * Acts::UnitConstants::MeV;
0157     nextP = std::max(minP, nextP);
0158     const double nextQOverP =
0159         std::copysign(particleHypothesis.qOverP(nextP, absQ), qOverP);
0160     // update track parameters and covariance
0161     stepper.update(state.stepping, pos, dir, nextQOverP, time);
0162     state.stepping.cov(eBoundPhi, eBoundPhi) = updateVariance(
0163         state.stepping.cov(eBoundPhi, eBoundPhi), variancePhi, updateMode);
0164     state.stepping.cov(eBoundTheta, eBoundTheta) =
0165         updateVariance(state.stepping.cov(eBoundTheta, eBoundTheta),
0166                        varianceTheta, updateMode);
0167     state.stepping.cov(eBoundQOverP, eBoundQOverP) =
0168         updateVariance(state.stepping.cov(eBoundQOverP, eBoundQOverP),
0169                        varianceQoverP, updateMode);
0170   }
0171
0172  private:
0173   /// @brief Evaluates the contributions to the covariance matrix
0174   ///
0175   /// @param [in] multipleScattering Boolean to indicate the application of
0176   /// multiple scattering
0177   /// @param [in] energyLoss Boolean to indicate the application of energy loss
0178   void covarianceContributions(bool multipleScattering, bool energyLoss);
0179
0180   /// @brief Convenience method for better readability
0181   ///
0182   /// @param [in] variance A diagonal entry of the covariance matrix
0183   /// @param [in] change The change that may be applied to it
0184   /// @param [in] updateMode The noise update mode (in default: add noise)
0185   ///
0186   /// @return The updated variance
0187   double updateVariance(double variance, double change,
0188                         NoiseUpdateMode updateMode = addNoise) const;
0189 };
0190
0191 }  // namespace Acts::detail