em/interactor/KleinNishinaInteractor.hh

0001 //----------------------------------*-C++-*----------------------------------//
0002 // Copyright 2021-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/KleinNishinaInteractor.hh
0007 //---------------------------------------------------------------------------//
0008 #pragma once
0009
0010 #include "corecel/Macros.hh"
0011 #include "corecel/Types.hh"
0012 #include "corecel/data/StackAllocator.hh"
0013 #include "corecel/math/ArrayUtils.hh"
0014 #include "celeritas/Constants.hh"
0015 #include "celeritas/Quantities.hh"
0016 #include "celeritas/em/data/KleinNishinaData.hh"
0017 #include "celeritas/phys/Interaction.hh"
0018 #include "celeritas/phys/InteractionUtils.hh"
0019 #include "celeritas/phys/ParticleTrackView.hh"
0020 #include "celeritas/phys/Secondary.hh"
0021 #include "celeritas/random/distribution/BernoulliDistribution.hh"
0022 #include "celeritas/random/distribution/ReciprocalDistribution.hh"
0023 #include "celeritas/random/distribution/UniformRealDistribution.hh"
0024
0025 namespace celeritas
0026 {
0027 //---------------------------------------------------------------------------//
0028 /*!
0029  * Perform Compton scattering, neglecting atomic binding energy.
0030  *
0031  * This is a model for the discrete Compton inelastic scattering process. Given
0032  * an incident gamma, it adds a single secondary (electron) to the secondary
0033  * stack and returns an interaction for the change to the incident gamma
0034  * direction and energy. No cutoffs are performed for the incident energy or
0035  * the exiting gamma energy. A secondary production cutoff is applied to the
0036  * outgoing electron.
0037  *
0038  * \note This performs the same sampling routine as in Geant4's
0039  *  G4KleinNishinaCompton, as documented in section 6.4.2 of the Geant4 Physics
0040  *  Reference (release 10.6).
0041  */
0042 class KleinNishinaInteractor
0043 {
0044   public:
0045     // Construct from shared and state data
0046     inline CELER_FUNCTION
0047     KleinNishinaInteractor(KleinNishinaData const& shared,
0048                            ParticleTrackView const& particle,
0049                            Real3 const& inc_direction,
0050                            StackAllocator<Secondary>& allocate);
0051
0052     // Sample an interaction with the given RNG
0053     template<class Engine>
0054     inline CELER_FUNCTION Interaction operator()(Engine& rng);
0055
0056     //! Energy threshold for secondary production [MeV]
0057     static CELER_CONSTEXPR_FUNCTION units::MevEnergy secondary_cutoff()
0058     {
0059         return units::MevEnergy{1e-4};
0060     }
0061
0062   private:
0063     // Constant data
0064     KleinNishinaData const& shared_;
0065     // Incident gamma energy
0066     units::MevEnergy const inc_energy_;
0067     // Incident direction
0068     Real3 const& inc_direction_;
0069     // Allocate space for a secondary particle
0070     StackAllocator<Secondary>& allocate_;
0071 };
0072
0073 //---------------------------------------------------------------------------//
0074 // INLINE DEFINITIONS
0075 //---------------------------------------------------------------------------//
0076 /*!
0077  * Construct with shared and state data.
0078  *
0079  * The incident particle must be above the energy threshold: this should be
0080  * handled in code *before* the interactor is constructed.
0081  */
0082 CELER_FUNCTION KleinNishinaInteractor::KleinNishinaInteractor(
0083     KleinNishinaData const& shared,
0084     ParticleTrackView const& particle,
0085     Real3 const& inc_direction,
0086     StackAllocator<Secondary>& allocate)
0087     : shared_(shared)
0088     , inc_energy_(particle.energy())
0089     , inc_direction_(inc_direction)
0090     , allocate_(allocate)
0091 {
0092     CELER_EXPECT(particle.particle_id() == shared_.ids.gamma);
0093     CELER_EXPECT(inc_energy_ > zero_quantity());
0094 }
0095
0096 //---------------------------------------------------------------------------//
0097 /*!
0098  * Sample Compton scattering using the Klein-Nishina model.
0099  *
0100  * See section 6.4.2 of the Geant physics reference. Epsilon is the ratio of
0101  * outgoing to incident gamma energy, bounded in [epsilon_0, 1].
0102  */
0103 template<class Engine>
0104 CELER_FUNCTION Interaction KleinNishinaInteractor::operator()(Engine& rng)
0105 {
0106     using Energy = units::MevEnergy;
0107
0108     // Allocate space for the single electron to be emitted
0109     Secondary* electron_secondary = allocate_(1);
0110     if (electron_secondary == nullptr)
0111     {
0112         // Failed to allocate space for a secondary
0113         return Interaction::from_failure();
0114     }
0115
0116     // Value of epsilon corresponding to minimum photon energy
0117     real_type const inc_energy_per_mecsq = value_as<Energy>(inc_energy_)
0118                                            * shared_.inv_electron_mass;
0119     real_type const epsilon_0 = 1 / (1 + 2 * inc_energy_per_mecsq);
0120
0121     // Probability of alpha_1 to choose f_1 (sample epsilon)
0122     BernoulliDistribution choose_f1(-std::log(epsilon_0),
0123                                     real_type(0.5) * (1 - ipow<2>(epsilon_0)));
0124     // Sample f_1(\eps) \propto 1/\eps on [\eps_0, 1]
0125     ReciprocalDistribution<real_type> sample_f1(epsilon_0);
0126     // Sample square of f_2(\eps^2) \propto 1 on [\eps_0^2, 1]
0127     UniformRealDistribution<real_type> sample_f2_sq(ipow<2>(epsilon_0), 1);
0128
0129     // Rejection loop: sample epsilon (energy change) and direction change
0130     real_type epsilon;
0131     real_type one_minus_costheta;
0132     // Temporary sample values used in rejection
0133     real_type reject_prob;
0134     do
0135     {
0136         // Sample epsilon and square
0137         real_type epsilon_sq;
0138         if (choose_f1(rng))
0139         {
0140             // Sample f_1(\eps) \propto 1/\eps on [\eps_0, 1]
0141             // => \eps \gets \eps_0^\xi = \exp(\xi \log \eps_0)
0142             epsilon = sample_f1(rng);
0143             epsilon_sq = epsilon * epsilon;
0144         }
0145         else
0146         {
0147             // Sample f_2(\eps) = 2 * \eps / (1 - epsilon_0 * epsilon_0)
0148             epsilon_sq = sample_f2_sq(rng);
0149             epsilon = std::sqrt(epsilon_sq);
0150         }
0151         CELER_ASSERT(epsilon >= epsilon_0 && epsilon <= 1);
0152
0153         // Calculate angles: need sin^2 \theta for rejection
0154         one_minus_costheta = (1 - epsilon) / (epsilon * inc_energy_per_mecsq);
0155         CELER_ASSERT(one_minus_costheta >= 0 && one_minus_costheta <= 2);
0156         real_type sintheta_sq = one_minus_costheta * (2 - one_minus_costheta);
0157         reject_prob = epsilon * sintheta_sq / (1 + epsilon_sq);
0158     } while (BernoulliDistribution(reject_prob)(rng));
0159
0160     // Construct interaction for change to primary (incident) particle
0161     Interaction result;
0162     result.energy = Energy{epsilon * inc_energy_.value()};
0163     result.direction = inc_direction_;
0164     result.secondaries = {electron_secondary, 1};
0165
0166     // Sample azimuthal direction and rotate the outgoing direction
0167     result.direction = ExitingDirectionSampler{1 - one_minus_costheta,
0168                                                result.direction}(rng);
0169
0170     // Construct secondary energy by neglecting electron binding energy
0171     electron_secondary->energy
0172         = Energy{inc_energy_.value() - result.energy.value()};
0173
0174     // Apply secondary production cutoff
0175     if (electron_secondary->energy < KleinNishinaInteractor::secondary_cutoff())
0176     {
0177         result.energy_deposition = electron_secondary->energy;
0178         *electron_secondary = {};
0179         return result;
0180     }
0181
0182     // Outgoing secondary is an electron
0183     electron_secondary->particle_id = shared_.ids.electron;
0184     // Calculate exiting electron direction via conservation of momentum
0185     electron_secondary->direction
0186         = calc_exiting_direction({inc_energy_.value(), inc_direction_},
0187                                  {result.energy.value(), result.direction});
0188
0189     return result;
0190 }
0191
0192 //---------------------------------------------------------------------------//
0193 }  // namespace celeritas