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