em/distribution/MollerEnergyDistribution.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/MollerEnergyDistribution.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 Moller scattering.
0023  */
0024 class MollerEnergyDistribution
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 MollerEnergyDistribution(Mass electron_mass,
0036                                                    Energy min_valid_energy,
0037                                                    Energy inc_energy);
0038
0039     // Sample the exiting energy fraction
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);
0055     // Maximum energy fraction transferred to free electron [MeV]
0056     static CELER_CONSTEXPR_FUNCTION real_type max_energy_fraction()
0057     {
0058         return 0.5;
0059     }
0060 };
0061
0062 //---------------------------------------------------------------------------//
0063 // INLINE DEFINITIONS
0064 //---------------------------------------------------------------------------//
0065 /*!
0066  * Construct with data from MollerBhabhaInteractor.
0067  */
0068 CELER_FUNCTION
0069 MollerEnergyDistribution::MollerEnergyDistribution(Mass electron_mass,
0070                                                    Energy min_valid_energy,
0071                                                    Energy inc_energy)
0072     : min_energy_fraction_(value_as<Energy>(min_valid_energy)
0073                            / value_as<Energy>(inc_energy))
0074     , gamma_(1 + value_as<Energy>(inc_energy) / value_as<Mass>(electron_mass))
0075 {
0076     CELER_EXPECT(electron_mass > zero_quantity()
0077                  && inc_energy > zero_quantity());
0078 }
0079
0080 //---------------------------------------------------------------------------//
0081 /*!
0082  * Sample epsilon for Moller scattering.
0083  */
0084 template<class Engine>
0085 CELER_FUNCTION real_type MollerEnergyDistribution::operator()(Engine& rng)
0086 {
0087     real_type const g_denominator
0088         = this->calc_g_fraction(this->max_energy_fraction());
0089
0090     UniformRealDistribution<> sample_inverse_epsilon(
0091         1 / this->max_energy_fraction(), 1 / min_energy_fraction_);
0092
0093     // Sample fraction of exiting energy
0094     real_type epsilon;
0095     do
0096     {
0097         epsilon = 1 / sample_inverse_epsilon(rng);
0098     } while (RejectionSampler<>(calc_g_fraction(epsilon), g_denominator)(rng));
0099
0100     return epsilon;
0101 }
0102
0103 //---------------------------------------------------------------------------//
0104 /*!
0105  * Evaluate the rejection function g.
0106  */
0107 CELER_FUNCTION real_type
0108 MollerEnergyDistribution::calc_g_fraction(real_type epsilon)
0109 {
0110     real_type const two_gamma_term = (2 * gamma_ - 1) / ipow<2>(gamma_);
0111     real_type const complement_frac = 1 - epsilon;
0112
0113     return 1 - two_gamma_term * epsilon
0114            + ipow<2>(epsilon)
0115                  * (1 - two_gamma_term
0116                     + (1 - two_gamma_term * complement_frac)
0117                           / ipow<2>(complement_frac));
0118 }
0119
0120 //---------------------------------------------------------------------------//
0121 }  // namespace celeritas