em/distribution/MuBBEnergyDistribution.hh

0001 //----------------------------------*-C++-*----------------------------------//
0002 // Copyright 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/MuBBEnergyDistribution.hh
0007 //---------------------------------------------------------------------------//
0008 #pragma once
0009
0010 #include <cmath>
0011
0012 #include "corecel/Macros.hh"
0013 #include "corecel/Types.hh"
0014 #include "corecel/math/Algorithms.hh"
0015 #include "celeritas/Quantities.hh"
0016 #include "celeritas/phys/ParticleTrackView.hh"
0017 #include "celeritas/random/distribution/InverseSquareDistribution.hh"
0018 #include "celeritas/random/distribution/RejectionSampler.hh"
0019
0020 #include "detail/Utils.hh"
0021
0022 namespace celeritas
0023 {
0024 //---------------------------------------------------------------------------//
0025 /*!
0026  * Sample delta ray energy for the muon Bethe-Bloch ionization model.
0027  8
0028  * This samples the energy according to the muon Bethe-Bloch model, as
0029  * described in the Geant4 Physics Reference Manual release 11.2 section 11.1.
0030  * At the higher energies for which this model is applied, leading radiative
0031  * corrections are taken into account. The differential cross section can be
0032  * written as
0033  * \f[
0034    \sigma(E, \epsilon) = \sigma_{BB}(E, \epsilon)\left[1 + \frac{\alpha}{2\pi}
0035    \log \left(1 + \frac{2\epsilon}{m_e} \log \left(\frac{4 m_e E(E -
0036    \epsilon}{m_{\mu}^2(2\epsilon + m_e)} \right) \right) \right].
0037  * \f]
0038  * \f$ \sigma_{BB}(E, \epsilon) \f$ is the Bethe-Bloch cross section, \f$ m_e
0039  * \f$ is the electron mass, \f$ m_{\mu} \f$ is the muon mass, \f$ E \f$ is the
0040  * total energy of the muon, and \f$ \epsilon = \omega + T \f$ is the energy
0041  * transfer, where \f$ T \f$ is the kinetic energy of the electron and \f$
0042  * \omega \f$ is the energy of the radiative gamma (which is neglected).
0043  *
0044  * As in the Bethe-Bloch model, the energy is sampled by factorizing the cross
0045  * section as \f$ \sigma = C f(T) g(T) \f$, where \f$ f(T) = \frac{1}{T^2} \f$
0046  * and \f$ T \in [T_{cut}, T_{max}] \f$. The energy is sampled from \f$ f(T)
0047  * \f$ and accepted with probability \f$ g(T) \f$.
0048  */
0049 class MuBBEnergyDistribution
0050 {
0051   public:
0052     //!@{
0053     //! \name Type aliases
0054     using Energy = units::MevEnergy;
0055     using Mass = units::MevMass;
0056     //!@}
0057
0058   public:
0059     // Construct with incident and exiting particle data
0060     inline CELER_FUNCTION
0061     MuBBEnergyDistribution(ParticleTrackView const& particle,
0062                            Energy electron_cutoff,
0063                            Mass electron_mass);
0064
0065     // Sample the exiting energy
0066     template<class Engine>
0067     inline CELER_FUNCTION Energy operator()(Engine& rng);
0068
0069     //! Minimum energy of the secondary electron [MeV].
0070     CELER_FUNCTION Energy min_secondary_energy() const { return min_energy_; }
0071
0072     //! Maximum energy of the secondary electron [MeV].
0073     CELER_FUNCTION Energy max_secondary_energy() const { return max_energy_; }
0074
0075   private:
0076     //// DATA ////
0077
0078     // Incident partcle mass
0079     real_type inc_mass_;
0080     // Total energy of the incident particle [MeV]
0081     real_type total_energy_;
0082     // Square of fractional speed of light for incident particle
0083     real_type beta_sq_;
0084     // Secondary electron mass
0085     real_type electron_mass_;
0086     // Secondary electron cutoff energy [MeV]
0087     Energy min_energy_;
0088     // Maximum energy of the secondary electron [MeV]
0089     Energy max_energy_;
0090     // Whether to apply the radiative correction
0091     bool use_rad_correction_;
0092     // Envelope distribution for rejection sampling
0093     real_type envelope_;
0094
0095     //// CONSTANTS ////
0096
0097     //! Incident energy above which the radiative correction is applied [MeV]
0098     static CELER_CONSTEXPR_FUNCTION Energy rad_correction_limit()
0099     {
0100         return units::MevEnergy{250};  //!< 250 MeV
0101     }
0102
0103     //! Lower limit of radiative correction integral [MeV]
0104     static CELER_CONSTEXPR_FUNCTION Energy kin_energy_limit()
0105     {
0106         // This is the lower limit of the energy transfer from the incident
0107         // muon to the delta ray and radiative gammas
0108         return units::MevEnergy{0.1};  //!< 100 keV
0109     }
0110
0111     //! Fine structure constant over two pi
0112     static CELER_CONSTEXPR_FUNCTION real_type alpha_over_twopi()
0113     {
0114         return constants::alpha_fine_structure / (2 * constants::pi);
0115     }
0116
0117     //// HELPER FUNCTIONS ////
0118
0119     inline CELER_FUNCTION real_type calc_envelope_distribution() const;
0120 };
0121
0122 //---------------------------------------------------------------------------//
0123 // INLINE DEFINITIONS
0124 //---------------------------------------------------------------------------//
0125 /*!
0126  * Construct with incident and exiting particle data.
0127  */
0128 CELER_FUNCTION
0129 MuBBEnergyDistribution::MuBBEnergyDistribution(ParticleTrackView const& particle,
0130                                                Energy electron_cutoff,
0131                                                Mass electron_mass)
0132     : inc_mass_(value_as<Mass>(particle.mass()))
0133     , total_energy_(value_as<Energy>(particle.total_energy()))
0134     , beta_sq_(particle.beta_sq())
0135     , electron_mass_(value_as<Mass>(electron_mass))
0136     , min_energy_(electron_cutoff)
0137     , max_energy_(detail::calc_max_secondary_energy(particle, electron_mass))
0138     , use_rad_correction_(particle.energy() > rad_correction_limit()
0139                           && max_energy_ > kin_energy_limit())
0140     , envelope_(this->calc_envelope_distribution())
0141 {
0142 }
0143
0144 //---------------------------------------------------------------------------//
0145 /*!
0146  * Sample secondary electron energy.
0147  */
0148 template<class Engine>
0149 CELER_FUNCTION auto MuBBEnergyDistribution::operator()(Engine& rng) -> Energy
0150 {
0151     InverseSquareDistribution sample_energy(value_as<Energy>(min_energy_),
0152                                             value_as<Energy>(max_energy_));
0153     real_type energy;
0154     real_type target;
0155     do
0156     {
0157         energy = sample_energy(rng);
0158         target = 1 - (beta_sq_ / value_as<Energy>(max_energy_)) * energy
0159                  + real_type(0.5) * ipow<2>(energy / total_energy_);
0160
0161         if (use_rad_correction_
0162             && energy > value_as<Energy>(kin_energy_limit()))
0163         {
0164             real_type a1 = std::log(1 + 2 * energy / electron_mass_);
0165             real_type a3 = std::log(4 * total_energy_ * (total_energy_ - energy)
0166                                     / ipow<2>(inc_mass_));
0167             target *= (1 + alpha_over_twopi() * a1 * (a3 - a1));
0168         }
0169     } while (RejectionSampler<>(target, envelope_)(rng));
0170
0171     CELER_ENSURE(energy > 0);
0172     return Energy{energy};
0173 }
0174
0175 //---------------------------------------------------------------------------//
0176 /*!
0177  * Calculate envelope maximum for rejection sampling of secondary energy.
0178  */
0179 CELER_FUNCTION real_type MuBBEnergyDistribution::calc_envelope_distribution() const
0180 {
0181     if (use_rad_correction_)
0182     {
0183         return 1
0184                + alpha_over_twopi()
0185                      * ipow<2>(std::log(2 * total_energy_ / inc_mass_));
0186     }
0187     return 1;
0188 }
0189
0190 //---------------------------------------------------------------------------//
0191 }  // namespace celeritas