em/xs/WentzelHelper.hh

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