em/xs/WentzelHelper.hh

0001 //----------------------------------*-C++-*----------------------------------//
0002 // Copyright 2023-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/WentzelHelper.hh
0007 //---------------------------------------------------------------------------//
0008 #pragma once
0009
0010 #include "corecel/Macros.hh"
0011 #include "corecel/Types.hh"
0012 #include "corecel/math/Algorithms.hh"
0013 #include "celeritas/em/data/CommonCoulombData.hh"
0014 #include "celeritas/em/data/WentzelOKVIData.hh"
0015 #include "celeritas/mat/MaterialView.hh"
0016 #include "celeritas/phys/AtomicNumber.hh"
0017 #include "celeritas/phys/ParticleTrackView.hh"
0018
0019 namespace celeritas
0020 {
0021 //---------------------------------------------------------------------------//
0022 /*!
0023  * Helper class for the Wentzel OK and VI Coulomb scattering model.
0024  *
0025  * This calculates the Moliere screening coefficient, the maximum scattering
0026  * angle off of electrons, and the ratio of the electron to total Wentzel cross
0027  * sections.
0028  *
0029  * References:
0030  * [PRM] Geant4 Physics Reference Manual (Release 11.1) section 8.5.
0031  */
0032 class WentzelHelper
0033 {
0034   public:
0035     //!@{
0036     //! \name Type aliases
0037     using Charge = units::ElementaryCharge;
0038     using Energy = units::MevEnergy;
0039     using Mass = units::MevMass;
0040     using MomentumSq = units::MevMomentumSq;
0041     //!@}
0042
0043   public:
0044     // Construct from particle and material properties
0045     inline CELER_FUNCTION
0046     WentzelHelper(ParticleTrackView const& particle,
0047                   MaterialView const& material,
0048                   AtomicNumber target_z,
0049                   NativeCRef<WentzelOKVIData> const& wentzel,
0050                   CoulombIds const& ids,
0051                   Energy cutoff);
0052
0053     //! Get the target atomic number
0054     CELER_FUNCTION AtomicNumber atomic_number() const { return target_z_; }
0055
0056     //! Get the Moliere screening coefficient
0057     CELER_FUNCTION real_type screening_coefficient() const
0058     {
0059         return screening_coefficient_;
0060     }
0061
0062     //! Get the Mott factor
0063     CELER_FUNCTION real_type mott_factor() const { return mott_factor_; }
0064
0065     //! Get the multiplicative factor for the cross section
0066     CELER_FUNCTION real_type kin_factor() const { return kin_factor_; }
0067
0068     //! Get the maximum scattering angle off of electrons
0069     CELER_FUNCTION real_type cos_thetamax_electron() const
0070     {
0071         return cos_thetamax_elec_;
0072     }
0073
0074     //! Get the maximum scattering angle off of a nucleus
0075     CELER_FUNCTION real_type cos_thetamax_nuclear() const
0076     {
0077         return cos_thetamax_nuc_;
0078     }
0079
0080     // Calculate the electron cross section for Coulomb scattering
0081     inline CELER_FUNCTION real_type
0082     calc_xs_electron(real_type cos_thetamin, real_type cos_thetamax) const;
0083
0084     // Calculate the nuclear cross section for Coulomb scattering
0085     inline CELER_FUNCTION real_type
0086     calc_xs_nuclear(real_type cos_thetamin, real_type cos_thetamax) const;
0087
0088   private:
0089     //// DATA ////
0090
0091     AtomicNumber const target_z_;
0092     real_type screening_coefficient_;
0093     real_type kin_factor_;
0094     real_type mott_factor_;
0095     real_type cos_thetamax_elec_;
0096     real_type cos_thetamax_nuc_;
0097
0098     //// HELPER FUNCTIONS ////
0099
0100     // Calculate the Moliere screening coefficient
0101     inline CELER_FUNCTION real_type
0102     calc_screening_coefficient(ParticleTrackView const& particle) const;
0103
0104     // Calculate the screening coefficient R^2 for electrons
0105     CELER_CONSTEXPR_FUNCTION real_type screen_r_sq_elec() const;
0106
0107     // Calculate the multiplicative factor for the cross section
0108     inline CELER_FUNCTION real_type
0109     calc_kin_factor(ParticleTrackView const&) const;
0110
0111     // Calculate the (cosine of) the maximum scattering angle off of electrons
0112     inline CELER_FUNCTION real_type calc_cos_thetamax_electron(
0113         ParticleTrackView const&, CoulombIds const&, Energy) const;
0114
0115     // Calculate the (cosine of) the maximum scattering angle off of a nucleus
0116     inline CELER_FUNCTION real_type calc_cos_thetamax_nuclear(
0117         ParticleTrackView const&,
0118         MaterialView const& material,
0119         NativeCRef<WentzelOKVIData> const& wentzel) const;
0120
0121     // Calculate the common factor in the electron and nuclear cross section
0122     inline CELER_FUNCTION real_type calc_xs_factor(real_type cos_thetamin,
0123                                                    real_type cos_thetamax) const;
0124 };
0125
0126 //---------------------------------------------------------------------------//
0127 // INLINE DEFINITIONS
0128 //---------------------------------------------------------------------------//
0129 /*!
0130  * Construct from particle and material properties.
0131  */
0132 CELER_FUNCTION
0133 WentzelHelper::WentzelHelper(ParticleTrackView const& particle,
0134                              MaterialView const& material,
0135                              AtomicNumber target_z,
0136                              NativeCRef<WentzelOKVIData> const& wentzel,
0137                              CoulombIds const& ids,
0138                              Energy cutoff)
0139     : target_z_(target_z)
0140     , screening_coefficient_(this->calc_screening_coefficient(particle)
0141                              * wentzel.params.screening_factor)
0142     , kin_factor_(this->calc_kin_factor(particle))
0143     , mott_factor_(particle.particle_id() == ids.electron
0144                        ? 1 + real_type(2e-4) * ipow<2>(target_z_.get())
0145                        : 1)
0146     , cos_thetamax_elec_(
0147           this->calc_cos_thetamax_electron(particle, ids, cutoff))
0148     , cos_thetamax_nuc_(
0149           this->calc_cos_thetamax_nuclear(particle, material, wentzel))
0150 {
0151     CELER_EXPECT(screening_coefficient_ > 0);
0152     CELER_EXPECT(cos_thetamax_elec_ >= -1 && cos_thetamax_elec_ <= 1);
0153     CELER_EXPECT(cos_thetamax_nuc_ >= -1 && cos_thetamax_nuc_ <= 1);
0154 }
0155
0156 //---------------------------------------------------------------------------//
0157 /*!
0158  * Calculate the electron cross section for Coulomb scattering.
0159  */
0160 CELER_FUNCTION real_type WentzelHelper::calc_xs_electron(
0161     real_type cos_thetamin, real_type cos_thetamax) const
0162 {
0163     cos_thetamin = max(cos_thetamin, cos_thetamax_elec_);
0164     cos_thetamax = max(cos_thetamax, cos_thetamax_elec_);
0165     if (cos_thetamin <= cos_thetamax)
0166     {
0167         return 0;
0168     }
0169     return this->calc_xs_factor(cos_thetamin, cos_thetamax);
0170 }
0171
0172 //---------------------------------------------------------------------------//
0173 /*!
0174  * Calculate the nuclear cross section for Coulomb scattering.
0175  */
0176 CELER_FUNCTION real_type WentzelHelper::calc_xs_nuclear(
0177     real_type cos_thetamin, real_type cos_thetamax) const
0178 {
0179     return target_z_.get() * this->calc_xs_factor(cos_thetamin, cos_thetamax);
0180 }
0181
0182 //---------------------------------------------------------------------------//
0183 /*!
0184  * Calculate the common factor in the electron and nuclear cross section.
0185  */
0186 CELER_FUNCTION real_type WentzelHelper::calc_xs_factor(
0187     real_type cos_thetamin, real_type cos_thetamax) const
0188 {
0189     return kin_factor_ * mott_factor_ * (cos_thetamin - cos_thetamax)
0190            / ((1 - cos_thetamin + 2 * screening_coefficient_)
0191               * (1 - cos_thetamax + 2 * screening_coefficient_));
0192 }
0193
0194 //---------------------------------------------------------------------------//
0195 /*!
0196  * Calculate the Moliere screening coefficient as in [PRM] eqn 8.51.
0197  *
0198  * \note The \c screenZ in Geant4 is equal to twice the screening coefficient.
0199  */
0200 CELER_FUNCTION real_type WentzelHelper::calc_screening_coefficient(
0201     ParticleTrackView const& particle) const
0202 {
0203     // TODO: Reference for just proton correction?
0204     real_type correction = 1;
0205     real_type sq_cbrt_z = fastpow(real_type(target_z_.get()), real_type{2} / 3);
0206     if (target_z_.get() > 1)
0207     {
0208         real_type tau = value_as<Energy>(particle.energy())
0209                         / value_as<Mass>(particle.mass());
0210         // TODO: Reference for this factor?
0211         real_type factor = std::sqrt(tau / (tau + sq_cbrt_z));
0212
0213         correction = min(target_z_.get() * real_type{1.13},
0214                          real_type{1.13}
0215                              + real_type{3.76}
0216                                    * ipow<2>(target_z_.get()
0217                                              * constants::alpha_fine_structure)
0218                                    * factor / particle.beta_sq());
0219     }
0220
0221     return correction * this->screen_r_sq_elec() * sq_cbrt_z
0222            / value_as<MomentumSq>(particle.momentum_sq());
0223 }
0224
0225 //---------------------------------------------------------------------------//
0226 /*!
0227  * Calculate the constant factor of the screening coefficient.
0228  *
0229  * This is the constant prefactor \f$ R^2 / Z^{2/3} \f$ of the screening
0230  * coefficient for incident electrons (equation 8.51 in [PRM]). The screening
0231  * radius \f$ R \f$ is given by:
0232  * \f[
0233    R = \frac{\hbar Z^{1/3}}{2C_{TF} a_0},
0234  * \f]
0235  * where the Thomas-Fermi constant \f$ C_{TF} \f$ is defined as
0236  * \f[
0237    C_{TF} = \frac{1}{2} \left(\frac{3\pi}{4}\right)^{2/3}.
0238  * \f]
0239  */
0240 CELER_CONSTEXPR_FUNCTION real_type WentzelHelper::screen_r_sq_elec() const
0241 {
0242     //! Thomas-Fermi constant \f$ C_{TF} \f$
0243     constexpr real_type ctf = 0.8853413770001135;
0244
0245     return native_value_to<MomentumSq>(
0246                ipow<2>(constants::hbar_planck / (2 * ctf * constants::a0_bohr)))
0247         .value();
0248 }
0249
0250 //---------------------------------------------------------------------------//
0251 /*!
0252  * Calculate the multiplicative factor for the cross section.
0253  *
0254  * This calculates the factor
0255  * \f[
0256    f = \frac{2 \pi m_e^2 r_e^2 Z q^2}{\beta^2 p^2},
0257  * \f]
0258  * where \f$ m_e, r_e, Z, q, \beta \f$, and \f$ p \f$ are the electron mass,
0259  * classical electron radius, atomic number of the target atom, charge,
0260  * relativistic speed, and momentum of the incident particle, respectively.
0261  */
0262 CELER_FUNCTION real_type
0263 WentzelHelper::calc_kin_factor(ParticleTrackView const& particle) const
0264 {
0265     real_type constexpr twopi_mrsq
0266         = 2 * constants::pi
0267           * ipow<2>(native_value_to<Mass>(constants::electron_mass).value()
0268                     * constants::r_electron);
0269
0270     return twopi_mrsq * target_z_.get()
0271            * ipow<2>(value_as<Charge>(particle.charge()))
0272            / (particle.beta_sq()
0273               * value_as<MomentumSq>(particle.momentum_sq()));
0274 }
0275
0276 //---------------------------------------------------------------------------//
0277 /*!
0278  * Calculate the maximum scattering angle off the target's electrons.
0279  *
0280  * This calculates the cosine of the maximum polar angle that the incident
0281  * particle can scatter off of the target's electrons.
0282  */
0283 CELER_FUNCTION real_type WentzelHelper::calc_cos_thetamax_electron(
0284     ParticleTrackView const& particle, CoulombIds const& ids, Energy cutoff) const
0285 {
0286     real_type inc_energy = value_as<Energy>(particle.energy());
0287     real_type mass = value_as<Mass>(particle.mass());
0288
0289     real_type max_energy = particle.particle_id() == ids.electron
0290                                ? real_type{0.5} * inc_energy
0291                                : inc_energy;
0292     real_type final_energy = inc_energy
0293                              - min(value_as<Energy>(cutoff), max_energy);
0294
0295     if (final_energy > 0)
0296     {
0297         real_type incident_ratio = 1 + 2 * mass / inc_energy;
0298         real_type final_ratio = 1 + 2 * mass / final_energy;
0299         real_type cos_t_max = std::sqrt(incident_ratio / final_ratio);
0300
0301         return clamp(cos_t_max, real_type{0}, real_type{1});
0302     }
0303     return 0;
0304 }
0305
0306 //---------------------------------------------------------------------------//
0307 /*!
0308  * Calculate the maximum scattering angle off the target nucleus.
0309  */
0310 CELER_FUNCTION real_type WentzelHelper::calc_cos_thetamax_nuclear(
0311     ParticleTrackView const& particle,
0312     MaterialView const& material,
0313     NativeCRef<WentzelOKVIData> const& wentzel) const
0314 {
0315     if (wentzel.params.is_combined)
0316     {
0317         CELER_ASSERT(material.material_id() < wentzel.inv_mass_cbrt_sq.size());
0318         return max(wentzel.params.costheta_limit,
0319                    1
0320                        - wentzel.params.a_sq_factor
0321                              * wentzel.inv_mass_cbrt_sq[material.material_id()]
0322                              / value_as<MomentumSq>(particle.momentum_sq()));
0323     }
0324     return wentzel.params.costheta_limit;
0325 }
0326
0327 //---------------------------------------------------------------------------//
0328 }  // namespace celeritas