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