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