em/distribution/BetheBlochEnergyDistribution.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/BetheBlochEnergyDistribution.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/Constants.hh"
0017 #include "celeritas/Quantities.hh"
0018 #include "celeritas/phys/ParticleTrackView.hh"
0019
0020 #include "detail/Utils.hh"
0021
0022 namespace celeritas
0023 {
0024 //---------------------------------------------------------------------------//
0025 /*!
0026  * Sample the energy of the delta ray for muon or hadron ionization.
0027  *
0028  * This samples the energy according to the Bethe-Bloch model, as described in
0029  * the Geant4 Physics Reference Manual release 11.2 section 12.1.5. The
0030  * Bethe-Bloch differential cross section can be written as
0031  * \f[
0032    \difd{\sigma}{T} = 2\pi r_e^2 mc^2 Z \frac{z_p^2}{\beta^2}\frac{1}{T^2}
0033    \left[1 - \beta^2 \frac{T}{T_{max}} + s \frac{T^2}{2E^2} \right]
0034  * \f]
0035  * and factorized as
0036  * \f[
0037    \difd{\sigma}{T} = C f(T) g(T)
0038  * \f]
0039  * with \f$ T \in [T_{cut}, T_{max}] \f$, where \f$ f(T) = \frac{1}{T^2} \f$,
0040  * \f$ g(T) = 1 - \beta^2 \frac{T}{T_max} + s \frac{T^2}{2 E^2} \f$, \f$ T \f$
0041  * is the kinetic energy of the electron, \f$ E \f$ is the total energy of the
0042  * incident particle, and \f$ s \f$ is 0 for spinless particles and 1
0043  * otherwise. The energy is sampled from \f$ f(T) \f$ and accepted with
0044  * probability \f$ g(T) \f$.
0045  */
0046 class BetheBlochEnergyDistribution
0047 {
0048   public:
0049     //!@{
0050     //! \name Type aliases
0051     using Energy = units::MevEnergy;
0052     using Mass = units::MevMass;
0053     //!@}
0054
0055   public:
0056     // Construct with incident and exiting particle data
0057     inline CELER_FUNCTION
0058     BetheBlochEnergyDistribution(ParticleTrackView const& particle,
0059                                  Energy electron_cutoff,
0060                                  Mass electron_mass);
0061
0062     // Sample the exiting energy
0063     template<class Engine>
0064     inline CELER_FUNCTION Energy operator()(Engine& rng);
0065
0066     //! Minimum energy of the secondary electron [MeV].
0067     CELER_FUNCTION Energy min_secondary_energy() const { return min_energy_; }
0068
0069     //! Maximum energy of the secondary electron [MeV].
0070     CELER_FUNCTION Energy max_secondary_energy() const { return max_energy_; }
0071
0072   private:
0073     // Incident partcle mass
0074     real_type inc_mass_;
0075     // Square of fractional speed of light for incident particle
0076     real_type beta_sq_;
0077     // Minimum energy of the secondary electron [MeV]
0078     Energy min_energy_;
0079     // Maximum energy of the secondary electron [MeV]
0080     Energy max_energy_;
0081 };
0082
0083 //---------------------------------------------------------------------------//
0084 // INLINE DEFINITIONS
0085 //---------------------------------------------------------------------------//
0086 /*!
0087  * Construct with incident and exiting particle data.
0088  */
0089 CELER_FUNCTION
0090 BetheBlochEnergyDistribution::BetheBlochEnergyDistribution(
0091     ParticleTrackView const& particle,
0092     Energy electron_cutoff,
0093     Mass electron_mass)
0094     : inc_mass_(value_as<Mass>(particle.mass()))
0095     , beta_sq_(particle.beta_sq())
0096     , min_energy_(electron_cutoff)
0097     , max_energy_(detail::calc_max_secondary_energy(particle, electron_mass))
0098 {
0099 }
0100
0101 //---------------------------------------------------------------------------//
0102 /*!
0103  * Sample secondary electron energy.
0104  */
0105 template<class Engine>
0106 CELER_FUNCTION auto
0107 BetheBlochEnergyDistribution::operator()(Engine& rng) -> Energy
0108 {
0109     InverseSquareDistribution sample_energy(value_as<Energy>(min_energy_),
0110                                             value_as<Energy>(max_energy_));
0111     real_type energy;
0112     do
0113     {
0114         // Sample 1/E^2 from Emin to Emax
0115         energy = sample_energy(rng);
0116         /*!
0117          * \todo Adjust rejection functions if particle has positive spin
0118          */
0119     } while (RejectionSampler<>(
0120         1 - (beta_sq_ / value_as<Energy>(max_energy_)) * energy)(rng));
0121
0122     /*!
0123      * \todo For hadrons, suppress high energy delta ray production with the
0124      * projectile form factor
0125      */
0126
0127     return Energy{energy};
0128 }
0129
0130 //---------------------------------------------------------------------------//
0131 }  // namespace celeritas