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