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