em/distribution/TsaiUrbanDistribution.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/TsaiUrbanDistribution.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/BernoulliDistribution.hh"
0013 #include "celeritas/Quantities.hh"
0014
0015 namespace celeritas
0016 {
0017 //---------------------------------------------------------------------------//
0018 /*!
0019  * Sample exiting polar angle for pair production and bremsstrahlung.
0020  *
0021  * For pair production, the polar angle of the electron (or positron) is
0022  * defined with respect to the direction of the parent photon. The energy-
0023  * angle distribution given by Tsai is quite complicated to
0024  * sample and can be approximated by a density function suggested by Urban.
0025  *
0026  * The angular distribution of the emitted photons is obtained from a
0027  * simplified formula based on the Tsai cross-section,
0028  * which is expected to become isotropic in the low energy limit.
0029  *
0030  * The sample result is the cosine of the exiting angle with respect to the
0031  * incident angle.
0032  *
0033  * \note This performs the same sampling routine as in Geant4's
0034  * ModifiedTsai class, based on derivation from \cite{tsai-1974}
0035  * and documented in section 6.5.2 (pair-production), and 10.2.1 and 10.2.4
0036  * (bremsstrahlung) of the Geant4 Physics Reference \cite{g4prm}
0037  * (release 10.6).
0038  */
0039 class TsaiUrbanDistribution
0040 {
0041   public:
0042     //!@{
0043     //! \name Type aliases
0044     using Energy = units::MevEnergy;
0045     using Mass = units::MevMass;
0046     using result_type = real_type;
0047     //!@}
0048
0049   public:
0050     // Construct with defaults
0051     inline CELER_FUNCTION TsaiUrbanDistribution(Energy energy, Mass mass);
0052
0053     // Sample cos(theta) using the given random number generator
0054     template<class Engine>
0055     inline CELER_FUNCTION result_type operator()(Engine& rng);
0056
0057   private:
0058     // Dimensionless ratio of energy [Mev] to  mass * c^2 [MevMass*c^2]
0059     real_type umax_;
0060 };
0061
0062 //---------------------------------------------------------------------------//
0063 // INLINE DEFINITIONS
0064 //---------------------------------------------------------------------------//
0065 /*!
0066  * Construct from input data.
0067  */
0068 CELER_FUNCTION
0069 TsaiUrbanDistribution::TsaiUrbanDistribution(Energy energy, Mass mass)
0070     : umax_(2 * (1 + energy.value() / mass.value()))
0071 {
0072 }
0073
0074 //---------------------------------------------------------------------------//
0075 /*!
0076  * Sample the cosine of the polar angle of the exiting gamma.
0077  */
0078 template<class Engine>
0079 CELER_FUNCTION real_type TsaiUrbanDistribution::operator()(Engine& rng)
0080 {
0081     real_type u;
0082     do
0083     {
0084         real_type uu
0085             = -std::log(generate_canonical(rng) * generate_canonical(rng));
0086         u = uu
0087             * (BernoulliDistribution(0.25)(rng) ? real_type(1.6)
0088                                                 : real_type(1.6 / 3));
0089     } while (u > umax_);
0090
0091     return 1 - 2 * ipow<2>(u / umax_);
0092 }
0093
0094 //---------------------------------------------------------------------------//
0095 }  // namespace celeritas