em/xs/RBDiffXsCalculator.hh

0001 //----------------------------------*-C++-*----------------------------------//
0002 // Copyright 2021-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/RBDiffXsCalculator.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/em/data/RelativisticBremData.hh"
0020 #include "celeritas/em/interactor/detail/PhysicsConstants.hh"
0021
0022 #include "LPMCalculator.hh"
0023
0024 namespace celeritas
0025 {
0026 //---------------------------------------------------------------------------//
0027 /*!
0028  * Calculate differential cross sections for relativistic bremsstrahlung.
0029  *
0030  * This accounts for the LPM effect if the option is enabled and the
0031  * electron energy is high enough.
0032  *
0033  * This is a shape function used for rejection, so as long as the resulting
0034  * cross section is scaled by the maximum value the units do not matter.
0035  */
0036 class RBDiffXsCalculator
0037 {
0038   public:
0039     //!@{
0040     //! \name Type aliases
0041     using Energy = units::MevEnergy;
0042     using ElementData = RelBremElementData;
0043     //!@}
0044
0045   public:
0046     // Construct with incident electron and current element
0047     inline CELER_FUNCTION RBDiffXsCalculator(RelativisticBremRef const& shared,
0048                                              ParticleTrackView const& particle,
0049                                              MaterialView const& material,
0050                                              ElementComponentId elcomp_id);
0051
0052     // Compute cross section of exiting gamma energy
0053     inline CELER_FUNCTION real_type operator()(Energy energy);
0054
0055     //! Density correction factor [Energy^2]
0056     CELER_FUNCTION real_type density_correction() const
0057     {
0058         return density_corr_;
0059     }
0060
0061     //! Return the maximum value of the differential cross section
0062     CELER_FUNCTION real_type maximum_value() const
0063     {
0064         return elem_data_.factor1 + elem_data_.factor2;
0065     }
0066
0067   private:
0068     //// TYPES ////
0069
0070     //! Intermediate data for screening functions
0071     struct ScreenFunctions
0072     {
0073         real_type phi1{0};
0074         real_type phi2{0};
0075         real_type psi1{0};
0076         real_type psi2{0};
0077     };
0078
0079     using R = real_type;
0080
0081     //// DATA ////
0082
0083     // Element data of the current material
0084     ElementData const& elem_data_;
0085     // Shared problem data for the current material
0086     MaterialView const& material_;
0087     // Shared problem data for the current element
0088     ElementView const element_;
0089     // Total energy of the incident particle
0090     real_type total_energy_;
0091     // Density correction for the current material
0092     real_type density_corr_;
0093     // Flag for the LPM effect
0094     bool enable_lpm_;
0095     // Flag for dialectric suppression effect in LPM
0096     bool dielectric_suppression_;
0097
0098     //// HELPER FUNCTIONS ////
0099
0100     //! Calculate the differential cross section per atom
0101     inline CELER_FUNCTION real_type dxsec_per_atom(real_type energy);
0102
0103     //! Calculate the differential cross section per atom with the LPM effect
0104     inline CELER_FUNCTION real_type dxsec_per_atom_lpm(real_type energy);
0105
0106     //! Compute screening functions
0107     inline CELER_FUNCTION ScreenFunctions
0108     compute_screen_functions(real_type gamma, real_type epsilon);
0109 };
0110
0111 //---------------------------------------------------------------------------//
0112 // INLINE DEFINITIONS
0113 //---------------------------------------------------------------------------//
0114 /*!
0115  * Construct with incident electron and current element.
0116  */
0117 CELER_FUNCTION
0118 RBDiffXsCalculator::RBDiffXsCalculator(RelativisticBremRef const& shared,
0119                                        ParticleTrackView const& particle,
0120                                        MaterialView const& material,
0121                                        ElementComponentId elcomp_id)
0122     : elem_data_(shared.elem_data[material.element_id(elcomp_id)])
0123     , material_(material)
0124     , element_(material.make_element_view(elcomp_id))
0125     , total_energy_(value_as<Energy>(particle.total_energy()))
0126 {
0127     real_type density_factor = material.electron_density()
0128                                * detail::migdal_constant();
0129     density_corr_ = density_factor * ipow<2>(total_energy_);
0130
0131     real_type lpm_energy
0132         = material.radiation_length()
0133           * value_as<detail::MevPerLen>(detail::lpm_constant());
0134     real_type lpm_threshold = lpm_energy * std::sqrt(density_factor);
0135     enable_lpm_ = (shared.enable_lpm && (total_energy_ > lpm_threshold));
0136     dielectric_suppression_ = shared.dielectric_suppression();
0137 }
0138
0139 //---------------------------------------------------------------------------//
0140 /*!
0141  * Compute the relativistic differential cross section per atom at the given
0142  * bremsstrahlung photon energy in MeV.
0143  */
0144 CELER_FUNCTION
0145 real_type RBDiffXsCalculator::operator()(Energy energy)
0146 {
0147     CELER_EXPECT(energy > zero_quantity());
0148     return enable_lpm_ ? this->dxsec_per_atom_lpm(energy.value())
0149                        : this->dxsec_per_atom(energy.value());
0150 }
0151
0152 //---------------------------------------------------------------------------//
0153 /*!
0154  * Compute the differential cross section without the LPM effect.
0155  */
0156 CELER_FUNCTION
0157 real_type RBDiffXsCalculator::dxsec_per_atom(real_type gamma_energy)
0158 {
0159     real_type dxsec{0};
0160
0161     real_type y = gamma_energy / total_energy_;
0162     real_type onemy = 1 - y;
0163     real_type term0 = onemy + R(0.75) * ipow<2>(y);
0164
0165     if (element_.atomic_number() < AtomicNumber{5})
0166     {
0167         // The Dirac-Fock model
0168         dxsec = term0 * elem_data_.factor1 + onemy * elem_data_.factor2;
0169     }
0170     else
0171     {
0172         // Tsai's analytical approximation.
0173         real_type invz = 1
0174                          / static_cast<real_type>(
0175                              element_.atomic_number().unchecked_get());
0176         real_type term1 = y / (total_energy_ - gamma_energy);
0177         real_type gamma = term1 * elem_data_.gamma_factor;
0178         real_type epsilon = term1 * elem_data_.epsilon_factor;
0179
0180         // Evaluate the screening functions
0181         auto sfunc = compute_screen_functions(gamma, epsilon);
0182
0183         dxsec = term0
0184                     * ((R(0.25) * sfunc.phi1 - elem_data_.fz)
0185                        + (R(0.25) * sfunc.psi1 - 2 * element_.log_z() / 3)
0186                              * invz)
0187                 + R(0.125) * onemy * (sfunc.phi2 + sfunc.psi2 * invz);
0188     }
0189
0190     return celeritas::max(dxsec, R(0));
0191 }
0192
0193 //---------------------------------------------------------------------------//
0194 /*!
0195  * Compute the differential cross section with the LPM effect.
0196  */
0197 CELER_FUNCTION
0198 real_type RBDiffXsCalculator::dxsec_per_atom_lpm(real_type gamma_energy)
0199 {
0200     // Evaluate LPM functions
0201     real_type epsilon = total_energy_ / gamma_energy;
0202     LPMCalculator calc_lpm_functions(
0203         material_, element_, dielectric_suppression_, Energy{gamma_energy});
0204     auto lpm = calc_lpm_functions(epsilon);
0205
0206     real_type y = gamma_energy / total_energy_;
0207     real_type onemy = 1 - y;
0208     real_type y2 = R(0.25) * ipow<2>(y);
0209     real_type term = lpm.xi * (y2 * lpm.g + (onemy + 2 * y2) * lpm.phi);
0210     real_type dxsec = term * elem_data_.factor1 + onemy * elem_data_.factor2;
0211
0212     return max(dxsec, R(0));
0213 }
0214
0215 //---------------------------------------------------------------------------//
0216 /*!
0217  * Compute screen_functions: Tsai's analytical approximations of coherent and
0218  * incoherent screening function to the numerical screening functions computed
0219  * by using the Thomas-Fermi model: Y.-S.Tsai, Rev. Mod. Phys. 49 (1977) 421.
0220  */
0221 CELER_FUNCTION auto
0222 RBDiffXsCalculator::compute_screen_functions(real_type gam,
0223                                              real_type eps) -> ScreenFunctions
0224 {
0225     ScreenFunctions func;
0226     real_type gam2 = ipow<2>(gam);
0227     real_type eps2 = ipow<2>(eps);
0228
0229     func.phi1 = R(16.863) - 2 * std::log(1 + R(0.311877) * gam2)
0230                 + R(2.4) * std::exp(R(-0.9) * gam)
0231                 + R(1.6) * std::exp(R(-1.5) * gam);
0232     func.phi2 = 2 / (3 + R(19.5) * gam + 18 * gam2);
0233
0234     func.psi1 = R(24.34) - 2 * std::log(1 + R(13.111641) * eps2)
0235                 + R(2.8) * std::exp(R(-8) * eps)
0236                 + R(1.2) * std::exp(R(-29.2) * eps);
0237     func.psi2 = 2 / (3 + 120 * eps + 1200 * eps2);
0238
0239     return func;
0240 }
0241
0242 //---------------------------------------------------------------------------//
0243 }  // namespace celeritas