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