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  * Polar angular distribution for pair-production and bremsstrahlung processes.
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  * \note This performs the same sampling routine as in Geant4's
0031  * ModifiedTsai class, based on derivation from \cite{tsai-1974}
0032  * and documented in section 6.5.2 (pair-production), and 10.2.1 and 10.2.4
0033  * (bremsstrahlung) of the Geant4 Physics Reference \cite{g4prm}
0034  * (release 10.6).
0035  */
0036 class TsaiUrbanDistribution
0037 {
0038   public:
0039     //!@{
0040     //! \name Type aliases
0041     using Energy = units::MevEnergy;
0042     using Mass = units::MevMass;
0043     using result_type = real_type;
0044     //!@}
0045
0046   public:
0047     // Construct with defaults
0048     inline CELER_FUNCTION TsaiUrbanDistribution(Energy energy, Mass mass);
0049
0050     // Sample cos(theta) using the given random number generator
0051     template<class Engine>
0052     inline CELER_FUNCTION result_type operator()(Engine& rng);
0053
0054   private:
0055     // Dimensionless ratio of energy [Mev] to  mass * c^2 [MevMass*c^2]
0056     real_type umax_;
0057 };
0058
0059 //---------------------------------------------------------------------------//
0060 // INLINE DEFINITIONS
0061 //---------------------------------------------------------------------------//
0062 /*!
0063  * Construct from input data.
0064  */
0065 CELER_FUNCTION
0066 TsaiUrbanDistribution::TsaiUrbanDistribution(Energy energy, Mass mass)
0067     : umax_(2 * (1 + energy.value() / mass.value()))
0068 {
0069 }
0070
0071 //---------------------------------------------------------------------------//
0072 /*!
0073  * Sample the cosine of the polar angle of the exiting gamma.
0074  *
0075  * The z-axis is with respect to the direction of the parent particle.
0076  */
0077 template<class Engine>
0078 CELER_FUNCTION real_type TsaiUrbanDistribution::operator()(Engine& rng)
0079 {
0080     real_type u;
0081     do
0082     {
0083         real_type uu
0084             = -std::log(generate_canonical(rng) * generate_canonical(rng));
0085         u = uu
0086             * (BernoulliDistribution(0.25)(rng) ? real_type(1.6)
0087                                                 : real_type(1.6 / 3));
0088     } while (u > umax_);
0089
0090     return 1 - 2 * ipow<2>(u / umax_);
0091 }
0092
0093 //---------------------------------------------------------------------------//
0094 }  // namespace celeritas