Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-09-17 08:53:34

0001 //------------------------------- -*- C++ -*- -------------------------------//
0002 // Copyright Celeritas contributors: see top-level COPYRIGHT file for details
0003 // SPDX-License-Identifier: (Apache-2.0 OR MIT)
0004 //---------------------------------------------------------------------------//
0005 //! \file celeritas/em/interactor/detail/SBPositronXsCorrector.hh
0006 //---------------------------------------------------------------------------//
0007 #pragma once
0008 
0009 #include <cmath>
0010 
0011 #include "corecel/Types.hh"
0012 #include "corecel/math/NumericLimits.hh"
0013 #include "celeritas/Constants.hh"
0014 #include "celeritas/Quantities.hh"
0015 #include "celeritas/mat/ElementView.hh"
0016 
0017 namespace celeritas
0018 {
0019 namespace detail
0020 {
0021 //---------------------------------------------------------------------------//
0022 /*!
0023  * Scale SB differential cross sections for positrons.
0024  *
0025  * This correction factor is the ratio of the positron to electron cross
0026  * sections. It appears in the bowels of \c
0027  * G4SeltzerBergerModel::SampleEnergyTransfer in Geant4 and scales the SB cross
0028  * section by a factor
0029  * \f[
0030    \frac{\sigma^+}{\sigma^-}
0031        = \exp( 2 \pi \alpha Z [ \beta^{-1}(E - k_c) - \beta^{-1}(E - k) ] )
0032  * \f]
0033  * where the inverse positron speed is : \f[
0034   \beta^{-1}(E) = \frac{c}{v} = \sqrt{1 - \left( \frac{m_e c^2}{E + m_e c^2}
0035   \right)^2}
0036   \,,
0037  * \f]
0038  * \f$ \alpha \f$ is the fine structure constant, \f$E\f$ is the
0039  * incident positron kinetic energy, \f$k_c\f$ is the gamma
0040  * production cutoff energy, and \f$ k \f$ is the provisionally sampled
0041  * energy of the emitted photon.
0042  *
0043  * \f$ \frac{\sigma^+}{\sigma^-} = 1 \f$ at the low end of the spectrum (where
0044  * \f$ k / E = 0 \f$) and is 0 at the tip of the spectrum where \f$ k / E \to 1
0045  * \f$.
0046  *
0047  * The correction factor is described in:
0048  *
0049  *   Kim, Longhuan, R. H. Pratt, S. M. Seltzer, and M. J. Berger. “Ratio of
0050  *   Positron to Electron Bremsstrahlung Energy Loss: An Approximate Scaling
0051  *   Law.” Physical Review A 33, no. 5 (May 1, 1986): 3002–9.
0052  *   https://doi.org/10.1103/PhysRevA.33.3002.
0053  */
0054 class SBPositronXsCorrector
0055 {
0056   public:
0057     //!@{
0058     using Energy = units::MevEnergy;
0059     using Mass = units::MevMass;
0060     using Xs = RealQuantity<SBElementTableData::XsUnits>;
0061     //!@}
0062 
0063   public:
0064     // Construct with positron data
0065     inline CELER_FUNCTION SBPositronXsCorrector(Mass positron_mass,
0066                                                 ElementView const& el,
0067                                                 Energy min_gamma_energy,
0068                                                 Energy inc_energy);
0069 
0070     // Calculate cross section scaling factor for the given exiting energy
0071     inline CELER_FUNCTION real_type operator()(Energy energy) const;
0072 
0073   private:
0074     //// DATA ////
0075 
0076     real_type const positron_mass_;
0077     real_type const alpha_z_;
0078     real_type const inc_energy_;
0079     real_type const cutoff_invbeta_;
0080 
0081     //// HELPER FUNCTIONS ////
0082 
0083     inline CELER_FUNCTION real_type calc_invbeta(real_type gamma_energy) const;
0084 };
0085 
0086 //---------------------------------------------------------------------------//
0087 // INLINE DEFINITIONS
0088 //---------------------------------------------------------------------------//
0089 /*!
0090  * Construct with positron data and energy range.
0091  */
0092 CELER_FUNCTION
0093 SBPositronXsCorrector::SBPositronXsCorrector(Mass positron_mass,
0094                                              ElementView const& el,
0095                                              Energy min_gamma_energy,
0096                                              Energy inc_energy)
0097     : positron_mass_{positron_mass.value()}
0098     , alpha_z_{2 * constants::pi * celeritas::constants::alpha_fine_structure
0099                * el.atomic_number().unchecked_get()}
0100     , inc_energy_(inc_energy.value())
0101     , cutoff_invbeta_{this->calc_invbeta(min_gamma_energy.value())}
0102 {
0103     CELER_EXPECT(inc_energy > min_gamma_energy);
0104 }
0105 
0106 //---------------------------------------------------------------------------//
0107 /*!
0108  * Calculate scaling factor for the given exiting gamma energy.
0109  *
0110  * Eq. 21 in Kim et al.
0111  */
0112 CELER_FUNCTION real_type SBPositronXsCorrector::operator()(Energy energy) const
0113 {
0114     CELER_EXPECT(energy > zero_quantity());
0115     CELER_EXPECT(energy.value() < inc_energy_);
0116     real_type delta = cutoff_invbeta_ - this->calc_invbeta(energy.value());
0117 
0118     // Avoid positive delta values due to floating point inaccuracies
0119     // See https://github.com/celeritas-project/celeritas/issues/617
0120     CELER_ASSERT(delta < 100 * numeric_limits<real_type>::epsilon());
0121     real_type result = std::exp(alpha_z_ * celeritas::min(delta, real_type{0}));
0122     CELER_ENSURE(result <= 1);
0123     return result;
0124 }
0125 
0126 //---------------------------------------------------------------------------//
0127 /*!
0128  * Calculate the inverse of the relativistic positron speed.
0129  *
0130  * The input here is the exiting gamma energy, so the positron energy is the
0131  * remainder from the incident energy. The relativistic speed \f$ \beta \f$
0132  * is:
0133  * \f[
0134   \beta^{-1}(K)
0135    = \frac{K + m c^2}{\sqrt{K (K + 2 m c^2)}}
0136    = \frac{K + mc^2}{\sqrt{K^2 + 2 K mc^2 + (mc^2)^2 - (mc^2)^2}}
0137    = \frac{K + mc^2}{\sqrt{(K + mc^2)^2 - mc^2}}
0138    = 1/\sqrt{1 - \left( \frac{mc^2}{K + mc^2} \right)^2}
0139    = 1 / \beta(K)
0140  * \f]
0141  */
0142 CELER_FUNCTION real_type
0143 SBPositronXsCorrector::calc_invbeta(real_type gamma_energy) const
0144 {
0145     CELER_EXPECT(gamma_energy > 0 && gamma_energy < inc_energy_);
0146     //! \todo Use local-data ParticleTrackView
0147     // Positron has all the energy except what it gave to the gamma
0148     real_type energy = inc_energy_ - gamma_energy;
0149     return (energy + positron_mass_)
0150            / std::sqrt(energy * (energy + 2 * positron_mass_));
0151 }
0152 
0153 //---------------------------------------------------------------------------//
0154 }  // namespace detail
0155 }  // namespace celeritas