em/distribution/WentzelDistribution.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/distribution/WentzelDistribution.hh
0006 //---------------------------------------------------------------------------//
0007 #pragma once
0008
0009 #include <cmath>
0010
0011 #include "corecel/Macros.hh"
0012 #include "corecel/Types.hh"
0013 #include "corecel/math/Algorithms.hh"
0014 #include "corecel/random/distribution/BernoulliDistribution.hh"
0015 #include "corecel/random/distribution/RejectionSampler.hh"
0016 #include "celeritas/em/interactor/detail/PhysicsConstants.hh"
0017 #include "celeritas/em/xs/MottRatioCalculator.hh"
0018 #include "celeritas/em/xs/NuclearFormFactors.hh"
0019 #include "celeritas/em/xs/WentzelHelper.hh"
0020 #include "celeritas/mat/IsotopeView.hh"
0021 #include "celeritas/phys/ParticleTrackView.hh"
0022
0023 namespace celeritas
0024 {
0025 //---------------------------------------------------------------------------//
0026 /*!
0027  * Sample the polar scattering angle cosine for Wentzel Coulomb scattering.
0028  *
0029  * This chooses between sampling scattering off an electron or nucleus based on
0030  * the relative cross sections. Electron scattering angle imposes a maximum
0031  * scattering angle (see WentzelHelper::cos_thetamax_electron), and nuclear
0032  * sattering rejects an angular change based on the Mott cross section (see
0033  * MottRatioCalculator). Nuclear scattering depends on the electronic and
0034  * nuclear cross sections (calculated by \c WentzelHelper) and nuclear form
0035  * factors (see \c ExpNuclearFormFactor, \c GaussianNuclearFormFactor, and \c
0036  * UUNuclearFormFactor ) that describe an approximate spatial charge
0037  * distribution of the nucleus. The class is used by \c
0038  * CoulombScatteringInteractor .
0039  *
0040  * The polar angle distribution is given in
0041  * \citet{fernandez-msc-1993, https://doi.org/10.1016/0168-583X(93)95827-R}
0042  * Eq. 88 and is normalized on the interval
0043  * \f$ cos\theta \in [\cos\theta_\mathrm{min}, \cos\theta_\mathrm{max}] \f$.
0044  * The sampling function for the angular deflection
0045  * \f[
0046  * \mu(\theta) \equiv \frac{1}{2}(1 - \cos\theta)
0047  * \f]
0048  * is
0049  * \f[
0050    \mu = \mu_1 + \frac{(A + \mu_1) \xi (\mu_2 - \mu_1)}{A + \mu_2 - \xi (\mu_2
0051    - \mu_1)},
0052  * \f]
0053  * where \f$ \mu_1 = \frac{1}{2}(1 - \cos\theta_\mathrm{min}) \f$,
0054  * \f$ \mu_2 = \frac{1}{2}(1 - \cos\theta_\mathrm{max}) \f$,
0055  * \f$ A \f$ is the screening coefficient, and
0056  * \f$ \xi \sim U(0,1) \f$.
0057  */
0058 class WentzelDistribution
0059 {
0060   public:
0061     //!@{
0062     //! \name Type aliases
0063     using MomentumSq = units::MevMomentumSq;
0064     using Energy = units::MevEnergy;
0065     using Mass = units::MevMass;
0066     using MottCoeffMatrix = MottElementData::MottCoeffMatrix;
0067     //!@}
0068
0069   public:
0070     // Construct with state and model data
0071     inline CELER_FUNCTION
0072     WentzelDistribution(NativeCRef<WentzelOKVIData> const& wentzel,
0073                         WentzelHelper const& helper,
0074                         ParticleTrackView const& particle,
0075                         IsotopeView const& target,
0076                         ElementId el_id,
0077                         real_type cos_thetamin,
0078                         real_type cos_thetamax);
0079
0080     // Sample the polar scattering angle
0081     template<class Engine>
0082     inline CELER_FUNCTION real_type operator()(Engine& rng) const;
0083
0084   private:
0085     //// DATA ////
0086
0087     // Shared Coulomb scattering data
0088     NativeCRef<WentzelOKVIData> const& wentzel_;
0089
0090     // Helper for calculating xs ratio and other quantities
0091     WentzelHelper const& helper_;
0092
0093     // Incident particle
0094     ParticleTrackView const& particle_;
0095
0096     // Target isotope
0097     IsotopeView const& target_;
0098
0099     // Matrix of Mott coefficients for the target element
0100     MottCoeffMatrix const& mott_coeffs_;
0101
0102     // Cosine of the minimum scattering angle
0103     real_type cos_thetamin_;
0104
0105     // Cosine of the maximum scattering angle
0106     real_type cos_thetamax_;
0107
0108     //// TYPES ////
0109
0110     using InvMomSq = RealQuantity<UnitInverse<units::MevMomentumSq::unit_type>>;
0111
0112     //// HELPER FUNCTIONS ////
0113
0114     // Calculates the form factor from the scattered polar angle
0115     inline CELER_FUNCTION real_type calculate_form_factor(real_type cos_t) const;
0116
0117     // Calculate the nuclear form momentum scale
0118     inline CELER_FUNCTION real_type nuclear_form_momentum_scale() const;
0119
0120     // Calculate the squared momentum transfer to the target
0121     inline CELER_FUNCTION real_type mom_transfer_sq(real_type cos_t) const;
0122
0123     // Momentum prefactor used in exponential and gaussian form factors
0124     inline CELER_FUNCTION InvMomSq nuclear_form_prefactor() const;
0125
0126     // Sample the scattered polar angle
0127     template<class Engine>
0128     inline CELER_FUNCTION real_type sample_costheta(real_type cos_thetamin,
0129                                                     real_type cos_thetamax,
0130                                                     Engine& rng) const;
0131
0132     // Helper function for calculating the flat form factor
0133     inline static CELER_FUNCTION real_type flat_form_factor(real_type x);
0134
0135     // Momentum coefficient used in the flat nuclear form factor model
0136     static CELER_CONSTEXPR_FUNCTION real_type flat_coeff();
0137 };
0138
0139 //---------------------------------------------------------------------------//
0140 // INLINE DEFINITIONS
0141 //---------------------------------------------------------------------------//
0142 /*!
0143  * Construct with state and model data.
0144  */
0145 CELER_FUNCTION
0146 WentzelDistribution::WentzelDistribution(
0147     NativeCRef<WentzelOKVIData> const& wentzel,
0148     WentzelHelper const& helper,
0149     ParticleTrackView const& particle,
0150     IsotopeView const& target,
0151     ElementId el_id,
0152     real_type cos_thetamin,
0153     real_type cos_thetamax)
0154     : wentzel_(wentzel)
0155     , helper_(helper)
0156     , particle_(particle)
0157     , target_(target)
0158     , mott_coeffs_(particle_.charge() < zero_quantity()
0159                        ? wentzel_.mott_coeffs[el_id].electron
0160                        : wentzel_.mott_coeffs[el_id].positron)
0161     , cos_thetamin_(cos_thetamin)
0162     , cos_thetamax_(cos_thetamax)
0163 {
0164     CELER_EXPECT(el_id < wentzel_.mott_coeffs.size());
0165     CELER_EXPECT(cos_thetamin_ >= -1 && cos_thetamin_ <= 1);
0166     CELER_EXPECT(cos_thetamax_ >= -1 && cos_thetamax_ <= 1);
0167     CELER_EXPECT(cos_thetamax_ <= cos_thetamin_);
0168 }
0169
0170 //---------------------------------------------------------------------------//
0171 /*!
0172  * Sample the cosine polar scattered angle of the incident particle.
0173  */
0174 template<class Engine>
0175 CELER_FUNCTION real_type WentzelDistribution::operator()(Engine& rng) const
0176 {
0177     real_type cos_theta = 1;
0178     if (BernoulliDistribution(
0179             helper_.calc_xs_electron(cos_thetamin_, cos_thetamax_),
0180             helper_.calc_xs_nuclear(cos_thetamin_, cos_thetamax_))(rng))
0181     {
0182         // Scattered off of electrons
0183         real_type const cos_thetamax_elec = helper_.cos_thetamax_electron();
0184         real_type cos_thetamin = max(cos_thetamin_, cos_thetamax_elec);
0185         real_type cos_thetamax = max(cos_thetamax_, cos_thetamax_elec);
0186         CELER_ASSERT(cos_thetamin > cos_thetamax);
0187         cos_theta = this->sample_costheta(cos_thetamin, cos_thetamax, rng);
0188     }
0189     else
0190     {
0191         // Scattered off of nucleus
0192         cos_theta = this->sample_costheta(cos_thetamin_, cos_thetamax_, rng);
0193
0194         // Calculate rejection for false scattering
0195         MottRatioCalculator calc_mott_ratio(mott_coeffs_,
0196                                             std::sqrt(particle_.beta_sq()));
0197         real_type xs = calc_mott_ratio(cos_theta)
0198                        * ipow<2>(this->calculate_form_factor(cos_theta));
0199         if (RejectionSampler(xs, helper_.mott_factor())(rng))
0200         {
0201             // Reject scattering event: no change in direction
0202             cos_theta = 1;
0203         }
0204     }
0205     return cos_theta;
0206 }
0207
0208 //---------------------------------------------------------------------------//
0209 /*!
0210  * Calculate the form factor based on the form factor model.
0211  *
0212  * The models, described in \citet{leroy-2016,
0213  * https://doi-org.ezproxy.cern.ch/10.1142/9167} section 2.4.2.1, parameterize
0214  * the charge distribution inside a nucleus.
0215  */
0216 CELER_FUNCTION real_type
0217 WentzelDistribution::calculate_form_factor(real_type cos_t) const
0218 {
0219     CELER_EXPECT(cos_t >= -1 && cos_t <= 1);
0220
0221     using MomSq = NuclearFormFactorTraits::MomentumSq;
0222
0223     if (wentzel_.params.form_factor_type == NuclearFormFactorType::none)
0224     {
0225         return 1;
0226     }
0227
0228     // Approximate squared momentum transfer to the target (assuming heavy
0229     // target and light projectile)
0230     auto mt_sq
0231         = MomSq{2 * value_as<MomSq>(particle_.momentum_sq()) * (1 - cos_t)};
0232
0233     auto calc_ff = [mt_sq](auto&& calc_form_factor) {
0234         real_type result = calc_form_factor(mt_sq);
0235         CELER_ENSURE(result >= 0 && result <= 1);
0236         return result;
0237     };
0238
0239     switch (wentzel_.params.form_factor_type)
0240     {
0241         case NuclearFormFactorType::flat:
0242             return calc_ff(UUNuclearFormFactor{target_.atomic_mass_number()});
0243         case NuclearFormFactorType::exponential:
0244             return calc_ff(
0245                 ExpNuclearFormFactor{this->nuclear_form_prefactor()});
0246         case NuclearFormFactorType::gaussian:
0247             return calc_ff(
0248                 GaussianNuclearFormFactor{this->nuclear_form_prefactor()});
0249         default:
0250             CELER_ASSERT_UNREACHABLE();
0251     }
0252 }
0253
0254 //---------------------------------------------------------------------------//
0255 /*!
0256  * Calculate the flat form factor.
0257  *
0258  * See \cite{leroy-2016} Eq. 2.265.
0259  */
0260 CELER_FUNCTION real_type WentzelDistribution::flat_form_factor(real_type x)
0261 {
0262     return 3 * (std::sin(x) - x * std::cos(x)) / ipow<3>(x);
0263 }
0264
0265 //---------------------------------------------------------------------------//
0266 /*!
0267  * Momentum coefficient used in the flat model for the nuclear form factors.
0268  *
0269  * This is the ratio of \f$ r_1 / \hbar \f$ where \f$ r_1 \f$ is defined in
0270  * Eq. 2.265 of \cite{leroy-2016}.
0271  */
0272 CELER_CONSTEXPR_FUNCTION real_type WentzelDistribution::flat_coeff()
0273 {
0274     return native_value_to<units::MevMomentum>(2 * units::femtometer
0275                                                / constants::hbar_planck)
0276         .value();
0277 }
0278
0279 //---------------------------------------------------------------------------//
0280 /*!
0281  * Get the constant prefactor of the squared momentum transfer.
0282  */
0283 CELER_FUNCTION auto
0284 WentzelDistribution::nuclear_form_prefactor() const -> InvMomSq
0285 {
0286     CELER_EXPECT(target_.isotope_id() < wentzel_.nuclear_form_prefactor.size());
0287     return InvMomSq{wentzel_.nuclear_form_prefactor[target_.isotope_id()]};
0288 }
0289
0290 //---------------------------------------------------------------------------//
0291 /*!
0292  * Sample the scattering polar angle of the incident particle.
0293  */
0294 template<class Engine>
0295 CELER_FUNCTION real_type WentzelDistribution::sample_costheta(
0296     real_type cos_thetamin, real_type cos_thetamax, Engine& rng) const
0297 {
0298     // Sample scattering angle [fern] eqn 92, where cos(theta) = 1 - 2*mu
0299     real_type const mu1 = real_type{0.5} * (1 - cos_thetamin);
0300     real_type const mu2 = real_type{0.5} * (1 - cos_thetamax);
0301     real_type const w = generate_canonical(rng) * (mu2 - mu1);
0302     real_type const sc = helper_.screening_coefficient();
0303
0304     real_type result = 1 - 2 * mu1 - 2 * (sc + mu1) * w / (sc + mu2 - w);
0305     CELER_ENSURE(result >= -1 && result <= 1);
0306     return result;
0307 }
0308
0309 //---------------------------------------------------------------------------//
0310 }  // namespace celeritas