EIC code displayed by LXR |
|
|
|
|||
File indexing completed on 2025-02-22 10:31:20
0001 //----------------------------------*-C++-*----------------------------------// 0002 // Copyright 2024 UT-Battelle, LLC, and other Celeritas developers. 0003 // See the top-level COPYRIGHT file for details. 0004 // SPDX-License-Identifier: (Apache-2.0 OR MIT) 0005 //---------------------------------------------------------------------------// 0006 //! \file celeritas/em/xs/MuBremsDiffXsCalculator.hh 0007 //---------------------------------------------------------------------------// 0008 #pragma once 0009 0010 #include <cmath> 0011 0012 #include "corecel/Macros.hh" 0013 #include "corecel/Types.hh" 0014 #include "corecel/math/Algorithms.hh" 0015 #include "corecel/math/Quantity.hh" 0016 #include "celeritas/Constants.hh" 0017 #include "celeritas/Quantities.hh" 0018 #include "celeritas/Types.hh" 0019 #include "celeritas/mat/ElementView.hh" 0020 0021 namespace celeritas 0022 { 0023 //---------------------------------------------------------------------------// 0024 /*! 0025 * Calculate the differential cross section for muon bremsstrahlung. 0026 * 0027 * The differential cross section can be written as 0028 * \f[ 0029 \difd{\sigma}{\epsilon} = \frac{16}{3} \alpha N_A (\frac{m}{\mu} 0030 r_e)^2 \frac{1}{\epsilon A} Z(Z \Phi_n + \Phi_e) (1 - v + \frac{3}{4} v^2), 0031 * \f] 0032 * where \f$ \epsilon \f$ is the photon energy, \f$ \alpha \f$ is the fine 0033 * structure constant, \f$ N_A \f$ is Avogadro's number, \f$ m \f$ is the 0034 * electron mass, \f$ \mu \f$ is the muon mass, \f$ r_e \f$ is the classical 0035 * electron radius, \f$ Z \f$ is the atomic number, and \f$ A \f$ is the atomic 0036 * mass. \f$ v = \epsilon / E \f$ is the relative energy transfer, where \f$ E 0037 * \f$ is the total energy of the incident muon. 0038 * 0039 * The contribution to the cross section from the nucleus is given by 0040 * \f[ 0041 \Phi_n = \ln \frac{B Z^{-1/3} (\mu + \delta(D'_n \sqrt{e} - 2))}{D'_n (m + 0042 \delta \sqrt{e} B Z^{-1/3})} \f$, 0043 * \f] 0044 * where \f$ \delta = \frac{\mu^2 v}{2(E - \epsilon)}\f$ is the minimum 0045 * momentum transfer and \f$ D'_n \f$ is the correction to the nuclear form 0046 * factor. 0047 * 0048 * The contribution to the cross section from electrons is given by 0049 * \f[ 0050 \Phi_e = \ln \frac{B' Z^{-2/3} \mu}{\left(1 + \frac{\delta \mu}{m^2 0051 \sqrt{e}}\right)(m + \delta \sqrt{e} B' Z^{-2/3})} \f$. 0052 * \f] 0053 * 0054 * The constants \f$ B \f$ and \f$ B' \f$ were calculated using the 0055 * Thomas-Fermi model. In the case of hydrogen, where the Thomas-Fermi model 0056 * does not serve as a good approximation, the exact values of the constants 0057 * were calculated analytically. 0058 * 0059 * This performs the same calculation as in Geant4's \c 0060 * G4MuBremsstrahlungModel::ComputeDMicroscopicCrossSection() and as described 0061 * in section 11.2.1 of the Physics Reference Manual. The formulae are taken 0062 * mainly from SR Kelner, RP Kokoulin, and AA Petrukhin. About cross section 0063 * for high-energy muon bremsstrahlung. Technical Report, MEphI, 1995. Preprint 0064 * MEPhI 024-95, Moscow, 1995, CERN SCAN-9510048. 0065 */ 0066 class MuBremsDiffXsCalculator 0067 { 0068 public: 0069 //!@{ 0070 //! \name Type aliases 0071 using Energy = units::MevEnergy; 0072 using Mass = units::MevMass; 0073 //!@} 0074 0075 public: 0076 // Construct with incident particle data and current element 0077 inline CELER_FUNCTION MuBremsDiffXsCalculator(ElementView const& element, 0078 Energy inc_energy, 0079 Mass inc_mass, 0080 Mass electron_mass); 0081 0082 // Compute cross section of exiting gamma energy 0083 inline CELER_FUNCTION real_type operator()(Energy energy); 0084 0085 private: 0086 //// DATA //// 0087 0088 // Atomic number of the current element 0089 int atomic_number_; 0090 // Atomic mass of the current element 0091 real_type atomic_mass_; 0092 // \f$ Z^{-1/3} \f$ 0093 real_type inv_cbrt_z_; 0094 // Energy of the incident particle 0095 real_type inc_energy_; 0096 // Mass of the incident particle 0097 real_type inc_mass_; 0098 // Square of the incident particle mass 0099 real_type inc_mass_sq_; 0100 // Total energy of the incident particle 0101 real_type total_energy_; 0102 // Mass of an electron 0103 real_type electron_mass_; 0104 // Correction to the nuclear form factor 0105 real_type d_n_; 0106 // Constant in the radiation logarithm 0107 real_type b_; 0108 // Constant in the inelastic radiation logarithm 0109 real_type b_prime_; 0110 }; 0111 0112 //---------------------------------------------------------------------------// 0113 // INLINE DEFINITIONS 0114 //---------------------------------------------------------------------------// 0115 /*! 0116 * Construct with incident particle data and current element. 0117 */ 0118 CELER_FUNCTION 0119 MuBremsDiffXsCalculator::MuBremsDiffXsCalculator(ElementView const& element, 0120 Energy inc_energy, 0121 Mass inc_mass, 0122 Mass electron_mass) 0123 : atomic_number_(element.atomic_number().unchecked_get()) 0124 , atomic_mass_(value_as<units::AmuMass>(element.atomic_mass())) 0125 , inv_cbrt_z_(1 / element.cbrt_z()) 0126 , inc_energy_(value_as<Energy>(inc_energy)) 0127 , inc_mass_(value_as<Mass>(inc_mass)) 0128 , inc_mass_sq_(ipow<2>(inc_mass_)) 0129 , total_energy_(inc_energy_ + inc_mass_) 0130 , electron_mass_(value_as<Mass>(electron_mass)) 0131 { 0132 CELER_EXPECT(inc_energy_ > 0); 0133 0134 d_n_ = real_type(1.54) * std::pow(atomic_mass_, real_type(0.27)); 0135 if (atomic_number_ == 1) 0136 { 0137 // Constants calculated calculated analytically 0138 b_ = real_type(202.4); 0139 b_prime_ = 446; 0140 } 0141 else 0142 { 0143 // Constants calculated using the Thomas-Fermi model 0144 b_ = 183; 0145 b_prime_ = 1429; 0146 d_n_ = std::pow(d_n_, 1 - real_type(1) / atomic_number_); 0147 } 0148 } 0149 0150 //---------------------------------------------------------------------------// 0151 /*! 0152 * Compute the differential cross section per atom at the given photon energy. 0153 */ 0154 CELER_FUNCTION 0155 real_type MuBremsDiffXsCalculator::operator()(Energy energy) 0156 { 0157 CELER_EXPECT(energy > zero_quantity()); 0158 0159 if (value_as<Energy>(energy) >= inc_energy_) 0160 { 0161 return 0; 0162 } 0163 0164 // Calculate the relative energy transfer 0165 real_type v = value_as<Energy>(energy) / total_energy_; 0166 0167 // Calculate the minimum momentum transfer 0168 real_type delta = real_type(0.5) * inc_mass_sq_ * v 0169 / (total_energy_ - value_as<Energy>(energy)); 0170 0171 // Calculate the contribution to the cross section from the nucleus 0172 real_type sqrt_euler = std::sqrt(constants::euler); 0173 real_type phi_n = clamp_to_nonneg(std::log( 0174 b_ * inv_cbrt_z_ * (inc_mass_ + delta * (d_n_ * sqrt_euler - 2)) 0175 / (d_n_ * (electron_mass_ + delta * sqrt_euler * b_ * inv_cbrt_z_)))); 0176 0177 // Photon energy above which there is no contribution from electrons 0178 real_type energy_max_prime = total_energy_ 0179 / (1 0180 + real_type(0.5) * inc_mass_sq_ 0181 / (electron_mass_ * total_energy_)); 0182 0183 // Calculate the contribution to the cross section from electrons 0184 real_type phi_e = 0; 0185 if (value_as<Energy>(energy) < energy_max_prime) 0186 { 0187 real_type inv_cbrt_z_sq = ipow<2>(inv_cbrt_z_); 0188 phi_e = clamp_to_nonneg(std::log( 0189 b_prime_ * inv_cbrt_z_sq * inc_mass_ 0190 / ((1 + delta * inc_mass_ / (ipow<2>(electron_mass_) * sqrt_euler)) 0191 * (electron_mass_ 0192 + delta * sqrt_euler * b_prime_ * inv_cbrt_z_sq)))); 0193 } 0194 0195 // Calculate the differential cross section 0196 return 16 * constants::alpha_fine_structure * constants::na_avogadro 0197 * ipow<2>(electron_mass_ * constants::r_electron) * atomic_number_ 0198 * (atomic_number_ * phi_n + phi_e) 0199 * (1 - v * (1 - real_type(0.75) * v)) 0200 / (3 * inc_mass_sq_ * value_as<Energy>(energy) * atomic_mass_); 0201 } 0202 0203 //---------------------------------------------------------------------------// 0204 } // namespace celeritas
| [ Source navigation ] | [ Diff markup ] | [ Identifier search ] | [ general search ] |
|
This page was automatically generated by the 2.3.7 LXR engine. The LXR team |
|