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