em/distribution/BraggICRU73QOEnergyDistribution.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/BraggICRU73QOEnergyDistribution.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 Bragg and ICRU73QO models, as
0030  * described in the Geant4 Physics Reference Manual release 11.2 section 11.1.
0031  */
0032 class BraggICRU73QOEnergyDistribution
0033 {
0034   public:
0035     //!@{
0036     //! \name Type aliases
0037     using Energy = units::MevEnergy;
0038     using Mass = units::MevMass;
0039     //!@}
0040
0041   public:
0042     // Construct with incident and exiting particle data
0043     inline CELER_FUNCTION
0044     BraggICRU73QOEnergyDistribution(ParticleTrackView const& particle,
0045                                     Energy electron_cutoff,
0046                                     Mass electron_mass);
0047
0048     // Sample the exiting energy
0049     template<class Engine>
0050     inline CELER_FUNCTION Energy operator()(Engine& rng);
0051
0052     //! Minimum energy of the secondary electron [MeV].
0053     CELER_FUNCTION Energy min_secondary_energy() const { return min_energy_; }
0054
0055     //! Maximum energy of the secondary electron [MeV].
0056     CELER_FUNCTION Energy max_secondary_energy() const { return max_energy_; }
0057
0058   private:
0059     //// DATA ////
0060
0061     // Incident partcle mass
0062     real_type inc_mass_;
0063     // Square of fractional speed of light for incident particle
0064     real_type beta_sq_;
0065     // Minimum energy of the secondary electron [MeV]
0066     Energy min_energy_;
0067     // Maximum energy of the secondary electron [MeV]
0068     Energy max_energy_;
0069
0070     //// CONSTANTS ////
0071
0072     //! Used in Bragg model to calculate minimum energy transfer to electron
0073     static CELER_CONSTEXPR_FUNCTION Energy bragg_lowest_kin_energy()
0074     {
0075         return units::MevEnergy{2.5e-4};  //! 0.25 keV
0076     }
0077
0078     //! Used in ICRU73QO model to calculate minimum energy transfer to electron
0079     static CELER_CONSTEXPR_FUNCTION Energy icru73qo_lowest_kin_energy()
0080     {
0081         return units::MevEnergy{5e-3};  //! 5 keV
0082     }
0083 };
0084
0085 //---------------------------------------------------------------------------//
0086 // INLINE DEFINITIONS
0087 //---------------------------------------------------------------------------//
0088 /*!
0089  * Construct with incident and exiting particle data.
0090  *
0091  * \todo Use proton mass from imported data instead of a constant
0092  */
0093 CELER_FUNCTION
0094 BraggICRU73QOEnergyDistribution::BraggICRU73QOEnergyDistribution(
0095     ParticleTrackView const& particle,
0096     Energy electron_cutoff,
0097     Mass electron_mass)
0098     : inc_mass_(value_as<Mass>(particle.mass()))
0099     , beta_sq_(particle.beta_sq())
0100     , min_energy_(
0101           min(electron_cutoff,
0102               Energy(value_as<Energy>(particle.charge() < zero_quantity()
0103                                           ? icru73qo_lowest_kin_energy()
0104                                           : bragg_lowest_kin_energy())
0105                      * inc_mass_
0106                      / native_value_to<Mass>(constants::proton_mass).value())))
0107     , max_energy_(detail::calc_max_secondary_energy(particle, electron_mass))
0108 {
0109 }
0110
0111 //---------------------------------------------------------------------------//
0112 /*!
0113  * Sample secondary electron energy.
0114  */
0115 template<class Engine>
0116 CELER_FUNCTION auto
0117 BraggICRU73QOEnergyDistribution::operator()(Engine& rng) -> Energy
0118 {
0119     InverseSquareDistribution sample_energy(value_as<Energy>(min_energy_),
0120                                             value_as<Energy>(max_energy_));
0121     real_type energy;
0122     do
0123     {
0124         // Sample 1/E^2 from Emin to Emax
0125         energy = sample_energy(rng);
0126     } while (RejectionSampler<>(
0127         1 - (beta_sq_ / value_as<Energy>(max_energy_)) * energy)(rng));
0128
0129     CELER_ENSURE(energy > 0);
0130     return Energy{energy};
0131 }
0132
0133 //---------------------------------------------------------------------------//
0134 }  // namespace celeritas