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