em/interactor/CoulombScatteringInteractor.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/interactor/CoulombScatteringInteractor.hh
0007 //---------------------------------------------------------------------------//
0008 #pragma once
0009
0010 #include "corecel/Macros.hh"
0011 #include "corecel/math/ArrayUtils.hh"
0012 #include "celeritas/Constants.hh"
0013 #include "celeritas/Quantities.hh"
0014 #include "celeritas/Types.hh"
0015 #include "celeritas/em/data/CoulombScatteringData.hh"
0016 #include "celeritas/em/data/WentzelOKVIData.hh"
0017 #include "celeritas/em/distribution/WentzelDistribution.hh"
0018 #include "celeritas/em/xs/WentzelHelper.hh"
0019 #include "celeritas/mat/ElementView.hh"
0020 #include "celeritas/mat/MaterialView.hh"
0021 #include "celeritas/phys/CutoffView.hh"
0022 #include "celeritas/phys/Interaction.hh"
0023 #include "celeritas/phys/InteractionUtils.hh"
0024 #include "celeritas/phys/ParticleTrackView.hh"
0025
0026 #include "detail/PhysicsConstants.hh"
0027
0028 namespace celeritas
0029 {
0030 //---------------------------------------------------------------------------//
0031 /*!
0032  * Applies the Wentzel single Coulomb scattering model.
0033  *
0034  * This models incident high-energy electrons and positrons elastically
0035  * scattering off of nuclei and atomic electrons. Scattering off of the nucleus
0036  * versus electrons is randomly sampled based on the relative cross-sections.
0037  * No secondaries are created in this process (in the future, with hadronic
0038  * transport support, secondary ions may be emitted), however production cuts
0039  * are used to determine the maximum scattering angle off of electrons.
0040  *
0041  * \note This performs the same sampling as in Geant4's
0042  *  G4eCoulombScatteringModel, as documented in section 8.2 of the Geant4
0043  *  Physics Reference Manual (release 11.1).
0044  */
0045 class CoulombScatteringInteractor
0046 {
0047   public:
0048     //!@{
0049     //! \name Type aliases
0050     using Energy = units::MevEnergy;
0051     using Mass = units::MevMass;
0052     using MomentumSq = units::MevMomentumSq;
0053     //!@}
0054
0055   public:
0056     //! Construct with shared and state data
0057     inline CELER_FUNCTION
0058     CoulombScatteringInteractor(CoulombScatteringData const& shared,
0059                                 NativeCRef<WentzelOKVIData> const& wentzel,
0060                                 ParticleTrackView const& particle,
0061                                 Real3 const& inc_direction,
0062                                 MaterialView const& material,
0063                                 IsotopeView const& target,
0064                                 ElementId el_id,
0065                                 CutoffView const& cutoffs);
0066
0067     //! Sample an interaction with the given RNG
0068     template<class Engine>
0069     inline CELER_FUNCTION Interaction operator()(Engine& rng);
0070
0071   private:
0072     //// DATA ////
0073
0074     // Incident direction
0075     Real3 const& inc_direction_;
0076
0077     // Incident particle
0078     ParticleTrackView const& particle_;
0079
0080     // Target isotope
0081     IsotopeView const& target_;
0082
0083     // Helper for calculating xs ratio and other quantities
0084     WentzelHelper const helper_;
0085
0086     // Scattering direction distribution of the Wentzel model
0087     WentzelDistribution const sample_angle_;
0088
0089     //// HELPER FUNCTIONS ////
0090
0091     // Calculates the recoil energy for the given scattering angle
0092     inline CELER_FUNCTION real_type calc_recoil_energy(real_type cos_theta) const;
0093 };
0094
0095 //---------------------------------------------------------------------------//
0096 // INLINE DEFINITIONS
0097 //---------------------------------------------------------------------------//
0098 /*!
0099  * Construct from shared and state data.
0100  *
0101  * \todo Use the proton production cutoff when the recoiled nucleus production
0102  * is supported.
0103  */
0104 CELER_FUNCTION
0105 CoulombScatteringInteractor::CoulombScatteringInteractor(
0106     CoulombScatteringData const& shared,
0107     NativeCRef<WentzelOKVIData> const& wentzel,
0108     ParticleTrackView const& particle,
0109     Real3 const& inc_direction,
0110     MaterialView const& material,
0111     IsotopeView const& target,
0112     ElementId el_id,
0113     CutoffView const& cutoffs)
0114     : inc_direction_(inc_direction)
0115     , particle_(particle)
0116     , target_(target)
0117     , helper_(particle,
0118               material,
0119               target.atomic_number(),
0120               wentzel,
0121               shared.ids,
0122               cutoffs.energy(shared.ids.electron))
0123     , sample_angle_(wentzel,
0124                     helper_,
0125                     particle,
0126                     target,
0127                     el_id,
0128                     helper_.cos_thetamax_nuclear(),
0129                     shared.cos_thetamax())
0130 {
0131     CELER_EXPECT(particle_.particle_id() == shared.ids.electron
0132                  || particle_.particle_id() == shared.ids.positron);
0133     CELER_EXPECT(particle_.energy() > zero_quantity()
0134                  && particle_.energy() < detail::high_energy_limit());
0135 }
0136
0137 //---------------------------------------------------------------------------//
0138 /*!
0139  * Sample the Coulomb scattering of the incident particle.
0140  */
0141 template<class Engine>
0142 CELER_FUNCTION Interaction CoulombScatteringInteractor::operator()(Engine& rng)
0143 {
0144     // Incident particle scatters
0145     Interaction result;
0146
0147     // Sample the new direction
0148     real_type cos_theta = sample_angle_(rng);
0149     result.direction = ExitingDirectionSampler{cos_theta, inc_direction_}(rng);
0150
0151     // Recoil energy is kinetic energy transfered to the atom
0152     real_type inc_energy = value_as<Energy>(particle_.energy());
0153     real_type recoil_energy = this->calc_recoil_energy(cos_theta);
0154     CELER_ASSERT(0 <= recoil_energy && recoil_energy <= inc_energy);
0155     result.energy = Energy{inc_energy - recoil_energy};
0156
0157     //! \todo For high enough recoil energies, ions are produced
0158     result.energy_deposition = Energy{recoil_energy};
0159
0160     return result;
0161 }
0162
0163 //---------------------------------------------------------------------------//
0164 /*!
0165  * Calculates the recoil energy for the given scattering angle sampled
0166  * by WentzelDistribution.
0167  */
0168 CELER_FUNCTION real_type
0169 CoulombScatteringInteractor::calc_recoil_energy(real_type cos_theta) const
0170 {
0171     real_type projectile_momsq = value_as<MomentumSq>(particle_.momentum_sq());
0172     real_type projectile_mass = value_as<Mass>(particle_.mass())
0173                                 + value_as<Energy>(particle_.energy());
0174     real_type target_mass = value_as<Mass>(target_.nuclear_mass());
0175
0176     return projectile_momsq * (1 - cos_theta)
0177            / (target_mass + projectile_mass * (1 - cos_theta));
0178 }
0179
0180 //---------------------------------------------------------------------------//
0181 }  // namespace celeritas