em/interactor/BetheHeitlerInteractor.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/BetheHeitlerInteractor.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/Algorithms.hh"
0014 #include "corecel/math/ArrayUtils.hh"
0015 #include "celeritas/Constants.hh"
0016 #include "celeritas/Quantities.hh"
0017 #include "celeritas/em/data/BetheHeitlerData.hh"
0018 #include "celeritas/em/distribution/TsaiUrbanDistribution.hh"
0019 #include "celeritas/em/xs/LPMCalculator.hh"
0020 #include "celeritas/mat/ElementView.hh"
0021 #include "celeritas/phys/Interaction.hh"
0022 #include "celeritas/phys/ParticleTrackView.hh"
0023 #include "celeritas/phys/Secondary.hh"
0024 #include "celeritas/random/distribution/BernoulliDistribution.hh"
0025 #include "celeritas/random/distribution/GenerateCanonical.hh"
0026 #include "celeritas/random/distribution/UniformRealDistribution.hh"
0027
0028 namespace celeritas
0029 {
0030 //---------------------------------------------------------------------------//
0031 /*!
0032  * Relativistic model for electron-positron pair production.
0033  *
0034  * The energies of the secondary electron and positron are sampled using the
0035  * Bethe-Heitler cross sections with a Coulomb correction. The LPM effect is
0036  * taken into account for incident gamma energies above 100 GeV. Exiting
0037  * particle directions are sampled with the \c TsaiUrbanDistribution . Note
0038  * that energy is not exactly conserved.
0039  *
0040  * \note This performs the same sampling routine as in Geant4's
0041  * G4PairProductionRelModel, as documented in sections 6.5 (gamma conversion)
0042  * and 10.2.2 (LPM effect) of the Geant4 Physics Reference Manual (release
0043  * 10.7)
0044  */
0045 class BetheHeitlerInteractor
0046 {
0047   public:
0048     //!@{
0049     //! \name Type aliases
0050     using Mass = units::MevMass;
0051     using Energy = units::MevEnergy;
0052     //!@}
0053
0054   public:
0055     //! Construct sampler from shared and state data
0056     inline CELER_FUNCTION
0057     BetheHeitlerInteractor(BetheHeitlerData const& shared,
0058                            ParticleTrackView const& particle,
0059                            Real3 const& inc_direction,
0060                            StackAllocator<Secondary>& allocate,
0061                            MaterialView const& material,
0062                            ElementView const& element);
0063
0064     // Sample an interaction with the given RNG
0065     template<class Engine>
0066     inline CELER_FUNCTION Interaction operator()(Engine& rng);
0067
0068   private:
0069     //// TYPES ////
0070
0071     //! Screening functions \f$ \Phi_1 \f$ and \f$ \Phi_2 \f$
0072     struct ScreeningFunctions
0073     {
0074         real_type phi1;
0075         real_type phi2;
0076     };
0077
0078     //// DATA ////
0079
0080     // Shared model data
0081     BetheHeitlerData const& shared_;
0082     // Incident gamma energy
0083     real_type const inc_energy_;
0084     // Incident direction
0085     Real3 const& inc_direction_;
0086     // Allocate space for a secondary particle
0087     StackAllocator<Secondary>& allocate_;
0088     // Element properties for calculating screening functions and variables
0089     ElementView const& element_;
0090     // Whether LPM supression is applied
0091     bool const enable_lpm_;
0092     // Used to calculate the LPM suppression functions
0093     LPMCalculator calc_lpm_functions_;
0094     // Cached minimum epsilon, m_e*c^2/E_gamma; kinematical limit for Y -> e+e-
0095     real_type epsilon0_;
0096
0097     //// CONSTANTS ////
0098
0099     //! Energy above which the Coulomb correction is applied [MeV]
0100     static CELER_CONSTEXPR_FUNCTION Energy coulomb_corr_threshold()
0101     {
0102         return units::MevEnergy{50};
0103     }
0104
0105     //! Energy above which LPM suppression is applied, if enabled [MeV]
0106     static CELER_CONSTEXPR_FUNCTION Energy lpm_threshold()
0107     {
0108         return units::MevEnergy{1e5};
0109     }
0110
0111     //// HELPER FUNCTIONS ////
0112
0113     // Calculate the screening variable \f$ \delta \f$
0114     inline CELER_FUNCTION real_type impact_parameter(real_type eps) const;
0115
0116     // Calculate the screening functions \f$ \Phi_1 \f$ and \f$ \Phi_2 \f$
0117     inline CELER_FUNCTION ScreeningFunctions
0118     screening_phi1_phi2(real_type delta) const;
0119
0120     // Calculate the auxiliary screening function \f$ F_1 \f$
0121     inline CELER_FUNCTION real_type screening_f1(real_type delta) const;
0122
0123     // Calculate the auxiliary screening function \f$ F_2 \f$
0124     inline CELER_FUNCTION real_type screening_f2(real_type delta) const;
0125 };
0126
0127 //---------------------------------------------------------------------------//
0128 // INLINE DEFINITIONS
0129 //---------------------------------------------------------------------------//
0130 /*!
0131  * Construct with shared and state data.
0132  *
0133  * The incident gamma energy must be at least twice the electron rest mass.
0134  */
0135 CELER_FUNCTION BetheHeitlerInteractor::BetheHeitlerInteractor(
0136     BetheHeitlerData const& shared,
0137     ParticleTrackView const& particle,
0138     Real3 const& inc_direction,
0139     StackAllocator<Secondary>& allocate,
0140     MaterialView const& material,
0141     ElementView const& element)
0142     : shared_(shared)
0143     , inc_energy_(value_as<Energy>(particle.energy()))
0144     , inc_direction_(inc_direction)
0145     , allocate_(allocate)
0146     , element_(element)
0147     , enable_lpm_(shared.enable_lpm
0148                   && inc_energy_ > value_as<Energy>(lpm_threshold()))
0149     , calc_lpm_functions_(material,
0150                           element_,
0151                           shared_.dielectric_suppression(),
0152                           Energy{inc_energy_})
0153 {
0154     CELER_EXPECT(particle.particle_id() == shared_.ids.gamma);
0155     CELER_EXPECT(inc_energy_ >= 2 * value_as<Mass>(shared_.electron_mass));
0156
0157     epsilon0_ = value_as<Mass>(shared_.electron_mass) / inc_energy_;
0158 }
0159
0160 //---------------------------------------------------------------------------//
0161 /*!
0162  * Sample the distribution.
0163  */
0164 template<class Engine>
0165 CELER_FUNCTION Interaction BetheHeitlerInteractor::operator()(Engine& rng)
0166 {
0167     // Allocate space for the electron/positron pair
0168     Secondary* secondaries = allocate_(2);
0169     if (secondaries == nullptr)
0170     {
0171         // Failed to allocate space for secondaries
0172         return Interaction::from_failure();
0173     }
0174
0175     constexpr real_type half = 0.5;
0176
0177     // If E_gamma < 2 MeV, rejection not needed -- just sample uniformly
0178     real_type epsilon;
0179     if (inc_energy_ < value_as<Energy>(units::MevEnergy{2}))
0180     {
0181         UniformRealDistribution<real_type> sample_eps(epsilon0_, half);
0182         epsilon = sample_eps(rng);
0183     }
0184     else
0185     {
0186         // Calculate the minimum (when \epsilon = 1/2) and maximum (when
0187         // \epsilon = \epsilon_1) values of screening variable, \delta. Above
0188         // 50 MeV, a Coulomb correction function is introduced.
0189         real_type const delta_min = 4 * 136 / element_.cbrt_z() * epsilon0_;
0190         real_type f_z = real_type(8) / real_type(3) * element_.log_z();
0191         if (inc_energy_ > value_as<Energy>(coulomb_corr_threshold()))
0192         {
0193             f_z += 8 * element_.coulomb_correction();
0194         }
0195         real_type const delta_max
0196             = std::exp((real_type(42.038) - f_z) / real_type(8.29))
0197               - real_type(0.958);
0198         CELER_ASSERT(delta_min <= delta_max);
0199
0200         // Calculate the lower limit of epsilon. Due to the Coulomb correction,
0201         // the cross section can become negative even at kinematically allowed
0202         // \epsilon > \epsilon_0 values. To exclude these negative cross
0203         // sections, an additional constraint that \epsilon > \epsilon_1 is
0204         // introduced, where \epsilon_1 is the solution to
0205         // \Phi(\delta(\epsilon)) - F(Z)/2 = 0.
0206         real_type const epsilon1
0207             = half - half * std::sqrt(1 - delta_min / delta_max);
0208         real_type const epsilon_min = celeritas::max(epsilon0_, epsilon1);
0209
0210         // Decide to choose f1, g1 or f2, g2 based on N1, N2 (factors from
0211         // corrected Bethe-Heitler cross section; c.f. Eq. 6.6 of Geant4
0212         // Physics Reference 10.6)
0213         real_type const f10 = this->screening_f1(delta_min) - f_z;
0214         real_type const f20 = this->screening_f2(delta_min) - f_z;
0215         BernoulliDistribution choose_f1g1(ipow<2>(half - epsilon_min) * f10,
0216                                           real_type(1.5) * f20);
0217
0218         // Rejection function g_1 or g_2. Note the it's possible for g to be
0219         // greater than one
0220         real_type g;
0221         do
0222         {
0223             if (choose_f1g1(rng))
0224             {
0225                 // Used to sample from f1
0226                 epsilon = half
0227                           - (half - epsilon_min)
0228                                 * std::cbrt(generate_canonical(rng));
0229                 CELER_ASSERT(epsilon >= epsilon_min && epsilon <= half);
0230
0231                 // Calculate delta from element atomic number and sampled
0232                 // epsilon
0233                 real_type delta = this->impact_parameter(epsilon);
0234                 CELER_ASSERT(delta <= delta_max && delta >= delta_min);
0235
0236                 // Calculate g_1 rejection function
0237                 if (enable_lpm_)
0238                 {
0239                     auto screening = screening_phi1_phi2(delta);
0240                     auto lpm = calc_lpm_functions_(epsilon);
0241                     g = lpm.xi
0242                         * ((2 * lpm.phi + lpm.g) * screening.phi1
0243                            - lpm.g * screening.phi2 - lpm.phi * f_z)
0244                         / f10;
0245                 }
0246                 else
0247                 {
0248                     g = (this->screening_f1(delta) - f_z) / f10;
0249                 }
0250                 CELER_ASSERT(g > 0);
0251             }
0252             else
0253             {
0254                 // Used to sample from f2
0255                 epsilon = epsilon_min
0256                           + (half - epsilon_min) * generate_canonical(rng);
0257                 CELER_ASSERT(epsilon >= epsilon_min && epsilon <= half);
0258
0259                 // Calculate delta given the element atomic number and sampled
0260                 // epsilon
0261                 real_type delta = this->impact_parameter(epsilon);
0262                 CELER_ASSERT(delta <= delta_max && delta >= delta_min);
0263
0264                 // Calculate g_2 rejection function
0265                 if (enable_lpm_)
0266                 {
0267                     auto screening = screening_phi1_phi2(delta);
0268                     auto lpm = calc_lpm_functions_(epsilon);
0269                     g = lpm.xi
0270                         * ((lpm.phi + half * lpm.g) * screening.phi1
0271                            + half * lpm.g * screening.phi2
0272                            - half * (lpm.g + lpm.phi) * f_z)
0273                         / f20;
0274                 }
0275                 else
0276                 {
0277                     g = (this->screening_f2(delta) - f_z) / f20;
0278                 }
0279                 CELER_ASSERT(g > 0);
0280             }
0281         } while (g < generate_canonical(rng));
0282     }
0283
0284     // Construct interaction for change to primary (incident) particle (gamma)
0285     Interaction result = Interaction::from_absorption();
0286     result.secondaries = {secondaries, 2};
0287
0288     // Outgoing secondaries are electron and positron
0289     secondaries[0].particle_id = shared_.ids.electron;
0290     secondaries[1].particle_id = shared_.ids.positron;
0291
0292     secondaries[0].energy = Energy{(1 - epsilon) * inc_energy_
0293                                    - value_as<Mass>(shared_.electron_mass)};
0294     secondaries[1].energy = Energy{epsilon * inc_energy_
0295                                    - value_as<Mass>(shared_.electron_mass)};
0296
0297     // Select charges for child particles (e-, e+) randomly
0298     if (BernoulliDistribution(half)(rng))
0299     {
0300         trivial_swap(secondaries[0].energy, secondaries[1].energy);
0301     }
0302
0303     // Sample secondary directions.  Note that momentum is not exactly
0304     // conserved.
0305     real_type phi
0306         = UniformRealDistribution<real_type>(0, 2 * constants::pi)(rng);
0307
0308     // Electron
0309     TsaiUrbanDistribution sample_electron_angle(secondaries[0].energy,
0310                                                 shared_.electron_mass);
0311     real_type cost = sample_electron_angle(rng);
0312     secondaries[0].direction
0313         = rotate(from_spherical(cost, phi), inc_direction_);
0314     // Positron
0315     TsaiUrbanDistribution sample_positron_angle(secondaries[1].energy,
0316                                                 shared_.electron_mass);
0317     cost = sample_positron_angle(rng);
0318     secondaries[1].direction
0319         = rotate(from_spherical(cost, phi + constants::pi), inc_direction_);
0320     return result;
0321 }
0322
0323 //---------------------------------------------------------------------------//
0324 /*!
0325  * Screening variable \f$ \delta \f$.
0326  *
0327  * \f$ \delta \f$ is a function of \f$ \epsilon \f$ and is a measure of the
0328  * "impact parameter" of the incident photon.
0329  */
0330 CELER_FUNCTION real_type
0331 BetheHeitlerInteractor::impact_parameter(real_type eps) const
0332 {
0333     return 136 / element_.cbrt_z() * epsilon0_ / (eps * (1 - eps));
0334 }
0335
0336 //---------------------------------------------------------------------------//
0337 /*!
0338  * Screening functions \f$ \Phi_1(\delta) \f$ and \f$ \Phi_2(\delta) \f$.
0339  */
0340 CELER_FUNCTION auto BetheHeitlerInteractor::screening_phi1_phi2(
0341     real_type delta) const -> ScreeningFunctions
0342 {
0343     using R = real_type;
0344
0345     ScreeningFunctions result;
0346     if (delta > R(1.4))
0347     {
0348         result.phi1 = R(21.0190) - R(4.145) * std::log(delta + R(0.958));
0349         result.phi2 = result.phi1;
0350     }
0351     else
0352     {
0353         result.phi1 = R(20.806) - delta * (R(3.190) - R(0.5710) * delta);
0354         result.phi2 = R(20.234) - delta * (R(2.126) - R(0.0903) * delta);
0355     }
0356     return result;
0357 }
0358
0359 //---------------------------------------------------------------------------//
0360 /*!
0361  * Auxiliary screening functions \f$ F_1(\delta) \f$ and \f$ F_2(\delta) \f$.
0362  *
0363  * The functions \f$ F_1 = 3 \Phi_1(\delta) - \Phi_2(\delta) \f$
0364  * and \f$ F_2 = 1.5\Phi_1(\delta) - 0.5\Phi_2(\delta) \f$
0365  * are decreasing functions of \f$ \delta \f$ for all \f$ \delta \f$
0366  * in \f$ [\delta_\textrm{min}, \delta_\textrm{max}] \f$.
0367  * They reach their maximum value at
0368  * \f$ \delta_\textrm{min} = \delta(\epsilon = 1/2)\f$. They are used in the
0369  * composition + rejection technique for sampling \f$ \epsilon \f$.
0370  */
0371 CELER_FUNCTION real_type BetheHeitlerInteractor::screening_f1(real_type delta) const
0372 {
0373     using R = real_type;
0374
0375     return delta > R(1.4) ? R(42.038) - R(8.29) * std::log(delta + R(0.958))
0376                           : R(42.184) - delta * (R(7.444) - R(1.623) * delta);
0377 }
0378
0379 CELER_FUNCTION real_type BetheHeitlerInteractor::screening_f2(real_type delta) const
0380 {
0381     using R = real_type;
0382
0383     return delta > R(1.4) ? R(42.038) - R(8.29) * std::log(delta + R(0.958))
0384                           : R(41.326) - delta * (R(5.848) - R(0.902) * delta);
0385 }
0386 //---------------------------------------------------------------------------//
0387 }  // namespace celeritas