em/interactor/EPlusGGInteractor.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/EPlusGGInteractor.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/ArrayOperators.hh"
0013 #include "corecel/math/ArrayUtils.hh"
0014 #include "corecel/random/distribution/BernoulliDistribution.hh"
0015 #include "corecel/random/distribution/ReciprocalDistribution.hh"
0016 #include "geocel/random/IsotropicDistribution.hh"
0017 #include "celeritas/Quantities.hh"
0018 #include "celeritas/em/data/EPlusGGData.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  * Annihilate a positron to create two gammas.
0029  *
0030  * This is a model for the discrete positron-electron annihilation process
0031  * which simulates the in-flight annihilation of a positron with an atomic
0032  * electron and produces into two photons. It is assumed that the atomic
0033  * electron is initially free and at rest.
0034  *
0035  * \note This performs the same sampling routine as in Geant4's
0036  * G4eeToTwoGammaModel class, as documented in section 10.3 of the Geant4
0037  * Physics Reference (release 10.6). The maximum energy for the model is
0038  * specified in Geant4.
0039  */
0040 class EPlusGGInteractor
0041 {
0042   public:
0043     //!@{
0044     //! \name Type aliases
0045     using Mass = units::MevMass;
0046     using Energy = units::MevEnergy;
0047     //!@}
0048
0049   public:
0050     // Construct with shared and state data
0051     inline CELER_FUNCTION
0052     EPlusGGInteractor(EPlusGGData const& shared,
0053                       ParticleTrackView const& particle,
0054                       Real3 const& inc_direction,
0055                       StackAllocator<Secondary>& allocate);
0056
0057     // Sample an interaction with the given RNG
0058     template<class Engine>
0059     inline CELER_FUNCTION Interaction operator()(Engine& rng);
0060
0061   private:
0062     // Shared constant physics properties
0063     EPlusGGData const& shared_;
0064     // Incident positron energy [MevEnergy]
0065     real_type const inc_energy_;
0066     // Incident direction
0067     Real3 const& inc_direction_;
0068     // Allocate space for secondary particles (two photons)
0069     StackAllocator<Secondary>& allocate_;
0070 };
0071
0072 //---------------------------------------------------------------------------//
0073 // INLINE DEFINITIONS
0074 //---------------------------------------------------------------------------//
0075 /*!
0076  * Construct with shared and state data.
0077  */
0078 CELER_FUNCTION
0079 EPlusGGInteractor::EPlusGGInteractor(EPlusGGData const& shared,
0080                                      ParticleTrackView const& particle,
0081                                      Real3 const& inc_direction,
0082                                      StackAllocator<Secondary>& allocate)
0083     : shared_(shared)
0084     , inc_energy_(value_as<Energy>(particle.energy()))
0085     , inc_direction_(inc_direction)
0086     , allocate_(allocate)
0087 {
0088     CELER_EXPECT(particle.particle_id() == shared_.positron);
0089 }
0090
0091 //---------------------------------------------------------------------------//
0092 /*!
0093  * Sample two gammas using the G4eeToTwoGammaModel model.
0094  *
0095  * Polarization is *not* implemented.
0096  */
0097 template<class Engine>
0098 CELER_FUNCTION Interaction EPlusGGInteractor::operator()(Engine& rng)
0099 {
0100     // Allocate space for two gammas
0101     Secondary* secondaries = allocate_(2);
0102     if (secondaries == nullptr)
0103     {
0104         // Failed to allocate space for two secondaries
0105         return Interaction::from_failure();
0106     }
0107
0108     // Construct an interaction with an absorbed process
0109     Interaction result = Interaction::from_absorption();
0110     result.secondaries = {secondaries, 2};
0111
0112     // Sample two gammas
0113     secondaries[0].particle_id = secondaries[1].particle_id = shared_.gamma;
0114
0115     if (inc_energy_ == 0)
0116     {
0117         // Save outgoing secondary data
0118         secondaries[0].energy = secondaries[1].energy
0119             = Energy{value_as<Mass>(shared_.electron_mass)};
0120
0121         IsotropicDistribution<real_type> gamma_dir;
0122         secondaries[0].direction = gamma_dir(rng);
0123         secondaries[1].direction = -secondaries[0].direction;
0124     }
0125     else
0126     {
0127         constexpr real_type half = 0.5;
0128         real_type const tau = inc_energy_
0129                               / value_as<Mass>(shared_.electron_mass);
0130         real_type const tau2 = tau + 2;
0131         real_type const sqgrate = std::sqrt(tau / tau2) * half;
0132
0133         // Evaluate limits of the energy sampling
0134         ReciprocalDistribution<real_type> sample_eps(half - sqgrate,
0135                                                      half + sqgrate);
0136
0137         // Sample the energy rate of the created gammas
0138         real_type epsil;
0139         do
0140         {
0141             epsil = sample_eps(rng);
0142         } while (BernoulliDistribution(
0143             epsil - (2 * (tau + 1) * epsil - 1) / (epsil * ipow<2>(tau2)))(rng));
0144
0145         // Scattered Gamma angles
0146         real_type const cost = (epsil * tau2 - 1)
0147                                / (epsil * std::sqrt(tau * tau2));
0148         CELER_ASSERT(std::fabs(cost) <= 1);
0149
0150         // Kinematic of the gamma pair
0151         real_type const total_energy
0152             = inc_energy_ + 2 * value_as<Mass>(shared_.electron_mass);
0153         real_type const gamma_energy = epsil * total_energy;
0154         real_type const eplus_moment = std::sqrt(inc_energy_ * total_energy);
0155
0156         // Sample and save outgoing secondary data
0157         secondaries[0].energy = Energy{gamma_energy};
0158         secondaries[0].direction
0159             = ExitingDirectionSampler{cost, inc_direction_}(rng);
0160         secondaries[1].energy = Energy{total_energy - gamma_energy};
0161         secondaries[1].direction = calc_exiting_direction(
0162             {eplus_moment, inc_direction_}, {inc_energy_, inc_direction_});
0163     }
0164
0165     return result;
0166 }
0167
0168 //---------------------------------------------------------------------------//
0169 }  // namespace celeritas