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
0060     BetheHeitlerInteractor(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 supression 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     CELER_EXPECT(particle.particle_id() == shared_.ids.gamma);
0152     CELER_EXPECT(value_as<Energy>(inc_energy_)
0153                  >= 2 * value_as<Mass>(shared_.electron_mass));
0154
0155     epsilon0_ = value_as<Mass>(shared_.electron_mass)
0156                 / value_as<Energy>(inc_energy_);
0157
0158     if (inc_energy_ < no_screening_threshold())
0159     {
0160         // Don't reject: just sample uniformly
0161         screen_delta_ = 0;
0162         f_z_ = 0;
0163     }
0164     else
0165     {
0166         screen_delta_ = epsilon0_ * 136 / element.cbrt_z();
0167         f_z_ = real_type(8) / real_type(3) * element.log_z();
0168         if (inc_energy_ > coulomb_corr_threshold())
0169         {
0170             // Apply Coulomb correction function
0171             f_z_ += 8 * element.coulomb_correction();
0172         }
0173     }
0174 }
0175
0176 //---------------------------------------------------------------------------//
0177 /*!
0178  * Sample the distribution.
0179  */
0180 template<class Engine>
0181 CELER_FUNCTION Interaction BetheHeitlerInteractor::operator()(Engine& rng)
0182 {
0183     // Allocate space for the electron/positron pair
0184     Secondary* secondaries = allocate_(2);
0185     if (secondaries == nullptr)
0186     {
0187         // Failed to allocate space for secondaries
0188         return Interaction::from_failure();
0189     }
0190
0191     constexpr real_type half = 0.5;
0192
0193     // Sample fraction of energy given to electron
0194     real_type epsilon;
0195     if (screen_delta_ == 0)
0196     {
0197         // If E_gamma < 2 MeV, rejection not needed -- just sample uniformly
0198         UniformRealDistribution<real_type> sample_eps(epsilon0_, half);
0199         epsilon = sample_eps(rng);
0200     }
0201     else
0202     {
0203         // Calculate the minimum (when \epsilon = 1/2) and maximum (when
0204         // \epsilon = \epsilon_1) values of screening variable, \delta. Above
0205         // 50 MeV, a Coulomb correction function is introduced.
0206         real_type const delta_min = 4 * screen_delta_;
0207         real_type const delta_max
0208             = std::exp((real_type(42.038) - f_z_) / real_type(8.29))
0209               - real_type(0.958);
0210         CELER_ASSERT(delta_min <= delta_max);
0211
0212         // Calculate the lower limit of epsilon. Due to the Coulomb correction,
0213         // the cross section can become negative even at kinematically allowed
0214         // \epsilon > \epsilon_0 values. To exclude these negative cross
0215         // sections, an additional constraint that \epsilon > \epsilon_1 is
0216         // introduced, where \epsilon_1 is the solution to
0217         // \Phi(\delta(\epsilon)) - F(Z)/2 = 0.
0218         real_type const epsilon1
0219             = half * (1 - std::sqrt(1 - delta_min / delta_max));
0220         real_type const epsilon_min = celeritas::max(epsilon0_, epsilon1);
0221
0222         // Decide to choose f1, g1 [brems] or f2, g2 [pair production]
0223         // based on N1, N2 (factors from corrected Bethe-Heitler cross section;
0224         // c.f. Eq. 6.6 of Geant4 Physics Reference 10.6)
0225         Real2 const fmin = this->screening_f(delta_min) - f_z_;
0226         BernoulliDistribution choose_f1g1(
0227             ipow<2>(half - epsilon_min) * fmin[0], real_type(1.5) * fmin[1]);
0228
0229         // Rejection function g_1 or g_2. Note the it's possible for g to be
0230         // greater than one
0231         real_type g;
0232         do
0233         {
0234             if (choose_f1g1(rng))
0235             {
0236                 // Used to sample from f1
0237                 epsilon = half
0238                           - (half - epsilon_min)
0239                                 * std::cbrt(generate_canonical(rng));
0240                 CELER_ASSERT(epsilon >= epsilon_min && epsilon <= half);
0241
0242                 // Calculate delta from element atomic number and sampled
0243                 // epsilon
0244                 real_type delta = this->impact_parameter(epsilon);
0245                 CELER_ASSERT(delta <= delta_max && delta >= delta_min);
0246
0247                 // Calculate g_1 rejection function
0248                 if (enable_lpm_)
0249                 {
0250                     auto screening = screening_phi(delta);
0251                     auto lpm = calc_lpm_functions_(epsilon);
0252                     g = lpm.xi
0253                         * ((2 * lpm.phi + lpm.g) * screening[0]
0254                            - lpm.g * screening[1] - lpm.phi * f_z_);
0255                 }
0256                 else
0257                 {
0258                     g = this->screening_f(delta)[0] - f_z_;
0259                 }
0260                 g /= fmin[0];
0261                 CELER_ASSERT(g > 0);
0262             }
0263             else
0264             {
0265                 // Used to sample from f2
0266                 epsilon = UniformRealDistribution{epsilon_min, half}(rng);
0267                 CELER_ASSERT(epsilon >= epsilon_min && epsilon <= half);
0268
0269                 // Calculate delta given the element atomic number and sampled
0270                 // epsilon
0271                 real_type delta = this->impact_parameter(epsilon);
0272                 CELER_ASSERT(delta <= delta_max && delta >= delta_min);
0273
0274                 // Calculate g_2 rejection function
0275                 if (enable_lpm_)
0276                 {
0277                     auto screening = screening_phi(delta);
0278                     auto lpm = calc_lpm_functions_(epsilon);
0279                     g = half * lpm.xi
0280                         * ((2 * lpm.phi + lpm.g) * screening[0]
0281                            + lpm.g * screening[1] - (lpm.g + lpm.phi) * f_z_);
0282                 }
0283                 else
0284                 {
0285                     g = this->screening_f(delta)[1] - f_z_;
0286                 }
0287                 g /= fmin[1];
0288                 CELER_ASSERT(g > 0);
0289             }
0290             // TODO: use rejection?
0291         } while (g < generate_canonical(rng));
0292     }
0293
0294     // Construct interaction for change to primary (incident) particle (gamma)
0295     Interaction result = Interaction::from_absorption();
0296     result.secondaries = {secondaries, 2};
0297
0298     // Outgoing secondaries are electron and positron with randomly selected
0299     // charges
0300     secondaries[0].particle_id = shared_.ids.electron;
0301     secondaries[1].particle_id = shared_.ids.positron;
0302
0303     // Incident energy split between the particles, with rest mass subtracted
0304     secondaries[0].energy = Energy{(1 - epsilon) * value_as<Energy>(inc_energy_)
0305                                    - value_as<Mass>(shared_.electron_mass)};
0306     secondaries[1].energy = Energy{epsilon * value_as<Energy>(inc_energy_)
0307                                    - value_as<Mass>(shared_.electron_mass)};
0308     if (BernoulliDistribution(half)(rng))
0309     {
0310         trivial_swap(secondaries[0].energy, secondaries[1].energy);
0311     }
0312
0313     // Sample secondary directions: note that momentum is not exactly conserved
0314     real_type const phi = UniformRealDistribution<real_type>(
0315         0, real_type(2 * constants::pi))(rng);
0316     auto sample_costheta = [&](Energy e) {
0317         return TsaiUrbanDistribution{e, shared_.electron_mass}(rng);
0318     };
0319
0320     // Note that positron has opposite azimuthal angle
0321     secondaries[0].direction
0322         = rotate(from_spherical(sample_costheta(secondaries[0].energy), phi),
0323                  inc_direction_);
0324     secondaries[1].direction
0325         = rotate(from_spherical(sample_costheta(secondaries[1].energy),
0326                                 phi + constants::pi),
0327                  inc_direction_);
0328
0329     return result;
0330 }
0331
0332 //---------------------------------------------------------------------------//
0333 /*!
0334  * Screening variable \f$ \delta \f$.
0335  *
0336  * \f$ \delta \f$ is a function of \f$ \epsilon \f$ and is a measure of the
0337  * "impact parameter" of the incident photon.
0338  */
0339 CELER_FUNCTION real_type
0340 BetheHeitlerInteractor::impact_parameter(real_type eps) const
0341 {
0342     return screen_delta_ / (eps * (1 - eps));
0343 }
0344
0345 //---------------------------------------------------------------------------//
0346 /*!
0347  * Screening functions \f$ \Phi_1(\delta) \f$ and \f$ \Phi_2(\delta) \f$.
0348  *
0349  * These correspond to \em f in Butcher: the first screening function is based
0350  * on the bremsstrahlung cross sections, and the second is due to pair
0351  * production. The values are improved function fits by M Novak (Geant4).
0352  */
0353 CELER_FUNCTION auto
0354 BetheHeitlerInteractor::screening_phi(real_type delta) -> Real2
0355 {
0356     using R = real_type;
0357
0358     Real2 result;
0359     if (delta > R(1.4))
0360     {
0361         result[0] = R(21.0190) - R(4.145) * std::log(delta + R(0.958));
0362         result[1] = result[0];
0363     }
0364     else
0365     {
0366         using PolyQuad = PolyEvaluator<real_type, 2>;
0367         result[0] = PolyQuad{20.806, -3.190, 0.5710}(delta);
0368         result[1] = PolyQuad{20.234, -2.126, 0.0903}(delta);
0369     }
0370     return result;
0371 }
0372
0373 //---------------------------------------------------------------------------//
0374 /*!
0375  * Auxiliary screening functions \f$ F_1(\delta) \f$ and \f$ F_2(\delta) \f$.
0376  *
0377  * The functions \f$ F_1 = 3 \Phi_1(\delta) - \Phi_2(\delta) \f$
0378  * and \f$ F_2 = 1.5\Phi_1(\delta) + 0.5\Phi_2(\delta) \f$
0379  * are decreasing functions of \f$ \delta \f$ for all \f$ \delta \f$
0380  * in \f$ [\delta_\textrm{min}, \delta_\textrm{max}] \f$.
0381  * They reach their maximum value at
0382  * \f$ \delta_\textrm{min} = \delta(\epsilon = 1/2)\f$. They are used in the
0383  * composition + rejection technique for sampling \f$ \epsilon \f$.
0384  *
0385  * Note that there's a typo in the Geant4 manual in the formula for F2:
0386  * subtraction should be addition.
0387  */
0388 CELER_FUNCTION auto
0389 BetheHeitlerInteractor::screening_f(real_type delta) -> Real2
0390 {
0391     using R = real_type;
0392     auto temp = screening_phi(delta);
0393     return {3 * temp[0] - temp[1], R{1.5} * temp[0] + R{0.5} * temp[1]};
0394 }
0395
0396 //---------------------------------------------------------------------------//
0397 }  // namespace celeritas