material/detail/relativistic_quantities.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 // Detray include(s)
0012 #include "detray/definitions/algebra.hpp"
0013 #include "detray/definitions/detail/qualifiers.hpp"
0014 #include "detray/definitions/math.hpp"
0015 #include "detray/definitions/pdg_particle.hpp"
0016 #include "detray/definitions/units.hpp"
0017 #include "detray/material/material.hpp"
0018
0019 // System include(s)
0020 #include <cassert>
0021
0022 namespace detray::detail {
0023
0024 template <concepts::scalar scalar_t>
0025 struct relativistic_quantities {
0026   using scalar_type = scalar_t;
0027
0028   // values from RPP2018 table 33.1
0029   // Bethe formular prefactor. 1/mol unit is just a factor 1 here.
0030   static constexpr auto K{static_cast<scalar_type>(
0031       0.307075 * unit<double>::MeV * unit<double>::cm2)};
0032   // Energy scale for plasma energy.
0033   static constexpr auto PlasmaEnergyScale{
0034       static_cast<scalar_type>(28.816 * unit<double>::eV)};
0035
0036   scalar_type m_qOverP{0.f};
0037   scalar_type m_q2OverBeta2{0.f};
0038   scalar_type m_beta{0.f};
0039   scalar_type m_beta2{0.f};
0040   scalar_type m_betaGamma{0.f};
0041   scalar_type m_gamma{0.f};
0042   scalar_type m_gamma2{0.f};
0043   scalar_type m_E{0.f};
0044   scalar_type m_Wmax{0.f};
0045
0046   DETRAY_HOST_DEVICE
0047   relativistic_quantities(const pdg_particle<scalar_type>& ptc,
0048                           const scalar_type qop)
0049       : relativistic_quantities(ptc.mass(), qop, ptc.charge()) {};
0050
0051   DETRAY_HOST_DEVICE relativistic_quantities(const scalar_type mass,
0052                                              const scalar_type qOverP,
0053                                              const scalar_type q)
0054       : m_qOverP{qOverP},
0055         // beta²/q² = (p/E)²/q² = p²/(q²m² + q²p²) = 1/(q² + (m²(q/p)²)
0056         // q²/beta² = q² + m²(q/p)²
0057         m_q2OverBeta2{q * q + (mass * qOverP) * (mass * qOverP)} {
0058     assert(m_qOverP != 0.f);
0059     assert(mass != 0.f);
0060
0061     // 1/p = q/(qp) = (q/p)/q
0062     const scalar_type mOverP{
0063         mass * ((q != 0.f) ? math::fabs(qOverP / q) : math::fabs(qOverP))};
0064     const scalar_type pOverM{1.f / mOverP};
0065     // beta² = p²/E² = p²/(m² + p²) = 1/(1 + (m/p)²)
0066     m_beta2 = 1.f / (1.f + mOverP * mOverP);
0067     m_beta = math::sqrt(m_beta2);
0068     // beta*gamma = (p/sqrt(m² + p²))*(sqrt(m² + p²)/m) = p/m
0069     m_betaGamma = pOverM;
0070     // gamma = sqrt(m² + p²)/m = sqrt(1 + (p/m)²)
0071     m_gamma = math::sqrt(1.f + pOverM * pOverM);
0072     m_gamma2 = m_gamma * m_gamma;
0073     // E = gamma * mass;
0074     m_E = m_gamma * mass;
0075
0076     const scalar_type mfrac{constant<scalar_type>::m_e / mass};
0077
0078     // Wmax = 2m_e c^2 beta^2 gamma^2 / (1+2gamma*m_e/M + (m_e/M)^2)
0079     m_Wmax = (2.f * constant<scalar_type>::m_e * m_betaGamma * m_betaGamma) /
0080              (1.f + 2.f * m_gamma * mfrac + mfrac * mfrac);
0081   }
0082
0083   /// @return 2 * mass * (beta * gamma)² mass term.
0084   DETRAY_HOST_DEVICE constexpr scalar_type compute_mass_term(
0085       const scalar_type mass) const {
0086     return 2.f * mass * m_betaGamma * m_betaGamma;
0087   }
0088
0089   /// @return [(K/2) * (Z/A) * z^2 / beta^2 * density] in [energy/length]
0090   /// @brief defined in 34.12 of 2023 PDG review
0091   DETRAY_HOST_DEVICE constexpr scalar_type compute_epsilon_per_length(
0092       const material<scalar_type>& mat) const {
0093     return 0.5f * K * mat.molar_electron_density() * m_q2OverBeta2;
0094   }
0095
0096   /// @return  (K/2) * (Z/A) * z^2 / beta^2 * density * path_length
0097   /// @brief defined in 34.12 of 2023 PDG review
0098   DETRAY_HOST_DEVICE constexpr scalar_type compute_epsilon(
0099       const material<scalar_type>& mat, const scalar_type path_length) const {
0100     return compute_epsilon_per_length(mat) * path_length;
0101   }
0102
0103   /// @return the bethe_log_term of bethe equation,
0104   /// @brief bethe_log_term = 1/2 ln ( 2m_e c^2 beta^2 gamma^2 W_max/ I^2)
0105   DETRAY_HOST_DEVICE scalar_type
0106   compute_bethe_bloch_log_term(const material<scalar_type>& mat) const {
0107     const scalar_type I = mat.mean_excitation_energy();
0108
0109     assert(I != 0.f);
0110
0111     // u = 2 * m_e c^2* beta^2 * gamma^2
0112     const scalar_t u{compute_mass_term(constant<scalar_t>::m_e)};
0113     const scalar_type A = 0.5f * math::log(u * m_Wmax / (I * I));
0114     return A;
0115   }
0116
0117   /// @return d(bethe_log_term)/dqop
0118   /// @brief  dA/dqop = - 1 / (2 * qop) * [4 - W_max/ (gamma M c^2) ]
0119   DETRAY_HOST_DEVICE scalar_type derive_bethe_bloch_log_term() const {
0120     assert(m_gamma != 0.f);
0121     assert(m_E != 0.f);
0122     const scalar_type dAdqop = -1.f / (2.f * m_qOverP) * (4.f - m_Wmax / m_E);
0123     return dAdqop;
0124   }
0125
0126   /// @return d(beta^2)/dqop = - 2beta^2 / (qop * gamma^2)
0127   DETRAY_HOST_DEVICE scalar_type derive_beta2() const {
0128     assert(m_qOverP != 0.f);
0129     assert(m_gamma2 != 0.f);
0130     return -2.f * m_beta2 / (m_qOverP * m_gamma2);
0131   }
0132
0133   /// @return the half of density correction factor (delta/2).
0134   DETRAY_HOST_DEVICE inline scalar_type compute_delta_half(
0135       const material<scalar_type>& mat) const {
0136     if (!mat.has_density_effect_data()) {
0137       // Apply the approximated density effector correction to tracks with
0138       // beta*gamma > 10
0139       //
0140       // @NOTE The approximated function is only accurate for beta*gamma >
0141       // 100 therefore, the cutoff with 10 might not be a good choice. The
0142       // cutoff also introduces a step where the energy loss function is
0143       // not smooth anymore.
0144       //
0145       // For further discussion, please follow the ATLAS JIRA Tickets:
0146       // ATLASRECTS-3144 and ATLASRECTS-7586 (ATLAS-restricted)
0147       if (m_betaGamma < 10.f) {
0148         return 0.f;
0149       }
0150       // Equation 34.6 of PDG2022
0151       // @NOTE A factor of 1000 is required to convert the unit of density
0152       // (mm^-3 to cm^-3)
0153       const scalar_type plasmaEnergy{
0154           PlasmaEnergyScale *
0155           math::sqrt(1000.f * mat.molar_electron_density())};
0156       return math::log(m_betaGamma * plasmaEnergy /
0157                        mat.mean_excitation_energy()) -
0158              0.5f;
0159     } else {
0160       const auto& density = mat.density_effect_data();
0161
0162       const scalar_type cden{density.get_C_density()};
0163       const scalar_type mden{density.get_M_density()};
0164       const scalar_type aden{density.get_A_density()};
0165       const scalar_type x0den{density.get_X0_density()};
0166       const scalar_type x1den{density.get_X1_density()};
0167
0168       const scalar_type x{math::log10(m_betaGamma)};
0169
0170       scalar_type delta{0.f};
0171
0172       // From Geant4
0173       // processes/electromagnetic/lowenergy/src/G4hBetheBlochModel.cc
0174       if (x < x0den) {
0175         delta = 0.f;
0176         // @TODO: Add a branch for conductors (Eq 34.7 of
0177         // https://pdg.lbl.gov/2023/reviews/rpp2023-rev-particle-detectors-accel.pdf)
0178       } else {
0179         delta = 2.f * constant<scalar_type>::ln10 * x - cden;
0180         if (x < x1den) {
0181           delta += aden * math::pow((x1den - x), mden);
0182         }
0183       }
0184
0185       return 0.5f * delta;
0186     }
0187   }
0188
0189   /// @return the derivation of the density correction factor delta/2.
0190   DETRAY_HOST_DEVICE inline scalar_type derive_delta_half(
0191       const material<scalar_type>& mat) const {
0192     assert(m_qOverP != 0.f);
0193
0194     if (!mat.has_density_effect_data()) {
0195       // d(ln(betagamma))/dqop = -1/qop
0196       return -1.f / m_qOverP;
0197     } else {
0198       const auto& density = mat.density_effect_data();
0199
0200       // const scalar_type cden{density.get_C_density()};
0201       const scalar_type mden{density.get_M_density()};
0202       const scalar_type aden{density.get_A_density()};
0203       const scalar_type x0den{density.get_X0_density()};
0204       const scalar_type x1den{density.get_X1_density()};
0205
0206       const scalar_type x{math::log10(m_betaGamma)};
0207
0208       scalar_type delta{0.f};
0209
0210       // From Geant4
0211       // processes/electromagnetic/lowenergy/src/G4hBetheBlochModel.cc
0212       if (x < x0den) {
0213         delta = 0.f;
0214         // @TODO: Add a branch for conductors (Eq 34.7 of
0215         // https://pdg.lbl.gov/2023/reviews/rpp2023-rev-particle-detectors-accel.pdf)
0216       } else {
0217         delta = -2.f / m_qOverP;
0218         if (x < x1den) {
0219           delta += aden * mden / (m_qOverP * constant<scalar_type>::ln10) *
0220                    math::pow(x1den - x, mden - 1.f);
0221         }
0222       }
0223
0224       return 0.5f * delta;
0225     }
0226   }
0227 };
0228
0229 }  // namespace detray::detail