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 WentzelHelper(
0052         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 calc_xs_electron(
0088         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 calc_xs_nuclear(
0092         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(
0129         real_type cos_thetamin, real_type cos_thetamax) const;
0130 };
0131
0132 //---------------------------------------------------------------------------//
0133 // INLINE DEFINITIONS
0134 //---------------------------------------------------------------------------//
0135 /*!
0136  * Construct from particle and material properties.
0137  */
0138 CELER_FUNCTION WentzelHelper::WentzelHelper(
0139     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     using namespace celeritas::literals;
0211
0212     // TODO: Reference for just proton correction?
0213     real_type correction = 1;
0214     real_type sq_cbrt_z = fastpow(real_type(target_z_.get()), 2_r / 3);
0215     if (target_z_.get() > 1)
0216     {
0217         // TODO: tau correction factor and "min" value are of unknown
0218         // provenance. The equation in Fernandez 1993 has factor=1, no special
0219         // casing for z=1, and no "min" for the correction
0220         real_type factor;
0221         real_type z_factor = 1;
0222         if (particle.particle_id() == ids.electron
0223             || particle.particle_id() == ids.positron)
0224         {
0225             // Electrons and positrons
0226             real_type tau = value_as<Energy>(particle.energy())
0227                             / value_as<Mass>(particle.mass());
0228             factor = std::sqrt(tau / (tau + sq_cbrt_z));
0229         }
0230         else
0231         {
0232             // Muons and hadrons
0233             factor = ipow<2>(value_as<Charge>(particle.charge()));
0234             z_factor += std::exp(-ipow<2>(target_z_.get()) * 0.001_r);
0235         }
0236         correction = min(target_z_.get() * 1.13_r,
0237                          1.13_r
0238                              + 3.76_r
0239                                    * ipow<2>(target_z_.get()
0240                                              * constants::alpha_fine_structure)
0241                                    * factor / particle.beta_sq())
0242                      * z_factor;
0243     }
0244
0245     return correction * sq_cbrt_z
0246            * value_as<MomentumSq>(this->screen_r_sq_elec())
0247            / value_as<MomentumSq>(particle.momentum_sq());
0248 }
0249
0250 //---------------------------------------------------------------------------//
0251 /*!
0252  * Calculate the constant factor of the screening coefficient.
0253  *
0254  * This is the constant prefactor \f$ R^2 / Z^{2/3} \f$ of the screening
0255  * coefficient for incident electrons (equation 8.51 in [PRM]). The screening
0256  * radius \f$ R \f$ is given by:
0257  * \f[
0258    R = \frac{\hbar Z^{1/3}}{2C_{TF} a_0},
0259  * \f]
0260  * where the Thomas-Fermi constant \f$ C_{TF} \f$ is defined as
0261  * \f[
0262    C_{TF} = \frac{1}{2} \left(\frac{3\pi}{4}\right)^{2/3}.
0263  * \f]
0264  */
0265 CELER_CONSTEXPR_FUNCTION auto WentzelHelper::screen_r_sq_elec() -> MomentumSq
0266 {
0267     //! Thomas-Fermi constant \f$ C_{TF} \f$
0268     constexpr Constant ctf{0.8853413770001135};
0269
0270     return native_value_to<MomentumSq>(
0271         ipow<2>(constants::hbar_planck / (2 * ctf * constants::a0_bohr)));
0272 }
0273
0274 //---------------------------------------------------------------------------//
0275 /*!
0276  * Calculate the multiplicative factor for the cross section.
0277  *
0278  * This calculates the factor
0279  * \f[
0280    f = \frac{2 \pi m_e^2 r_e^2 Z q^2}{\beta^2 p^2},
0281  * \f]
0282  * where \f$ m_e, r_e, Z, q, \beta \f$, and \f$ p \f$ are the electron mass,
0283  * classical electron radius, atomic number of the target atom, charge,
0284  * relativistic speed, and momentum of the incident particle, respectively.
0285  */
0286 CELER_FUNCTION real_type WentzelHelper::calc_kin_factor(
0287     ParticleTrackView const& particle, Mass electron_mass) const
0288 {
0289     constexpr Constant twopirsq = 2 * constants::pi
0290                                   * ipow<2>(constants::r_electron);
0291     return twopirsq * target_z_.get()
0292            * ipow<2>(value_as<Mass>(electron_mass)
0293                      * value_as<Charge>(particle.charge()))
0294            / (particle.beta_sq()
0295               * value_as<MomentumSq>(particle.momentum_sq()));
0296 }
0297
0298 //---------------------------------------------------------------------------//
0299 /*!
0300  * Calculate the maximum scattering angle off the target's electrons.
0301  *
0302  * This calculates the cosine of the maximum polar angle that the incident
0303  * particle can scatter off of the target's electrons.
0304  */
0305 CELER_FUNCTION real_type WentzelHelper::calc_cos_thetamax_electron(
0306     ParticleTrackView const& particle,
0307     CoulombIds const& ids,
0308     Energy cutoff,
0309     Mass electron_mass)
0310 {
0311     using namespace celeritas::literals;
0312
0313     real_type result = 0;
0314     real_type inc_energy = value_as<Energy>(particle.energy());
0315     real_type mass = value_as<Mass>(particle.mass());
0316
0317     if (particle.particle_id() == ids.electron
0318         || particle.particle_id() == ids.positron)
0319     {
0320         // Electrons and positrons
0321         real_type max_energy = particle.particle_id() == ids.electron
0322                                    ? 0.5_r * inc_energy
0323                                    : inc_energy;
0324         real_type final_energy = inc_energy
0325                                  - min(value_as<Energy>(cutoff), max_energy);
0326         if (final_energy > 0)
0327         {
0328             real_type inc_ratio = 1 + 2 * mass / inc_energy;
0329             real_type final_ratio = 1 + 2 * mass / final_energy;
0330             result = clamp<real_type>(std::sqrt(inc_ratio / final_ratio), 0, 1);
0331         }
0332     }
0333     else
0334     {
0335         // Muons and hadrons
0336         real_type mass_ratio = value_as<Mass>(electron_mass) / mass;
0337         real_type tau = inc_energy / mass;
0338         real_type max_energy
0339             = 2 * value_as<Mass>(electron_mass) * tau * (tau + 2)
0340               / (1 + 2 * mass_ratio * (tau + 1) + ipow<2>(mass_ratio));
0341         result = -min(value_as<Energy>(cutoff), max_energy)
0342                  * value_as<Mass>(electron_mass)
0343                  / value_as<MomentumSq>(particle.momentum_sq());
0344     }
0345     CELER_ENSURE(result >= 0 && result <= 1);
0346     return result;
0347 }
0348
0349 //---------------------------------------------------------------------------//
0350 /*!
0351  * Calculate the maximum scattering angle off the target nucleus.
0352  */
0353 CELER_FUNCTION real_type WentzelHelper::calc_cos_thetamax_nuclear(
0354     ParticleTrackView const& particle,
0355     MaterialView const& material,
0356     NativeCRef<WentzelOKVIData> const& wentzel) const
0357 {
0358     if (wentzel.params.is_combined)
0359     {
0360         CELER_ASSERT(material.material_id() < wentzel.inv_mass_cbrt_sq.size());
0361         return max(wentzel.params.costheta_limit,
0362                    1
0363                        - wentzel.params.a_sq_factor
0364                              * wentzel.inv_mass_cbrt_sq[material.material_id()]
0365                              / value_as<MomentumSq>(particle.momentum_sq()));
0366     }
0367     return wentzel.params.costheta_limit;
0368 }
0369
0370 //---------------------------------------------------------------------------//
0371 }  // namespace celeritas