em/distribution/TsaiUrbanDistribution.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/TsaiUrbanDistribution.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/BernoulliDistribution.hh"
0015
0016 namespace celeritas
0017 {
0018 //---------------------------------------------------------------------------//
0019 /*!
0020  * Polar angular distribution for pair-production and bremsstrahlung processes.
0021  *
0022  * For pair production, the polar angle of the electron (or positron) is
0023  * defined with respect to the direction of the parent photon. The energy-
0024  * angle distribution given by Tsai is quite complicated to
0025  * sample and can be approximated by a density function suggested by Urban.
0026  *
0027  * The angular distribution of the emitted photons is obtained from a
0028  * simplified formula based on the Tsai cross-section,
0029  * which is expected to become isotropic in the low energy limit.
0030  *
0031  * \note This performs the same sampling routine as in Geant4's
0032  * ModifiedTsai class, based on derivation from Tsai (Rev Mod Phys 49,421(1977)
0033  * and documented in section 6.5.2 (pair-production), and 10.2.1 and 10.2.4
0034  * (bremsstrahlung) of the Geant4 Physics Reference (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