em/distribution/BhabhaEnergyDistribution.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/distribution/BhabhaEnergyDistribution.hh
0007 //---------------------------------------------------------------------------//
0008 #pragma once
0009
0010 #include "corecel/Macros.hh"
0011 #include "corecel/Types.hh"
0012 #include "corecel/math/Algorithms.hh"
0013 #include "celeritas/Quantities.hh"
0014 #include "celeritas/random/distribution/RejectionSampler.hh"
0015 #include "celeritas/random/distribution/UniformRealDistribution.hh"
0016
0017 namespace celeritas
0018 {
0019 //---------------------------------------------------------------------------//
0020 /*!
0021  * Helper class for \c MollerBhabhaInteractor .
0022  *
0023  * Sample the exiting energy fraction for Bhabha scattering.
0024  */
0025 class BhabhaEnergyDistribution
0026 {
0027   public:
0028     //!@{
0029     //! \name Type aliases
0030     using Mass = units::MevMass;
0031     using Energy = units::MevEnergy;
0032     //!@}
0033
0034   public:
0035     // Construct with data from MollerBhabhaInteractor
0036     inline CELER_FUNCTION BhabhaEnergyDistribution(Mass electron_mass,
0037                                                    Energy min_valid_energy,
0038                                                    Energy inc_energy);
0039
0040     // Sample the exiting energy
0041     template<class Engine>
0042     inline CELER_FUNCTION real_type operator()(Engine& rng);
0043
0044   private:
0045     //// DATA ////
0046
0047     // Minimum energy fraction transferred to free electron
0048     real_type min_energy_fraction_;
0049     // Sampling parameter
0050     real_type gamma_;
0051
0052     //// HELPER FUNCTIONS ////
0053
0054     // Helper function for calculating rejection function g
0055     inline CELER_FUNCTION real_type calc_g_fraction(real_type epsilon_min,
0056                                                     real_type epsilon_max);
0057
0058     // Maximum energy fraction transferred to free electron [MeV]
0059     static CELER_CONSTEXPR_FUNCTION real_type max_energy_fraction()
0060     {
0061         return 1;
0062     }
0063
0064 };  // namespace BhabhaEnergyDistribution
0065
0066 //---------------------------------------------------------------------------//
0067 // INLINE DEFINITIONS
0068 //---------------------------------------------------------------------------//
0069 /*!
0070  * Construct with data from MollerBhabhaInteractor.
0071  */
0072 CELER_FUNCTION
0073 BhabhaEnergyDistribution::BhabhaEnergyDistribution(Mass electron_mass,
0074                                                    Energy min_valid_energy,
0075                                                    Energy inc_energy)
0076     : min_energy_fraction_(value_as<Energy>(min_valid_energy)
0077                            / value_as<Energy>(inc_energy))
0078     , gamma_(1 + value_as<Energy>(inc_energy) / value_as<Mass>(electron_mass))
0079 {
0080     CELER_EXPECT(electron_mass > zero_quantity()
0081                  && inc_energy > zero_quantity());
0082 }
0083
0084 //---------------------------------------------------------------------------//
0085 /*!
0086  * Sample epsilon for Bhabha scattering.
0087  */
0088 template<class Engine>
0089 CELER_FUNCTION real_type BhabhaEnergyDistribution::operator()(Engine& rng)
0090 {
0091     real_type const g_denominator = this->calc_g_fraction(
0092         min_energy_fraction_, this->max_energy_fraction());
0093
0094     UniformRealDistribution<> sample_inverse_epsilon(
0095         1 / this->max_energy_fraction(), 1 / min_energy_fraction_);
0096
0097     // Sample epsilon
0098     real_type epsilon;
0099     do
0100     {
0101         epsilon = 1 / sample_inverse_epsilon(rng);
0102     } while (RejectionSampler<>(this->calc_g_fraction(epsilon, epsilon),
0103                                 g_denominator)(rng));
0104
0105     return epsilon;
0106 }
0107
0108 //---------------------------------------------------------------------------//
0109 /*!
0110  * Evaluate the rejection function g.
0111  */
0112 CELER_FUNCTION real_type BhabhaEnergyDistribution::calc_g_fraction(
0113     real_type epsilon_min, real_type epsilon_max)
0114 {
0115     real_type const y = 1 / (1 + gamma_);
0116     real_type const y_sq = ipow<2>(y);
0117     real_type const one_minus_2y = 1 - 2 * y;
0118
0119     real_type const b1 = 2 - y_sq;
0120     real_type const b2 = one_minus_2y * (3 + y_sq);
0121     real_type const b4 = ipow<3>(one_minus_2y);
0122     real_type const b3 = ipow<2>(one_minus_2y) + b4;
0123     real_type const beta_sq = 1 - (1 / ipow<2>(gamma_));
0124
0125     return 1
0126            + (ipow<4>(epsilon_max) * b4 - ipow<3>(epsilon_min) * b3
0127               + ipow<2>(epsilon_max) * b2 - epsilon_min * b1)
0128                  * beta_sq;
0129 }
0130
0131 //---------------------------------------------------------------------------//
0132 }  // namespace celeritas