em/distribution/WentzelDistribution.hh

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