Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-01-18 09:27:44

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