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