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