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