From 7bbf401cf4bed862b18d83fbb0f14cfdc1723b19 Mon Sep 17 00:00:00 2001 From: Amanda Lund Date: Fri, 27 Sep 2024 07:04:29 -0500 Subject: [PATCH] Add Bethe-Bloch model and refactor muon ionization (#1425) * Add Bethe-Bloch ionization model * Refactor muon ionization - Use a single interactor/executor templated on energy distribution - Construct models with applicability to allow different incident particles and energy bounds - Add a helper function to calculate max delta ray energy * Remove incident particle IDs from muon/hadron ionization data * Clean up documentation * Update documentation --- doc/implementation/core-physics/stepping.rst | 2 +- doc/implementation/em-physics.rst | 32 +- src/celeritas/CMakeLists.txt | 1 + src/celeritas/em/data/BraggICRU73QOData.hh | 47 --- ...theBlochData.hh => MuHadIonizationData.hh} | 18 +- .../BetheBlochEnergyDistribution.hh | 132 +++++++ .../BraggICRU73QOEnergyDistribution.hh | 134 +++++++ .../em/distribution/MuBBEnergyDistribution.hh | 82 +++-- src/celeritas/em/distribution/detail/Utils.hh | 45 +++ .../em/executor/MuBetheBlochExecutor.hh | 48 --- ...Executor.hh => MuHadIonizationExecutor.hh} | 18 +- .../em/interactor/MuBetheBlochInteractor.hh | 179 ---------- ...ractor.hh => MuHadIonizationInteractor.hh} | 98 ++--- .../em/interactor/detail/PhysicsConstants.hh | 18 - src/celeritas/em/model/BetheBlochModel.cc | 89 +++++ src/celeritas/em/model/BetheBlochModel.cu | 38 ++ src/celeritas/em/model/BetheBlochModel.hh | 53 +++ src/celeritas/em/model/BraggModel.cc | 39 +- src/celeritas/em/model/BraggModel.cu | 7 +- src/celeritas/em/model/BraggModel.hh | 13 +- src/celeritas/em/model/ICRU73QOModel.cc | 39 +- src/celeritas/em/model/ICRU73QOModel.cu | 7 +- src/celeritas/em/model/ICRU73QOModel.hh | 13 +- src/celeritas/em/model/MuBetheBlochModel.cc | 41 +-- src/celeritas/em/model/MuBetheBlochModel.cu | 6 +- src/celeritas/em/model/MuBetheBlochModel.hh | 13 +- .../em/model/detail/MuHadIonizationBuilder.hh | 98 +++++ .../em/process/MuIonizationProcess.cc | 51 ++- .../em/process/MuIonizationProcess.hh | 9 +- .../ext/detail/CelerEmStandardPhysics.cc | 2 +- src/celeritas/phys/ImportedProcessAdapter.hh | 18 + src/orange/detail/OrientedBoundingZone.hh | 8 +- test/celeritas/CMakeLists.txt | 1 + test/celeritas/em/BetheBloch.test.cc | 334 ++++++++++++++++++ test/celeritas/em/BraggICRU73QO.test.cc | 147 +++++--- test/celeritas/em/MuBetheBloch.test.cc | 103 +++++- test/celeritas/phys/ProcessBuilder.test.cc | 11 +- 37 files changed, 1406 insertions(+), 588 deletions(-) delete mode 100644 src/celeritas/em/data/BraggICRU73QOData.hh rename src/celeritas/em/data/{MuBetheBlochData.hh => MuHadIonizationData.hh} (72%) create mode 100644 src/celeritas/em/distribution/BetheBlochEnergyDistribution.hh create mode 100644 src/celeritas/em/distribution/BraggICRU73QOEnergyDistribution.hh create mode 100644 src/celeritas/em/distribution/detail/Utils.hh delete mode 100644 src/celeritas/em/executor/MuBetheBlochExecutor.hh rename src/celeritas/em/executor/{BraggICRU73QOExecutor.hh => MuHadIonizationExecutor.hh} (74%) delete mode 100644 src/celeritas/em/interactor/MuBetheBlochInteractor.hh rename src/celeritas/em/interactor/{BraggICRU73QOInteractor.hh => MuHadIonizationInteractor.hh} (55%) create mode 100644 src/celeritas/em/model/BetheBlochModel.cc create mode 100644 src/celeritas/em/model/BetheBlochModel.cu create mode 100644 src/celeritas/em/model/BetheBlochModel.hh create mode 100644 src/celeritas/em/model/detail/MuHadIonizationBuilder.hh create mode 100644 test/celeritas/em/BetheBloch.test.cc diff --git a/doc/implementation/core-physics/stepping.rst b/doc/implementation/core-physics/stepping.rst index 959850b967..618e62b892 100644 --- a/doc/implementation/core-physics/stepping.rst +++ b/doc/implementation/core-physics/stepping.rst @@ -64,7 +64,7 @@ stepping loop: :cpp:class:`CoreBeginRunActionInterface` is called once per event is called once per step, ordered using :cpp:enum:`celeritas::StepActionOrder`. .. doxygenclass:: celeritas::ActionInterface -.. doxygenclass:: celeritas::CoreBeginRunActionInterface +.. doxygentypedef:: celeritas::CoreBeginRunActionInterface .. doxygenclass:: celeritas::StepActionInterface .. doxygenenum:: celeritas::StepActionOrder diff --git a/doc/implementation/em-physics.rst b/doc/implementation/em-physics.rst index 52020530df..e65dcd6f8f 100644 --- a/doc/implementation/em-physics.rst +++ b/doc/implementation/em-physics.rst @@ -54,15 +54,19 @@ The following table summarizes the EM processes and models in Celeritas. | +---------------------+-----------------------------+-----------------------------------------------------+--------------------------+ | | Rayleigh scattering | Livermore | :cpp:class:`celeritas::RayleighInteractor` | 0--100 TeV | +----------------+---------------------+-----------------------------+-----------------------------------------------------+--------------------------+ - | :math:`\mu^-` | Ionization | ICRU73QO | :cpp:class:`celeritas::BraggICRU73QOInteractor` | 0--200 keV | + | :math:`\mu^-` | Ionization | ICRU73QO | :cpp:class:`celeritas::MuHadIonizationInteractor` | 0--200 keV | | + +-----------------------------+-----------------------------------------------------+--------------------------+ - | | | Bethe--Bloch | :cpp:class:`celeritas::MuBetheBlochInteractor` | 200 keV--100 TeV | + | | | Bethe--Bloch | :cpp:class:`celeritas::MuHadIonizationInteractor` | 200 keV--1 GeV | + | + +-----------------------------+-----------------------------------------------------+--------------------------+ + | | | Mu Bethe--Bloch | :cpp:class:`celeritas::MuHadIonizationInteractor` | 200 keV--100 TeV | | +---------------------+-----------------------------+-----------------------------------------------------+--------------------------+ | | Bremsstrahlung | Mu bremsstrahlung | :cpp:class:`celeritas::MuBremsstrahlungInteractor` | 0--100 TeV | +----------------+---------------------+-----------------------------+-----------------------------------------------------+--------------------------+ - | :math:`\mu^+` | Ionization | Bragg | :cpp:class:`celeritas::BraggICRU73QOInteractor` | 0--200 keV | + | :math:`\mu^+` | Ionization | Bragg | :cpp:class:`celeritas::MuHadIonizationInteractor` | 0--200 keV | + | + +-----------------------------+-----------------------------------------------------+--------------------------+ + | | | Bethe--Bloch | :cpp:class:`celeritas::MuHadIonizationInteractor` | 200 keV--1 GeV | | + +-----------------------------+-----------------------------------------------------+--------------------------+ - | | | Bethe--Bloch | :cpp:class:`celeritas::MuBetheBlochInteractor` | 200 keV--100 TeV | + | | | Mu Bethe--Bloch | :cpp:class:`celeritas::MuHadIonizationInteractor` | 200 keV--100 TeV | | +---------------------+-----------------------------+-----------------------------------------------------+--------------------------+ | | Bremsstrahlung | Mu bremsstrahlung | :cpp:class:`celeritas::MuBremsstrahlungInteractor` | 0--100 TeV | +----------------+---------------------+-----------------------------+-----------------------------------------------------+--------------------------+ @@ -109,15 +113,19 @@ The following table summarizes the EM processes and models in Celeritas. \cline{2-5} & Rayleigh scattering & Livermore & \texttt{\scriptsize celeritas::RayleighInteractor} & 0--100 TeV \\ \hline - \multirow{3}{*}{$\mu^-$} & \multirow{2}{*}{Ionization} & ICRU73QO & \texttt{\scriptsize celeritas::BraggICRU73QOInteractor} & 0--200 keV \\ + \multirow{3}{*}{$\mu^-$} & \multirow{2}{*}{Ionization} & ICRU73QO & \texttt{\scriptsize celeritas::MuHadIonizationInteractor} & 0--200 keV \\ \cline{3-5} - & & Bethe--Bloch & \texttt{\scriptsize celeritas::MuBetheBlochInteractor} & 200 keV -- 100 TeV \\ + & & Bethe--Bloch & \texttt{\scriptsize celeritas::MuHadIonizationInteractor} & 200 keV -- 1 GeV \\ + \cline{3-5} + & & Mu Bethe--Bloch & \texttt{\scriptsize celeritas::MuHadIonizationInteractor} & 200 keV -- 100 TeV \\ \cline{2-5} & Bremsstrahlung & Mu bremsstrahlung & \texttt{\scriptsize celeritas::MuBremsstrahlungInteractor} & 0--100 TeV \\ \hline - \multirow{3}{*}{$\mu^+$} & \multirow{2}{*}{Ionization} & Bragg & \texttt{\scriptsize celeritas::BraggICRU73QOInteractor} & 0--200 keV \\ + \multirow{3}{*}{$\mu^+$} & \multirow{2}{*}{Ionization} & Bragg & \texttt{\scriptsize celeritas::MuHadIonizationInteractor} & 0--200 keV \\ + \cline{3-5} + & & Bethe--Bloch & \texttt{\scriptsize celeritas::MuHadIonizationInteractor} & 200 keV -- 1 GeV \\ \cline{3-5} - & & Bethe--Bloch & \texttt{\scriptsize celeritas::MuBetheBlochInteractor} & 200 keV -- 100 TeV \\ + & & Mu Bethe--Bloch & \texttt{\scriptsize celeritas::MuHadIonizationInteractor} & 200 keV -- 100 TeV \\ \cline{2-5} & Bremsstrahlung & Mu bremsstrahlung & \texttt{\scriptsize celeritas::MuBremsstrahlungInteractor} & 0--100 TeV \\ \hline @@ -154,13 +162,14 @@ as properties of any secondary particles produced. Ionization ---------- -.. doxygenclass:: celeritas::BraggICRU73QOInteractor .. doxygenclass:: celeritas::MollerBhabhaInteractor -.. doxygenclass:: celeritas::MuBetheBlochInteractor +.. doxygenclass:: celeritas::MuHadIonizationInteractor The exiting energy distribution from most of these ionization models are sampled using external helper distributions. +.. doxygenclass:: celeritas::BetheBlochEnergyDistribution +.. doxygenclass:: celeritas::BraggICRU73QOEnergyDistribution .. doxygenclass:: celeritas::BhabhaEnergyDistribution .. doxygenclass:: celeritas::MollerEnergyDistribution .. doxygenclass:: celeritas::MuBBEnergyDistribution @@ -280,4 +289,5 @@ to provide cross sections, setup options, and other data to the EM physics. .. doxygenstruct:: celeritas::ImportLivermoreSubshell .. doxygenstruct:: celeritas::ImportLivermorePE -.. doxygenstruct:: celeritas::ImportSBTable +.. doxygenstruct:: celeritas::ImportMuPairProductionTable +.. doxygentypedef:: celeritas::ImportSBTable diff --git a/src/celeritas/CMakeLists.txt b/src/celeritas/CMakeLists.txt index f4a8e400ec..eb5c8caf0b 100644 --- a/src/celeritas/CMakeLists.txt +++ b/src/celeritas/CMakeLists.txt @@ -253,6 +253,7 @@ macro(celeritas_polysource filename_we) endmacro() celeritas_polysource(em/model/BetheHeitlerModel) +celeritas_polysource(em/model/BetheBlochModel) celeritas_polysource(em/model/BraggModel) celeritas_polysource(em/model/CombinedBremModel) celeritas_polysource(em/model/EPlusGGModel) diff --git a/src/celeritas/em/data/BraggICRU73QOData.hh b/src/celeritas/em/data/BraggICRU73QOData.hh deleted file mode 100644 index bb3b0b9f18..0000000000 --- a/src/celeritas/em/data/BraggICRU73QOData.hh +++ /dev/null @@ -1,47 +0,0 @@ -//----------------------------------*-C++-*----------------------------------// -// Copyright 2024 UT-Battelle, LLC, and other Celeritas developers. -// See the top-level COPYRIGHT file for details. -// SPDX-License-Identifier: (Apache-2.0 OR MIT) -//---------------------------------------------------------------------------// -//! \file celeritas/em/data/BraggICRU73QOData.hh -//---------------------------------------------------------------------------// -#pragma once - -#include "corecel/Macros.hh" -#include "corecel/Types.hh" -#include "celeritas/Quantities.hh" -#include "celeritas/Types.hh" - -namespace celeritas -{ -//---------------------------------------------------------------------------// -/*! - * Data for the Bragg and ICRU73QO ionization models. - * - * The ICRU73QO and Bragg models apply to negatively and positively charged - * incident particles, respectively. Both models apply to muons and hadrons. - * Because the sampling of the secondary and final state is identical in the - * two models, the same data and interactor are used for both. - */ -struct BraggICRU73QOData -{ - //! Particle IDs - ParticleId inc_particle; //!< Model-dependent incident particle - ParticleId electron; - - //! Electron mass [MeV / c^2] - units::MevMass electron_mass; - - //! Model-dependent kinetic energy limit [MeV] - units::MevEnergy lowest_kin_energy; - - //! Whether all data are assigned and valid - explicit CELER_FUNCTION operator bool() const - { - return inc_particle && electron && electron_mass > zero_quantity() - && lowest_kin_energy > zero_quantity(); - } -}; - -//---------------------------------------------------------------------------// -} // namespace celeritas diff --git a/src/celeritas/em/data/MuBetheBlochData.hh b/src/celeritas/em/data/MuHadIonizationData.hh similarity index 72% rename from src/celeritas/em/data/MuBetheBlochData.hh rename to src/celeritas/em/data/MuHadIonizationData.hh index 4637d68453..f7e5b962ae 100644 --- a/src/celeritas/em/data/MuBetheBlochData.hh +++ b/src/celeritas/em/data/MuHadIonizationData.hh @@ -3,11 +3,10 @@ // See the top-level COPYRIGHT file for details. // SPDX-License-Identifier: (Apache-2.0 OR MIT) //---------------------------------------------------------------------------// -//! \file celeritas/em/data/MuBetheBlochData.hh +//! \file celeritas/em/data/MuHadIonizationData.hh //---------------------------------------------------------------------------// #pragma once -#include "corecel/Macros.hh" #include "corecel/Types.hh" #include "celeritas/Quantities.hh" #include "celeritas/Types.hh" @@ -16,23 +15,22 @@ namespace celeritas { //---------------------------------------------------------------------------// /*! - * Device data for creating an interactor. + * Data for muon and hadron ionization. + * + * This data is used for the Bragg, ICRU73QO, Bethe-Bloch, and muon Bethe-Bloch + * models and can be reused for different incident particle types. */ -struct MuBetheBlochData +struct MuHadIonizationData { - //! Particle IDs + //! Secondary particle ID ParticleId electron; - ParticleId mu_minus; - ParticleId mu_plus; - //! Electron mass [MeV / c^2] units::MevMass electron_mass; //! Whether all data are assigned and valid explicit CELER_FUNCTION operator bool() const { - return electron && mu_minus && mu_plus - && electron_mass > zero_quantity(); + return electron && electron_mass > zero_quantity(); } }; diff --git a/src/celeritas/em/distribution/BetheBlochEnergyDistribution.hh b/src/celeritas/em/distribution/BetheBlochEnergyDistribution.hh new file mode 100644 index 0000000000..423845b8e0 --- /dev/null +++ b/src/celeritas/em/distribution/BetheBlochEnergyDistribution.hh @@ -0,0 +1,132 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2024 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/em/distribution/BetheBlochEnergyDistribution.hh +//---------------------------------------------------------------------------// +#pragma once + +#include + +#include "corecel/Macros.hh" +#include "corecel/Types.hh" +#include "corecel/math/Algorithms.hh" +#include "celeritas/Constants.hh" +#include "celeritas/Quantities.hh" +#include "celeritas/phys/ParticleTrackView.hh" +#include "celeritas/random/distribution/InverseSquareDistribution.hh" +#include "celeritas/random/distribution/RejectionSampler.hh" + +#include "detail/Utils.hh" + +namespace celeritas +{ +//---------------------------------------------------------------------------// +/*! + * Sample the energy of the delta ray for muon or hadron ionization. + * + * This samples the energy according to the Bethe-Bloch model, as described in + * the Geant4 Physics Reference Manual release 11.2 section 12.1.5. The + * Bethe-Bloch differential cross section can be written as + * \f[ + \difd{\sigma}{T} = 2\pi r_e^2 mc^2 Z \frac{z_p^2}{\beta^2}\frac{1}{T^2} + \left[1 - \beta^2 \frac{T}{T_{max}} + s \frac{T^2}{2E^2} \right] + * \f] + * and factorized as + * \f[ + \difd{\sigma}{T} = C f(T) g(T) + * \f] + * with \f$ T \in [T_{cut}, T_{max}] \f$, where \f$ f(T) = \frac{1}{T^2} \f$, + * \f$ g(T) = 1 - \beta^2 \frac{T}{T_max} + s \frac{T^2}{2 E^2} \f$, \f$ T \f$ + * is the kinetic energy of the electron, \f$ E \f$ is the total energy of the + * incident particle, and \f$ s \f$ is 0 for spinless particles and 1 + * otherwise. The energy is sampled from \f$ f(T) \f$ and accepted with + * probability \f$ g(T) \f$. + */ +class BetheBlochEnergyDistribution +{ + public: + //!@{ + //! \name Type aliases + using Energy = units::MevEnergy; + using Mass = units::MevMass; + //!@} + + public: + // Construct with incident and exiting particle data + inline CELER_FUNCTION + BetheBlochEnergyDistribution(ParticleTrackView const& particle, + Energy electron_cutoff, + Mass electron_mass); + + // Sample the exiting energy + template + inline CELER_FUNCTION Energy operator()(Engine& rng); + + //! Minimum energy of the secondary electron [MeV]. + CELER_FUNCTION Energy min_secondary_energy() const { return min_energy_; } + + //! Maximum energy of the secondary electron [MeV]. + CELER_FUNCTION Energy max_secondary_energy() const { return max_energy_; } + + private: + // Incident partcle mass + real_type inc_mass_; + // Square of fractional speed of light for incident particle + real_type beta_sq_; + // Minimum energy of the secondary electron [MeV] + Energy min_energy_; + // Maximum energy of the secondary electron [MeV] + Energy max_energy_; +}; + +//---------------------------------------------------------------------------// +// INLINE DEFINITIONS +//---------------------------------------------------------------------------// +/*! + * Construct with incident and exiting particle data. + */ +CELER_FUNCTION +BetheBlochEnergyDistribution::BetheBlochEnergyDistribution( + ParticleTrackView const& particle, + Energy electron_cutoff, + Mass electron_mass) + : inc_mass_(value_as(particle.mass())) + , beta_sq_(particle.beta_sq()) + , min_energy_(electron_cutoff) + , max_energy_(detail::calc_max_secondary_energy(particle, electron_mass)) +{ +} + +//---------------------------------------------------------------------------// +/*! + * Sample secondary electron energy. + */ +template +CELER_FUNCTION auto BetheBlochEnergyDistribution::operator()(Engine& rng) + -> Energy +{ + InverseSquareDistribution sample_energy(value_as(min_energy_), + value_as(max_energy_)); + real_type energy; + do + { + // Sample 1/E^2 from Emin to Emax + energy = sample_energy(rng); + /*! + * \todo Adjust rejection functions if particle has positive spin + */ + } while (RejectionSampler<>( + 1 - (beta_sq_ / value_as(max_energy_)) * energy)(rng)); + + /*! + * \todo For hadrons, suppress high energy delta ray production with the + * projectile form factor + */ + + return Energy{energy}; +} + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/celeritas/em/distribution/BraggICRU73QOEnergyDistribution.hh b/src/celeritas/em/distribution/BraggICRU73QOEnergyDistribution.hh new file mode 100644 index 0000000000..d4a190c2e5 --- /dev/null +++ b/src/celeritas/em/distribution/BraggICRU73QOEnergyDistribution.hh @@ -0,0 +1,134 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2024 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/em/distribution/BraggICRU73QOEnergyDistribution.hh +//---------------------------------------------------------------------------// +#pragma once + +#include + +#include "corecel/Macros.hh" +#include "corecel/Types.hh" +#include "corecel/math/Algorithms.hh" +#include "celeritas/Constants.hh" +#include "celeritas/Quantities.hh" +#include "celeritas/phys/ParticleTrackView.hh" +#include "celeritas/random/distribution/InverseSquareDistribution.hh" +#include "celeritas/random/distribution/RejectionSampler.hh" + +#include "detail/Utils.hh" + +namespace celeritas +{ +//---------------------------------------------------------------------------// +/*! + * Sample the energy of the delta ray for muon or hadron ionization. + * + * This samples the energy according to the Bragg and ICRU73QO models, as + * described in the Geant4 Physics Reference Manual release 11.2 section 11.1. + */ +class BraggICRU73QOEnergyDistribution +{ + public: + //!@{ + //! \name Type aliases + using Energy = units::MevEnergy; + using Mass = units::MevMass; + //!@} + + public: + // Construct with incident and exiting particle data + inline CELER_FUNCTION + BraggICRU73QOEnergyDistribution(ParticleTrackView const& particle, + Energy electron_cutoff, + Mass electron_mass); + + // Sample the exiting energy + template + inline CELER_FUNCTION Energy operator()(Engine& rng); + + //! Minimum energy of the secondary electron [MeV]. + CELER_FUNCTION Energy min_secondary_energy() const { return min_energy_; } + + //! Maximum energy of the secondary electron [MeV]. + CELER_FUNCTION Energy max_secondary_energy() const { return max_energy_; } + + private: + //// DATA //// + + // Incident partcle mass + real_type inc_mass_; + // Square of fractional speed of light for incident particle + real_type beta_sq_; + // Minimum energy of the secondary electron [MeV] + Energy min_energy_; + // Maximum energy of the secondary electron [MeV] + Energy max_energy_; + + //// CONSTANTS //// + + //! Used in Bragg model to calculate minimum energy transfer to electron + static CELER_CONSTEXPR_FUNCTION Energy bragg_lowest_kin_energy() + { + return units::MevEnergy{2.5e-4}; //! 0.25 keV + } + + //! Used in ICRU73QO model to calculate minimum energy transfer to electron + static CELER_CONSTEXPR_FUNCTION Energy icru73qo_lowest_kin_energy() + { + return units::MevEnergy{5e-3}; //! 5 keV + } +}; + +//---------------------------------------------------------------------------// +// INLINE DEFINITIONS +//---------------------------------------------------------------------------// +/*! + * Construct with incident and exiting particle data. + * + * \todo Use proton mass from imported data instead of a constant + */ +CELER_FUNCTION +BraggICRU73QOEnergyDistribution::BraggICRU73QOEnergyDistribution( + ParticleTrackView const& particle, + Energy electron_cutoff, + Mass electron_mass) + : inc_mass_(value_as(particle.mass())) + , beta_sq_(particle.beta_sq()) + , min_energy_( + min(electron_cutoff, + Energy(value_as(particle.charge() < zero_quantity() + ? icru73qo_lowest_kin_energy() + : bragg_lowest_kin_energy()) + * inc_mass_ + / native_value_to(constants::proton_mass).value()))) + , max_energy_(detail::calc_max_secondary_energy(particle, electron_mass)) +{ +} + +//---------------------------------------------------------------------------// +/*! + * Sample secondary electron energy. + */ +template +CELER_FUNCTION auto BraggICRU73QOEnergyDistribution::operator()(Engine& rng) + -> Energy +{ + InverseSquareDistribution sample_energy(value_as(min_energy_), + value_as(max_energy_)); + real_type energy; + do + { + // Sample 1/E^2 from Emin to Emax + energy = sample_energy(rng); + } while (RejectionSampler<>( + 1 - (beta_sq_ / value_as(max_energy_)) * energy)(rng)); + + CELER_ENSURE(energy > 0); + return Energy{energy}; +} + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/celeritas/em/distribution/MuBBEnergyDistribution.hh b/src/celeritas/em/distribution/MuBBEnergyDistribution.hh index b897784629..f8ecd9c3cb 100644 --- a/src/celeritas/em/distribution/MuBBEnergyDistribution.hh +++ b/src/celeritas/em/distribution/MuBBEnergyDistribution.hh @@ -13,14 +13,38 @@ #include "corecel/Types.hh" #include "corecel/math/Algorithms.hh" #include "celeritas/Quantities.hh" +#include "celeritas/phys/ParticleTrackView.hh" #include "celeritas/random/distribution/InverseSquareDistribution.hh" #include "celeritas/random/distribution/RejectionSampler.hh" +#include "detail/Utils.hh" + namespace celeritas { //---------------------------------------------------------------------------// /*! - * Sample the energy of the delta ray for the MuBetheBloch ionization model. + * Sample delta ray energy for the muon Bethe-Bloch ionization model. + 8 + * This samples the energy according to the muon Bethe-Bloch model, as + * described in the Geant4 Physics Reference Manual release 11.2 section 11.1. + * At the higher energies for which this model is applied, leading radiative + * corrections are taken into account. The differential cross section can be + * written as + * \f[ + \sigma(E, \epsilon) = \sigma_{BB}(E, \epsilon)\left[1 + \frac{\alpha}{2\pi} + \log \left(1 + \frac{2\epsilon}{m_e} \log \left(\frac{4 m_e E(E - + \epsilon}{m_{\mu}^2(2\epsilon + m_e)} \right) \right) \right]. + * \f] + * \f$ \sigma_{BB}(E, \epsilon) \f$ is the Bethe-Bloch cross section, \f$ m_e + * \f$ is the electron mass, \f$ m_{\mu} \f$ is the muon mass, \f$ E \f$ is the + * total energy of the muon, and \f$ \epsilon = \omega + T \f$ is the energy + * transfer, where \f$ T \f$ is the kinetic energy of the electron and \f$ + * \omega \f$ is the energy of the radiative gamma (which is neglected). + * + * As in the Bethe-Bloch model, the energy is sampled by factorizing the cross + * section as \f$ \sigma = C f(T) g(T) \f$, where \f$ f(T) = \frac{1}{T^2} \f$ + * and \f$ T \in [T_{cut}, T_{max}] \f$. The energy is sampled from \f$ f(T) + * \f$ and accepted with probability \f$ g(T) \f$. */ class MuBBEnergyDistribution { @@ -32,23 +56,25 @@ class MuBBEnergyDistribution //!@} public: - // Construct with data from MuBetheBlochInteractor - inline CELER_FUNCTION MuBBEnergyDistribution(Energy inc_energy, - Mass inc_mass, - real_type beta_sq, - Mass electron_mass, - Energy electron_cutoff, - Energy max_secondary_energy); + // Construct with incident and exiting particle data + inline CELER_FUNCTION + MuBBEnergyDistribution(ParticleTrackView const& particle, + Energy electron_cutoff, + Mass electron_mass); // Sample the exiting energy template inline CELER_FUNCTION Energy operator()(Engine& rng); + //! Minimum energy of the secondary electron [MeV]. + CELER_FUNCTION Energy min_secondary_energy() const { return min_energy_; } + + //! Maximum energy of the secondary electron [MeV]. + CELER_FUNCTION Energy max_secondary_energy() const { return max_energy_; } + private: //// DATA //// - // Electron incident energy [MeV] - real_type inc_energy_; // Incident partcle mass real_type inc_mass_; // Total energy of the incident particle [MeV] @@ -58,9 +84,9 @@ class MuBBEnergyDistribution // Secondary electron mass real_type electron_mass_; // Secondary electron cutoff energy [MeV] - real_type electron_cutoff_; + Energy min_energy_; // Maximum energy of the secondary electron [MeV] - real_type max_secondary_energy_; + Energy max_energy_; // Whether to apply the radiative correction bool use_rad_correction_; // Envelope distribution for rejection sampling @@ -100,22 +126,17 @@ class MuBBEnergyDistribution * Construct with incident and exiting particle data. */ CELER_FUNCTION -MuBBEnergyDistribution::MuBBEnergyDistribution(Energy inc_energy, - Mass inc_mass, - real_type beta_sq, - Mass electron_mass, +MuBBEnergyDistribution::MuBBEnergyDistribution(ParticleTrackView const& particle, Energy electron_cutoff, - Energy max_secondary_energy) - : inc_energy_(value_as(inc_energy)) - , inc_mass_(value_as(inc_mass)) - , total_energy_(inc_energy_ + inc_mass_) - , beta_sq_(beta_sq) + Mass electron_mass) + : inc_mass_(value_as(particle.mass())) + , total_energy_(value_as(particle.energy()) + inc_mass_) + , beta_sq_(particle.beta_sq()) , electron_mass_(value_as(electron_mass)) - , electron_cutoff_(value_as(electron_cutoff)) - , max_secondary_energy_(value_as(max_secondary_energy)) - , use_rad_correction_( - inc_energy_ > value_as(rad_correction_limit()) - && max_secondary_energy_ > value_as(kin_energy_limit())) + , min_energy_(electron_cutoff) + , max_energy_(detail::calc_max_secondary_energy(particle, electron_mass)) + , use_rad_correction_(particle.energy() > rad_correction_limit() + && max_energy_ > kin_energy_limit()) , envelope_(this->calc_envelope_distribution()) { } @@ -127,16 +148,14 @@ MuBBEnergyDistribution::MuBBEnergyDistribution(Energy inc_energy, template CELER_FUNCTION auto MuBBEnergyDistribution::operator()(Engine& rng) -> Energy { - // Sample from electron cutoff to max secondary energy - InverseSquareDistribution sample_energy{electron_cutoff_, - max_secondary_energy_}; - + InverseSquareDistribution sample_energy(value_as(min_energy_), + value_as(max_energy_)); real_type energy; real_type target; do { energy = sample_energy(rng); - target = 1 - (beta_sq_ / max_secondary_energy_) * energy + target = 1 - (beta_sq_ / value_as(max_energy_)) * energy + real_type(0.5) * ipow<2>(energy / total_energy_); if (use_rad_correction_ @@ -149,6 +168,7 @@ CELER_FUNCTION auto MuBBEnergyDistribution::operator()(Engine& rng) -> Energy } } while (RejectionSampler<>(target, envelope_)(rng)); + CELER_ENSURE(energy > 0); return Energy{energy}; } diff --git a/src/celeritas/em/distribution/detail/Utils.hh b/src/celeritas/em/distribution/detail/Utils.hh new file mode 100644 index 0000000000..f2357c42a4 --- /dev/null +++ b/src/celeritas/em/distribution/detail/Utils.hh @@ -0,0 +1,45 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2024 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/em/distribution/detail/Utils.hh +//---------------------------------------------------------------------------// +#pragma once + +#include "corecel/Types.hh" +#include "celeritas/Constants.hh" +#include "celeritas/Quantities.hh" +#include "celeritas/phys/ParticleTrackView.hh" + +namespace celeritas +{ +namespace detail +{ +//---------------------------------------------------------------------------// +/*! + * Calculate the maximum energy transferable to a free electron in ionizaation. + * + * This calculates the maximum kinematically allowed kinetic energy of the + * delta ray produced in muon or hadron ionization, + * \f[ + T_{max} = \frac{2 m_e c^2 (\gamma^2 - 1)}{1 + 2\gamma (m_e/M) + (m_e/M)^2}, + * \f] + * where \f$ m_e \f$ is the electron mass and \f$ M \f$ is the mass of the + * incident particle. + */ +inline CELER_FUNCTION units::MevEnergy +calc_max_secondary_energy(ParticleTrackView const& particle, + units::MevMass electron_mass) +{ + real_type inc_mass = value_as(particle.mass()); + real_type mass_ratio = value_as(electron_mass) / inc_mass; + real_type tau = value_as(particle.energy()) / inc_mass; + return units::MevEnergy{ + 2 * value_as(electron_mass) * tau * (tau + 2) + / (1 + 2 * (tau + 1) * mass_ratio + ipow<2>(mass_ratio))}; +} + +//---------------------------------------------------------------------------// +} // namespace detail +} // namespace celeritas diff --git a/src/celeritas/em/executor/MuBetheBlochExecutor.hh b/src/celeritas/em/executor/MuBetheBlochExecutor.hh deleted file mode 100644 index 5d8ee564ce..0000000000 --- a/src/celeritas/em/executor/MuBetheBlochExecutor.hh +++ /dev/null @@ -1,48 +0,0 @@ -//----------------------------------*-C++-*----------------------------------// -// Copyright 2024 UT-Battelle, LLC, and other Celeritas developers. -// See the top-level COPYRIGHT file for details. -// SPDX-License-Identifier: (Apache-2.0 OR MIT) -//---------------------------------------------------------------------------// -//! \file celeritas/em/executor/MuBetheBlochExecutor.hh -//---------------------------------------------------------------------------// -#pragma once - -#include "corecel/Assert.hh" -#include "corecel/Macros.hh" -#include "celeritas/em/data/MuBetheBlochData.hh" -#include "celeritas/em/interactor/MuBetheBlochInteractor.hh" -#include "celeritas/global/CoreTrackView.hh" -#include "celeritas/phys/Interaction.hh" - -namespace celeritas -{ -//---------------------------------------------------------------------------// -struct MuBetheBlochExecutor -{ - inline CELER_FUNCTION Interaction - operator()(celeritas::CoreTrackView const& track); - - MuBetheBlochData params; -}; - -//---------------------------------------------------------------------------// -/*! - * Apply the MuBetheBlochInteractor to the current track. - */ -CELER_FUNCTION Interaction -MuBetheBlochExecutor::operator()(CoreTrackView const& track) -{ - auto particle = track.make_particle_view(); - auto cutoff = track.make_cutoff_view(); - auto const& dir = track.make_geo_view().dir(); - auto allocate_secondaries - = track.make_physics_step_view().make_secondary_allocator(); - - MuBetheBlochInteractor interact( - params, particle, cutoff, dir, allocate_secondaries); - auto rng = track.make_rng_engine(); - return interact(rng); -} - -//---------------------------------------------------------------------------// -} // namespace celeritas diff --git a/src/celeritas/em/executor/BraggICRU73QOExecutor.hh b/src/celeritas/em/executor/MuHadIonizationExecutor.hh similarity index 74% rename from src/celeritas/em/executor/BraggICRU73QOExecutor.hh rename to src/celeritas/em/executor/MuHadIonizationExecutor.hh index ef5193a745..602dc9855e 100644 --- a/src/celeritas/em/executor/BraggICRU73QOExecutor.hh +++ b/src/celeritas/em/executor/MuHadIonizationExecutor.hh @@ -3,34 +3,36 @@ // See the top-level COPYRIGHT file for details. // SPDX-License-Identifier: (Apache-2.0 OR MIT) //---------------------------------------------------------------------------// -//! \file celeritas/em/executor/BraggICRU73QOExecutor.hh +//! \file celeritas/em/executor/MuHadIonizationExecutor.hh //---------------------------------------------------------------------------// #pragma once #include "corecel/Assert.hh" #include "corecel/Macros.hh" -#include "celeritas/em/data/BraggICRU73QOData.hh" -#include "celeritas/em/interactor/BraggICRU73QOInteractor.hh" +#include "celeritas/em/data/MuHadIonizationData.hh" +#include "celeritas/em/interactor/MuHadIonizationInteractor.hh" #include "celeritas/global/CoreTrackView.hh" #include "celeritas/phys/Interaction.hh" namespace celeritas { //---------------------------------------------------------------------------// -struct BraggICRU73QOExecutor +template +struct MuHadIonizationExecutor { inline CELER_FUNCTION Interaction operator()(celeritas::CoreTrackView const& track); - BraggICRU73QOData params; + MuHadIonizationData params; }; //---------------------------------------------------------------------------// /*! - * Apply the BraggICRU73QOInteractor to the current track. + * Apply the muon or hadron ionization interactor to the current track. */ +template CELER_FUNCTION Interaction -BraggICRU73QOExecutor::operator()(CoreTrackView const& track) +MuHadIonizationExecutor::operator()(CoreTrackView const& track) { auto particle = track.make_particle_view(); auto cutoff = track.make_cutoff_view(); @@ -38,7 +40,7 @@ BraggICRU73QOExecutor::operator()(CoreTrackView const& track) auto allocate_secondaries = track.make_physics_step_view().make_secondary_allocator(); - BraggICRU73QOInteractor interact( + MuHadIonizationInteractor interact( params, particle, cutoff, dir, allocate_secondaries); auto rng = track.make_rng_engine(); return interact(rng); diff --git a/src/celeritas/em/interactor/MuBetheBlochInteractor.hh b/src/celeritas/em/interactor/MuBetheBlochInteractor.hh deleted file mode 100644 index 5f5360b5fc..0000000000 --- a/src/celeritas/em/interactor/MuBetheBlochInteractor.hh +++ /dev/null @@ -1,179 +0,0 @@ -//----------------------------------*-C++-*----------------------------------// -// Copyright 2024 UT-Battelle, LLC, and other Celeritas developers. -// See the top-level COPYRIGHT file for details. -// SPDX-License-Identifier: (Apache-2.0 OR MIT) -//---------------------------------------------------------------------------// -//! \file celeritas/em/interactor/MuBetheBlochInteractor.hh -//---------------------------------------------------------------------------// -#pragma once - -#include "corecel/Macros.hh" -#include "corecel/Types.hh" -#include "corecel/data/StackAllocator.hh" -#include "corecel/math/Algorithms.hh" -#include "corecel/math/ArrayUtils.hh" -#include "celeritas/Constants.hh" -#include "celeritas/Quantities.hh" -#include "celeritas/em/data/MuBetheBlochData.hh" -#include "celeritas/em/distribution/MuBBEnergyDistribution.hh" -#include "celeritas/phys/CutoffView.hh" -#include "celeritas/phys/Interaction.hh" -#include "celeritas/phys/ParticleTrackView.hh" -#include "celeritas/phys/Secondary.hh" - -#include "detail/IoniFinalStateHelper.hh" -#include "detail/PhysicsConstants.hh" - -namespace celeritas -{ -//---------------------------------------------------------------------------// -/*! - * Perform the discrete part of the muon ionization process. - * - * This simulates the production of delta rays by incident mu- and mu+ with - * energies above 200 keV. - * - * \note This performs the same sampling routine as in Geant4's - * G4MuBetheBlochModel and as documented in the Geant4 Physics Reference Manual - * (Release 11.1) section 11.1. - */ -class MuBetheBlochInteractor -{ - public: - //!@{ - //! \name Type aliases - using Mass = units::MevMass; - using Energy = units::MevEnergy; - using Momentum = units::MevMomentum; - //!@} - - public: - // Construct with shared and state data - inline CELER_FUNCTION - MuBetheBlochInteractor(MuBetheBlochData const& shared, - ParticleTrackView const& particle, - CutoffView const& cutoffs, - Real3 const& inc_direction, - StackAllocator& allocate); - - // Sample an interaction with the given RNG - template - inline CELER_FUNCTION Interaction operator()(Engine& rng); - - private: - //// DATA //// - - // Allocate space for the secondary particle - StackAllocator& allocate_; - // Incident direction - Real3 const& inc_direction_; - // Incident particle energy [MeV] - Energy inc_energy_; - // Incident particle momentum [MeV / c] - Momentum inc_momentum_; - // Muon mass - Mass inc_mass_; - // Electron mass - Mass electron_mass_; - // Secondary electron ID - ParticleId electron_id_; - // Secondary electron cutoff for current material [MeV] - Energy electron_cutoff_; - // Maximum energy of the secondary electron [MeV] - Energy max_secondary_energy_; - // Secondary electron energy distribution - MuBBEnergyDistribution sample_energy_; - - //// HELPER FUNCTIONS //// - - inline CELER_FUNCTION Energy calc_max_secondary_energy() const; -}; - -//---------------------------------------------------------------------------// -// INLINE DEFINITIONS -//---------------------------------------------------------------------------// -/*! - * Construct with shared and state data. - */ -CELER_FUNCTION MuBetheBlochInteractor::MuBetheBlochInteractor( - MuBetheBlochData const& shared, - ParticleTrackView const& particle, - CutoffView const& cutoffs, - Real3 const& inc_direction, - StackAllocator& allocate) - : allocate_(allocate) - , inc_direction_(inc_direction) - , inc_energy_(particle.energy()) - , inc_momentum_(particle.momentum()) - , inc_mass_(particle.mass()) - , electron_mass_(shared.electron_mass) - , electron_id_(shared.electron) - , electron_cutoff_(cutoffs.energy(electron_id_)) - , max_secondary_energy_(this->calc_max_secondary_energy()) - , sample_energy_(inc_energy_, - inc_mass_, - particle.beta_sq(), - electron_mass_, - electron_cutoff_, - max_secondary_energy_) -{ - CELER_EXPECT(particle.particle_id() == shared.mu_minus - || particle.particle_id() == shared.mu_plus); - CELER_EXPECT(inc_energy_ > electron_cutoff_); - CELER_EXPECT(inc_energy_ >= detail::mu_bethe_bloch_lower_limit()); -} - -//---------------------------------------------------------------------------// -/*! - * Simulate discrete muon ionization. - */ -template -CELER_FUNCTION Interaction MuBetheBlochInteractor::operator()(Engine& rng) -{ - if (electron_cutoff_ > max_secondary_energy_) - { - // No interaction if the maximum secondary energy is below the cutoff - return Interaction::from_unchanged(); - } - - // Allocate secondary electron - Secondary* secondary = allocate_(1); - if (secondary == nullptr) - { - // Failed to allocate space for a secondary - return Interaction::from_failure(); - } - - // Sample the delta ray energy to construct the final sampler - detail::IoniFinalStateHelper sample_interaction(inc_energy_, - inc_direction_, - inc_momentum_, - inc_mass_, - sample_energy_(rng), - electron_mass_, - electron_id_, - secondary); - - // Update kinematics of the final state and return this interaction - return sample_interaction(rng); -} - -//---------------------------------------------------------------------------// -/*! - * Calculate maximum kinetic energy of the secondary electron. - * - * TODO this should be precomputed and used as part of the "applicability" per - * particle type, not included in the interactor. See celeritas#907 . - */ -CELER_FUNCTION auto -MuBetheBlochInteractor::calc_max_secondary_energy() const -> Energy -{ - real_type mass_ratio = value_as(electron_mass_) - / value_as(inc_mass_); - real_type tau = value_as(inc_energy_) / value_as(inc_mass_); - return Energy{2 * value_as(electron_mass_) * tau * (tau + 2) - / (1 + 2 * (tau + 1) * mass_ratio + ipow<2>(mass_ratio))}; -} - -//---------------------------------------------------------------------------// -} // namespace celeritas diff --git a/src/celeritas/em/interactor/BraggICRU73QOInteractor.hh b/src/celeritas/em/interactor/MuHadIonizationInteractor.hh similarity index 55% rename from src/celeritas/em/interactor/BraggICRU73QOInteractor.hh rename to src/celeritas/em/interactor/MuHadIonizationInteractor.hh index e9c2c44066..d29aad2aef 100644 --- a/src/celeritas/em/interactor/BraggICRU73QOInteractor.hh +++ b/src/celeritas/em/interactor/MuHadIonizationInteractor.hh @@ -3,7 +3,7 @@ // See the top-level COPYRIGHT file for details. // SPDX-License-Identifier: (Apache-2.0 OR MIT) //---------------------------------------------------------------------------// -//! \file celeritas/em/interactor/BraggICRU73QOInteractor.hh +//! \file celeritas/em/interactor/MuHadIonizationInteractor.hh //---------------------------------------------------------------------------// #pragma once @@ -14,7 +14,7 @@ #include "corecel/math/ArrayUtils.hh" #include "celeritas/Constants.hh" #include "celeritas/Quantities.hh" -#include "celeritas/em/data/BraggICRU73QOData.hh" +#include "celeritas/em/data/MuHadIonizationData.hh" #include "celeritas/phys/CutoffView.hh" #include "celeritas/phys/Interaction.hh" #include "celeritas/phys/ParticleTrackView.hh" @@ -31,18 +31,17 @@ namespace celeritas /*! * Perform the discrete part of the muon or hadron ionization process. * - * This simulates the production of delta rays by incident muons or hadrons - * in the low energy region according to the Bragg and ICRU73QO models. The - * ICRU73QO model applies to negatively charged incident particles, and the - * Bragg model applies to positively charged particles. The sampling of the - * secondary and the final state is identical in the two models, so a single - * interactor is used for both. + * This simulates the production of delta rays by incident muons or hadrons. + * The same basic sampling routine is used by multiple models, but the energy + * of the secondary is sampled from a distribution unique to the model. * - * \note This performs the same sampling routine as in Geant4's G4ICRU73QOModel - * and G4BraggModel and as documented in the Geant4 Physics Reference Manual - * (Release 11.1) section 12.2.1 and section 12.1. + * \note This performs the same sampling routine as in Geant4's \c + * G4BetheBlochModel, \c G4MuBetheBlochModel, \c G4BraggModel, and \c + * G4ICRU73QOModel, as documented in the Geant4 Physics Reference Manual + * release 11.2 sections 11.1 and 12.1.5. */ -class BraggICRU73QOInteractor +template +class MuHadIonizationInteractor { public: //!@{ @@ -55,11 +54,11 @@ class BraggICRU73QOInteractor public: //! Construct with shared and state data inline CELER_FUNCTION - BraggICRU73QOInteractor(BraggICRU73QOData const& shared, - ParticleTrackView const& particle, - CutoffView const& cutoffs, - Real3 const& inc_direction, - StackAllocator& allocate); + MuHadIonizationInteractor(MuHadIonizationData const& shared, + ParticleTrackView const& particle, + CutoffView const& cutoffs, + Real3 const& inc_direction, + StackAllocator& allocate); // Sample an interaction with the given RNG template @@ -78,20 +77,14 @@ class BraggICRU73QOInteractor Momentum inc_momentum_; // Muon mass [MeV / c^2] Mass inc_mass_; - // Square of fractional speed of light for incident particle - real_type beta_sq_; // Electron mass [MeV / c^2] Mass electron_mass_; // Secondary electron ID ParticleId electron_id_; - // Minimum energy of the secondary electron [MeV] - real_type min_secondary_energy_; // Maximum energy of the secondary electron [MeV] real_type max_secondary_energy_; - - //// HELPER FUNCTIONS //// - - inline CELER_FUNCTION Energy calc_max_secondary_energy() const; + // Secondary electron energy distribution + EnergySampler sample_energy_; }; //---------------------------------------------------------------------------// @@ -99,12 +92,10 @@ class BraggICRU73QOInteractor //---------------------------------------------------------------------------// /*! * Construct with shared and state data. - * - * \todo Use proton mass from imported data instead of a constant. */ -CELER_FUNCTION -BraggICRU73QOInteractor::BraggICRU73QOInteractor( - BraggICRU73QOData const& shared, +template +CELER_FUNCTION MuHadIonizationInteractor::MuHadIonizationInteractor( + MuHadIonizationData const& shared, ParticleTrackView const& particle, CutoffView const& cutoffs, Real3 const& inc_direction, @@ -114,28 +105,23 @@ BraggICRU73QOInteractor::BraggICRU73QOInteractor( , inc_energy_(particle.energy()) , inc_momentum_(particle.momentum()) , inc_mass_(particle.mass()) - , beta_sq_(particle.beta_sq()) , electron_mass_(shared.electron_mass) , electron_id_(shared.electron) - , min_secondary_energy_(min( - value_as(cutoffs.energy(electron_id_)), - value_as(shared.lowest_kin_energy) * value_as(inc_mass_) - / native_value_to(constants::proton_mass).value())) - , max_secondary_energy_(value_as(this->calc_max_secondary_energy())) + , sample_energy_(particle, cutoffs.energy(electron_id_), electron_mass_) { - CELER_EXPECT(particle.particle_id() == shared.inc_particle); - CELER_EXPECT(inc_energy_ > cutoffs.energy(electron_id_)); - CELER_EXPECT(inc_energy_ < detail::mu_bethe_bloch_lower_limit()); + CELER_EXPECT(inc_energy_ > sample_energy_.min_secondary_energy()); } //---------------------------------------------------------------------------// /*! * Simulate discrete muon ionization. */ +template template -CELER_FUNCTION Interaction BraggICRU73QOInteractor::operator()(Engine& rng) +CELER_FUNCTION Interaction MuHadIonizationInteractor::operator()(Engine& rng) { - if (min_secondary_energy_ >= max_secondary_energy_) + if (sample_energy_.min_secondary_energy() + >= sample_energy_.max_secondary_energy()) { // No interaction if the maximum secondary energy is below the limit return Interaction::from_unchanged(); @@ -149,44 +135,16 @@ CELER_FUNCTION Interaction BraggICRU73QOInteractor::operator()(Engine& rng) return Interaction::from_failure(); } - // Sample the delta ray energy - InverseSquareDistribution sample_energy{min_secondary_energy_, - max_secondary_energy_}; - - real_type energy; - do - { - // Sample 1/E^2 from Emin to Emax - energy = sample_energy(rng); - } while (RejectionSampler<>( - 1 - (beta_sq_ / max_secondary_energy_) * energy)(rng)); - // Update kinematics of the final state and return the interaction return detail::IoniFinalStateHelper(inc_energy_, inc_direction_, inc_momentum_, inc_mass_, - Energy{energy}, + sample_energy_(rng), electron_mass_, electron_id_, secondary)(rng); } -//---------------------------------------------------------------------------// -/*! - * Calculate maximum kinematically allowed kinetic energy of the secondary. - * - * TODO: Duplicated in \c MuBetheBlochInteractor. - */ -CELER_FUNCTION auto -BraggICRU73QOInteractor::calc_max_secondary_energy() const -> Energy -{ - real_type mass_ratio = value_as(electron_mass_) - / value_as(inc_mass_); - real_type tau = value_as(inc_energy_) / value_as(inc_mass_); - return Energy{2 * value_as(electron_mass_) * tau * (tau + 2) - / (1 + 2 * (tau + 1) * mass_ratio + ipow<2>(mass_ratio))}; -} - //---------------------------------------------------------------------------// } // namespace celeritas diff --git a/src/celeritas/em/interactor/detail/PhysicsConstants.hh b/src/celeritas/em/interactor/detail/PhysicsConstants.hh index 2dc780d591..9c2af57fcb 100644 --- a/src/celeritas/em/interactor/detail/PhysicsConstants.hh +++ b/src/celeritas/em/interactor/detail/PhysicsConstants.hh @@ -64,30 +64,12 @@ CELER_CONSTEXPR_FUNCTION units::MevEnergy seltzer_berger_upper_limit() return units::MevEnergy{1e3}; //! 1 GeV } -//! Energy threshold between muon ionization models. TODO: make configurable -CELER_CONSTEXPR_FUNCTION units::MevEnergy mu_bethe_bloch_lower_limit() -{ - return units::MevEnergy{0.2}; //! 200 keV -} - //! Maximum energy for EM models to be valid CELER_CONSTEXPR_FUNCTION units::MevEnergy high_energy_limit() { return units::MevEnergy{1e8}; //! 100 TeV } -//! Used in Bragg model to calculate minimum energy transfer to delta ray -CELER_CONSTEXPR_FUNCTION units::MevEnergy bragg_lowest_kin_energy() -{ - return units::MevEnergy{2.5e-4}; //! 0.25 keV -} - -//! Used in ICRU73QO model to calculate minimum energy transfer to delta ray -CELER_CONSTEXPR_FUNCTION units::MevEnergy icru73qo_lowest_kin_energy() -{ - return units::MevEnergy{5e-3}; //! 5 keV -} - //---------------------------------------------------------------------------// } // namespace detail } // namespace celeritas diff --git a/src/celeritas/em/model/BetheBlochModel.cc b/src/celeritas/em/model/BetheBlochModel.cc new file mode 100644 index 0000000000..b5a1e78365 --- /dev/null +++ b/src/celeritas/em/model/BetheBlochModel.cc @@ -0,0 +1,89 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2024 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/em/model/BetheBlochModel.cc +//---------------------------------------------------------------------------// +#include "BetheBlochModel.hh" + +#include "celeritas/Quantities.hh" +#include "celeritas/em/data/MuHadIonizationData.hh" +#include "celeritas/em/distribution/BetheBlochEnergyDistribution.hh" +#include "celeritas/em/executor/MuHadIonizationExecutor.hh" +#include "celeritas/em/interactor/detail/PhysicsConstants.hh" +#include "celeritas/global/ActionLauncher.hh" +#include "celeritas/global/CoreParams.hh" +#include "celeritas/global/CoreState.hh" +#include "celeritas/global/TrackExecutor.hh" +#include "celeritas/phys/InteractionApplier.hh" +#include "celeritas/phys/PDGNumber.hh" +#include "celeritas/phys/ParticleParams.hh" +#include "celeritas/phys/ParticleView.hh" + +#include "detail/MuHadIonizationBuilder.hh" + +namespace celeritas +{ +//---------------------------------------------------------------------------// +/*! + * Construct from model ID and other necessary data. + */ +BetheBlochModel::BetheBlochModel(ActionId id, + ParticleParams const& particles, + SetApplicability applicability) + : StaticConcreteAction( + id, "ioni-bethe-bloch", "interact by ionization (Bethe-Bloch)") + , applicability_(applicability) + , data_(detail::MuHadIonizationBuilder(particles, + this->label())(applicability_)) +{ + CELER_EXPECT(id); + CELER_ENSURE(data_); +} + +//---------------------------------------------------------------------------// +/*! + * Particle types and energy ranges that this model applies to. + */ +auto BetheBlochModel::applicability() const -> SetApplicability +{ + return applicability_; +} + +//---------------------------------------------------------------------------// +/*! + * Get the microscopic cross sections for the given particle and material. + */ +auto BetheBlochModel::micro_xs(Applicability) const -> MicroXsBuilders +{ + // Aside from the production cut, the discrete interaction is material + // independent, so no element is sampled + return {}; +} + +//---------------------------------------------------------------------------// +/*! + * Interact with host data. + */ +void BetheBlochModel::step(CoreParams const& params, CoreStateHost& state) const +{ + auto execute = make_action_track_executor( + params.ptr(), + state.ptr(), + this->action_id(), + InteractionApplier{MuHadIonizationExecutor{ + this->host_ref()}}); + return launch_action(*this, params, state, execute); +} + +//---------------------------------------------------------------------------// +#if !CELER_USE_DEVICE +void BetheBlochModel::step(CoreParams const&, CoreStateDevice&) const +{ + CELER_NOT_CONFIGURED("CUDA OR HIP"); +} +#endif + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/celeritas/em/model/BetheBlochModel.cu b/src/celeritas/em/model/BetheBlochModel.cu new file mode 100644 index 0000000000..8377399fde --- /dev/null +++ b/src/celeritas/em/model/BetheBlochModel.cu @@ -0,0 +1,38 @@ +//---------------------------------*-CUDA-*----------------------------------// +// Copyright 2024 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/em/model/BetheBlochModel.cu +//---------------------------------------------------------------------------// +#include "BetheBlochModel.hh" + +#include "celeritas/em/distribution/BetheBlochEnergyDistribution.hh" +#include "celeritas/em/executor/MuHadIonizationExecutor.hh" +#include "celeritas/global/ActionLauncher.device.hh" +#include "celeritas/global/CoreParams.hh" +#include "celeritas/global/CoreState.hh" +#include "celeritas/global/TrackExecutor.hh" +#include "celeritas/phys/InteractionApplier.hh" + +namespace celeritas +{ +//---------------------------------------------------------------------------// +/*! + * Interact with device data. + */ +void BetheBlochModel::step(CoreParams const& params, + CoreStateDevice& state) const +{ + auto execute = make_action_track_executor( + params.ptr(), + state.ptr(), + this->action_id(), + InteractionApplier{MuHadIonizationExecutor{ + this->device_ref()}}); + static ActionLauncher const launch_kernel(*this); + launch_kernel(*this, params, state, execute); +} + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/celeritas/em/model/BetheBlochModel.hh b/src/celeritas/em/model/BetheBlochModel.hh new file mode 100644 index 0000000000..1d0ffae21d --- /dev/null +++ b/src/celeritas/em/model/BetheBlochModel.hh @@ -0,0 +1,53 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2024 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/em/model/BetheBlochModel.hh +//---------------------------------------------------------------------------// +#pragma once + +#include "celeritas/em/data/MuHadIonizationData.hh" +#include "celeritas/phys/Model.hh" + +namespace celeritas +{ +class ParticleParams; + +//---------------------------------------------------------------------------// +/*! + * Set up and launch the Bethe-Bloch ionization model interaction. + */ +class BetheBlochModel final : public Model, public StaticConcreteAction +{ + public: + // Construct from model ID and other necessary data + BetheBlochModel(ActionId, ParticleParams const&, SetApplicability); + + // Particle types and energy ranges that this model applies to + SetApplicability applicability() const final; + + // Get the microscopic cross sections for the given particle and material + MicroXsBuilders micro_xs(Applicability) const final; + + // Apply the interaction kernel on host + void step(CoreParams const&, CoreStateHost&) const final; + + // Apply the interaction kernel on device + void step(CoreParams const&, CoreStateDevice&) const final; + + //!@{ + //! Access model data + MuHadIonizationData const& host_ref() const { return data_; } + MuHadIonizationData const& device_ref() const { return data_; } + //!@} + + private: + // Particle types and energy ranges that this model applies to + SetApplicability applicability_; + // Model data + MuHadIonizationData data_; +}; + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/celeritas/em/model/BraggModel.cc b/src/celeritas/em/model/BraggModel.cc index 2ec103578e..f324c5cebf 100644 --- a/src/celeritas/em/model/BraggModel.cc +++ b/src/celeritas/em/model/BraggModel.cc @@ -8,8 +8,8 @@ #include "BraggModel.hh" #include "celeritas/Quantities.hh" -#include "celeritas/em/data/BraggICRU73QOData.hh" -#include "celeritas/em/executor/BraggICRU73QOExecutor.hh" +#include "celeritas/em/distribution/BraggICRU73QOEnergyDistribution.hh" +#include "celeritas/em/executor/MuHadIonizationExecutor.hh" #include "celeritas/em/interactor/detail/PhysicsConstants.hh" #include "celeritas/global/ActionLauncher.hh" #include "celeritas/global/CoreParams.hh" @@ -20,30 +20,24 @@ #include "celeritas/phys/ParticleParams.hh" #include "celeritas/phys/ParticleView.hh" +#include "detail/MuHadIonizationBuilder.hh" + namespace celeritas { //---------------------------------------------------------------------------// /*! * Construct from model ID and other necessary data. - * - * TODO: This model also applies to hadrons. */ -BraggModel::BraggModel(ActionId id, ParticleParams const& particles) +BraggModel::BraggModel(ActionId id, + ParticleParams const& particles, + SetApplicability applicability) : StaticConcreteAction( - id, "ioni-bragg", "interact by muon ionization (Bragg)") + id, "ioni-bragg", "interact by muon ionization (Bragg)") + , applicability_(applicability) + , data_(detail::MuHadIonizationBuilder(particles, + this->label())(applicability_)) { CELER_EXPECT(id); - - data_.inc_particle = particles.find(pdg::mu_plus()); - data_.electron = particles.find(pdg::electron()); - - CELER_VALIDATE(data_.electron && data_.inc_particle, - << "missing electron and/or mu+ particles (required for " - << this->description() << ")"); - - data_.electron_mass = particles.get(data_.electron).mass(); - data_.lowest_kin_energy = detail::bragg_lowest_kin_energy(); - CELER_ENSURE(data_); } @@ -53,12 +47,7 @@ BraggModel::BraggModel(ActionId id, ParticleParams const& particles) */ auto BraggModel::applicability() const -> SetApplicability { - Applicability applic; - applic.particle = data_.inc_particle; - applic.lower = zero_quantity(); - applic.upper = detail::mu_bethe_bloch_lower_limit(); - - return {applic}; + return applicability_; } //---------------------------------------------------------------------------// @@ -82,7 +71,9 @@ void BraggModel::step(CoreParams const& params, CoreStateHost& state) const params.ptr(), state.ptr(), this->action_id(), - InteractionApplier{BraggICRU73QOExecutor{this->host_ref()}}); + InteractionApplier{ + MuHadIonizationExecutor{ + this->host_ref()}}); return launch_action(*this, params, state, execute); } diff --git a/src/celeritas/em/model/BraggModel.cu b/src/celeritas/em/model/BraggModel.cu index 1785ef745b..a84a328920 100644 --- a/src/celeritas/em/model/BraggModel.cu +++ b/src/celeritas/em/model/BraggModel.cu @@ -7,7 +7,8 @@ //---------------------------------------------------------------------------// #include "BraggModel.hh" -#include "celeritas/em/executor/BraggICRU73QOExecutor.hh" +#include "celeritas/em/distribution/BraggICRU73QOEnergyDistribution.hh" +#include "celeritas/em/executor/MuHadIonizationExecutor.hh" #include "celeritas/global/ActionLauncher.device.hh" #include "celeritas/global/CoreParams.hh" #include "celeritas/global/CoreState.hh" @@ -26,7 +27,9 @@ void BraggModel::step(CoreParams const& params, CoreStateDevice& state) const params.ptr(), state.ptr(), this->action_id(), - InteractionApplier{BraggICRU73QOExecutor{this->device_ref()}}); + InteractionApplier{ + MuHadIonizationExecutor{ + this->device_ref()}}); static ActionLauncher const launch_kernel(*this); launch_kernel(*this, params, state, execute); } diff --git a/src/celeritas/em/model/BraggModel.hh b/src/celeritas/em/model/BraggModel.hh index 704ce7a06e..b32b219a0d 100644 --- a/src/celeritas/em/model/BraggModel.hh +++ b/src/celeritas/em/model/BraggModel.hh @@ -7,7 +7,7 @@ //---------------------------------------------------------------------------// #pragma once -#include "celeritas/em/data/BraggICRU73QOData.hh" +#include "celeritas/em/data/MuHadIonizationData.hh" #include "celeritas/phys/Model.hh" namespace celeritas @@ -22,7 +22,7 @@ class BraggModel final : public Model, public StaticConcreteAction { public: // Construct from model ID and other necessary data - BraggModel(ActionId id, ParticleParams const& particles); + BraggModel(ActionId, ParticleParams const&, SetApplicability); // Particle types and energy ranges that this model applies to SetApplicability applicability() const final; @@ -38,12 +38,15 @@ class BraggModel final : public Model, public StaticConcreteAction //!@{ //! Access model data - BraggICRU73QOData const& host_ref() const { return data_; } - BraggICRU73QOData const& device_ref() const { return data_; } + MuHadIonizationData const& host_ref() const { return data_; } + MuHadIonizationData const& device_ref() const { return data_; } //!@} private: - BraggICRU73QOData data_; + // Particle types and energy ranges that this model applies to + SetApplicability applicability_; + // Model data + MuHadIonizationData data_; }; //---------------------------------------------------------------------------// diff --git a/src/celeritas/em/model/ICRU73QOModel.cc b/src/celeritas/em/model/ICRU73QOModel.cc index 07936280cc..5b9e00e311 100644 --- a/src/celeritas/em/model/ICRU73QOModel.cc +++ b/src/celeritas/em/model/ICRU73QOModel.cc @@ -8,8 +8,8 @@ #include "ICRU73QOModel.hh" #include "celeritas/Quantities.hh" -#include "celeritas/em/data/BraggICRU73QOData.hh" -#include "celeritas/em/executor/BraggICRU73QOExecutor.hh" +#include "celeritas/em/distribution/BraggICRU73QOEnergyDistribution.hh" +#include "celeritas/em/executor/MuHadIonizationExecutor.hh" #include "celeritas/em/interactor/detail/PhysicsConstants.hh" #include "celeritas/global/ActionLauncher.hh" #include "celeritas/global/CoreParams.hh" @@ -20,30 +20,24 @@ #include "celeritas/phys/ParticleParams.hh" #include "celeritas/phys/ParticleView.hh" +#include "detail/MuHadIonizationBuilder.hh" + namespace celeritas { //---------------------------------------------------------------------------// /*! * Construct from model ID and other necessary data. - * - * TODO: This model also applies to hadrons. */ -ICRU73QOModel::ICRU73QOModel(ActionId id, ParticleParams const& particles) +ICRU73QOModel::ICRU73QOModel(ActionId id, + ParticleParams const& particles, + SetApplicability applicability) : StaticConcreteAction( - id, "ioni-icru73qo", "interact by muon ionization (ICRU73QO)") + id, "ioni-icru73qo", "interact by muon ionization (ICRU73QO)") + , applicability_(applicability) + , data_(detail::MuHadIonizationBuilder(particles, + this->label())(applicability_)) { CELER_EXPECT(id); - - data_.inc_particle = particles.find(pdg::mu_minus()); - data_.electron = particles.find(pdg::electron()); - - CELER_VALIDATE(data_.electron && data_.inc_particle, - << "missing electron and/or mu- particles (required for " - << this->description() << ")"); - - data_.electron_mass = particles.get(data_.electron).mass(); - data_.lowest_kin_energy = detail::icru73qo_lowest_kin_energy(); - CELER_ENSURE(data_); } @@ -53,12 +47,7 @@ ICRU73QOModel::ICRU73QOModel(ActionId id, ParticleParams const& particles) */ auto ICRU73QOModel::applicability() const -> SetApplicability { - Applicability applic; - applic.particle = data_.inc_particle; - applic.lower = zero_quantity(); - applic.upper = detail::mu_bethe_bloch_lower_limit(); - - return {applic}; + return applicability_; } //---------------------------------------------------------------------------// @@ -82,7 +71,9 @@ void ICRU73QOModel::step(CoreParams const& params, CoreStateHost& state) const params.ptr(), state.ptr(), this->action_id(), - InteractionApplier{BraggICRU73QOExecutor{this->host_ref()}}); + InteractionApplier{ + MuHadIonizationExecutor{ + this->host_ref()}}); return launch_action(*this, params, state, execute); } diff --git a/src/celeritas/em/model/ICRU73QOModel.cu b/src/celeritas/em/model/ICRU73QOModel.cu index 8403974d0f..adfcedcbc5 100644 --- a/src/celeritas/em/model/ICRU73QOModel.cu +++ b/src/celeritas/em/model/ICRU73QOModel.cu @@ -7,7 +7,8 @@ //---------------------------------------------------------------------------// #include "ICRU73QOModel.hh" -#include "celeritas/em/executor/BraggICRU73QOExecutor.hh" +#include "celeritas/em/distribution/BraggICRU73QOEnergyDistribution.hh" +#include "celeritas/em/executor/MuHadIonizationExecutor.hh" #include "celeritas/global/ActionLauncher.device.hh" #include "celeritas/global/CoreParams.hh" #include "celeritas/global/CoreState.hh" @@ -26,7 +27,9 @@ void ICRU73QOModel::step(CoreParams const& params, CoreStateDevice& state) const params.ptr(), state.ptr(), this->action_id(), - InteractionApplier{BraggICRU73QOExecutor{this->device_ref()}}); + InteractionApplier{ + MuHadIonizationExecutor{ + this->device_ref()}}); static ActionLauncher const launch_kernel(*this); launch_kernel(*this, params, state, execute); } diff --git a/src/celeritas/em/model/ICRU73QOModel.hh b/src/celeritas/em/model/ICRU73QOModel.hh index 1ec3fce150..3994057513 100644 --- a/src/celeritas/em/model/ICRU73QOModel.hh +++ b/src/celeritas/em/model/ICRU73QOModel.hh @@ -7,7 +7,7 @@ //---------------------------------------------------------------------------// #pragma once -#include "celeritas/em/data/BraggICRU73QOData.hh" +#include "celeritas/em/data/MuHadIonizationData.hh" #include "celeritas/phys/Model.hh" namespace celeritas @@ -22,7 +22,7 @@ class ICRU73QOModel final : public Model, public StaticConcreteAction { public: // Construct from model ID and other necessary data - ICRU73QOModel(ActionId id, ParticleParams const& particles); + ICRU73QOModel(ActionId, ParticleParams const&, SetApplicability); // Particle types and energy ranges that this model applies to SetApplicability applicability() const final; @@ -38,12 +38,15 @@ class ICRU73QOModel final : public Model, public StaticConcreteAction //!@{ //! Access model data - BraggICRU73QOData const& host_ref() const { return data_; } - BraggICRU73QOData const& device_ref() const { return data_; } + MuHadIonizationData const& host_ref() const { return data_; } + MuHadIonizationData const& device_ref() const { return data_; } //!@} private: - BraggICRU73QOData data_; + // Particle types and energy ranges that this model applies to + SetApplicability applicability_; + // Model data + MuHadIonizationData data_; }; //---------------------------------------------------------------------------// diff --git a/src/celeritas/em/model/MuBetheBlochModel.cc b/src/celeritas/em/model/MuBetheBlochModel.cc index 1435c817d8..e23a52e325 100644 --- a/src/celeritas/em/model/MuBetheBlochModel.cc +++ b/src/celeritas/em/model/MuBetheBlochModel.cc @@ -8,8 +8,9 @@ #include "MuBetheBlochModel.hh" #include "celeritas/Quantities.hh" -#include "celeritas/em/data/MuBetheBlochData.hh" -#include "celeritas/em/executor/MuBetheBlochExecutor.hh" +#include "celeritas/em/data/MuHadIonizationData.hh" +#include "celeritas/em/distribution/MuBBEnergyDistribution.hh" +#include "celeritas/em/executor/MuHadIonizationExecutor.hh" #include "celeritas/em/interactor/detail/PhysicsConstants.hh" #include "celeritas/global/ActionLauncher.hh" #include "celeritas/global/CoreParams.hh" @@ -20,6 +21,8 @@ #include "celeritas/phys/ParticleParams.hh" #include "celeritas/phys/ParticleView.hh" +#include "detail/MuHadIonizationBuilder.hh" + namespace celeritas { //---------------------------------------------------------------------------// @@ -27,22 +30,15 @@ namespace celeritas * Construct from model ID and other necessary data. */ MuBetheBlochModel::MuBetheBlochModel(ActionId id, - ParticleParams const& particles) - : StaticConcreteAction(id, - "ioni-mu-bethe-bloch", - "interact by muon ionization (Bethe-Bloch)") + ParticleParams const& particles, + SetApplicability applicability) + : StaticConcreteAction( + id, "ioni-mu-bethe-bloch", "interact by muon ionization (Bethe-Bloch)") + , applicability_(applicability) + , data_(detail::MuHadIonizationBuilder(particles, + this->label())(applicability_)) { CELER_EXPECT(id); - data_.electron = particles.find(pdg::electron()); - data_.mu_minus = particles.find(pdg::mu_minus()); - data_.mu_plus = particles.find(pdg::mu_plus()); - - CELER_VALIDATE(data_.electron && data_.mu_minus && data_.mu_plus, - << "missing electron and/or muon particles (required for " - << this->description() << ")"); - - data_.electron_mass = particles.get(data_.electron).mass(); - CELER_ENSURE(data_); } @@ -52,15 +48,7 @@ MuBetheBlochModel::MuBetheBlochModel(ActionId id, */ auto MuBetheBlochModel::applicability() const -> SetApplicability { - Applicability mu_minus_applic; - mu_minus_applic.particle = data_.mu_minus; - mu_minus_applic.lower = detail::mu_bethe_bloch_lower_limit(); - mu_minus_applic.upper = detail::high_energy_limit(); - - Applicability mu_plus_applic = mu_minus_applic; - mu_plus_applic.particle = data_.mu_plus; - - return {mu_minus_applic, mu_plus_applic}; + return applicability_; } //---------------------------------------------------------------------------// @@ -85,7 +73,8 @@ void MuBetheBlochModel::step(CoreParams const& params, params.ptr(), state.ptr(), this->action_id(), - InteractionApplier{MuBetheBlochExecutor{this->host_ref()}}); + InteractionApplier{MuHadIonizationExecutor{ + this->host_ref()}}); return launch_action(*this, params, state, execute); } diff --git a/src/celeritas/em/model/MuBetheBlochModel.cu b/src/celeritas/em/model/MuBetheBlochModel.cu index a34be3abc1..ad4053a910 100644 --- a/src/celeritas/em/model/MuBetheBlochModel.cu +++ b/src/celeritas/em/model/MuBetheBlochModel.cu @@ -7,7 +7,8 @@ //---------------------------------------------------------------------------// #include "MuBetheBlochModel.hh" -#include "celeritas/em/executor/MuBetheBlochExecutor.hh" +#include "celeritas/em/distribution/MuBBEnergyDistribution.hh" +#include "celeritas/em/executor/MuHadIonizationExecutor.hh" #include "celeritas/global/ActionLauncher.device.hh" #include "celeritas/global/CoreParams.hh" #include "celeritas/global/CoreState.hh" @@ -27,7 +28,8 @@ void MuBetheBlochModel::step(CoreParams const& params, params.ptr(), state.ptr(), this->action_id(), - InteractionApplier{MuBetheBlochExecutor{this->device_ref()}}); + InteractionApplier{MuHadIonizationExecutor{ + this->device_ref()}}); static ActionLauncher const launch_kernel(*this); launch_kernel(*this, params, state, execute); } diff --git a/src/celeritas/em/model/MuBetheBlochModel.hh b/src/celeritas/em/model/MuBetheBlochModel.hh index 49634434a6..6c6759efa5 100644 --- a/src/celeritas/em/model/MuBetheBlochModel.hh +++ b/src/celeritas/em/model/MuBetheBlochModel.hh @@ -7,7 +7,7 @@ //---------------------------------------------------------------------------// #pragma once -#include "celeritas/em/data/MuBetheBlochData.hh" +#include "celeritas/em/data/MuHadIonizationData.hh" #include "celeritas/phys/Model.hh" namespace celeritas @@ -22,7 +22,7 @@ class MuBetheBlochModel final : public Model, public StaticConcreteAction { public: // Construct from model ID and other necessary data - MuBetheBlochModel(ActionId id, ParticleParams const& particles); + MuBetheBlochModel(ActionId, ParticleParams const&, SetApplicability); // Particle types and energy ranges that this model applies to SetApplicability applicability() const final; @@ -38,12 +38,15 @@ class MuBetheBlochModel final : public Model, public StaticConcreteAction //!@{ //! Access model data - MuBetheBlochData const& host_ref() const { return data_; } - MuBetheBlochData const& device_ref() const { return data_; } + MuHadIonizationData const& host_ref() const { return data_; } + MuHadIonizationData const& device_ref() const { return data_; } //!@} private: - MuBetheBlochData data_; + // Particle types and energy ranges that this model applies to + SetApplicability applicability_; + // Model data + MuHadIonizationData data_; }; //---------------------------------------------------------------------------// diff --git a/src/celeritas/em/model/detail/MuHadIonizationBuilder.hh b/src/celeritas/em/model/detail/MuHadIonizationBuilder.hh new file mode 100644 index 0000000000..c22fc5651c --- /dev/null +++ b/src/celeritas/em/model/detail/MuHadIonizationBuilder.hh @@ -0,0 +1,98 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2024 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/em/model/detail/MuHadIonizationBuilder.hh +//---------------------------------------------------------------------------// +#pragma once + +#include +#include + +#include "celeritas/Quantities.hh" +#include "celeritas/em/data/MuHadIonizationData.hh" +#include "celeritas/phys/Applicability.hh" +#include "celeritas/phys/ParticleParams.hh" + +namespace celeritas +{ +namespace detail +{ +//---------------------------------------------------------------------------// +/*! + * Construct muon and hadron ionization data. + * + * This small helper class constructs and validates the data for the muon and + * hadron ionization models. + */ +class MuHadIonizationBuilder +{ + public: + //!@{ + //! \name Type aliases + using Energy = units::MevEnergy; + using SetApplicability = std::set; + //!@} + + public: + // Construct with shared particle data and model label + inline MuHadIonizationBuilder(ParticleParams const& particles, + std::string_view label); + + // Construct model data from applicability + inline MuHadIonizationData operator()(SetApplicability const&) const; + + private: + ParticleParams const& particles_; + std::string_view label_; +}; + +//---------------------------------------------------------------------------// +// INLINE DEFINITIONS +//---------------------------------------------------------------------------// +/*! + * Construct with shared particle data and model label. + */ +MuHadIonizationBuilder::MuHadIonizationBuilder(ParticleParams const& particles, + std::string_view label) + : particles_(particles), label_(label) +{ +} + +//---------------------------------------------------------------------------// +/*! + * Construct model data from applicability. + */ +MuHadIonizationData +MuHadIonizationBuilder::operator()(SetApplicability const& applicability) const +{ + CELER_EXPECT(!applicability.empty()); + + MuHadIonizationData data; + + for (auto const& applic : applicability) + { + CELER_VALIDATE( + applic, + << "invalid applicability with particle `" + << (applic.particle ? particles_.id_to_label(applic.particle) : "") + << "' and energy limits (" << value_as(applic.lower) + << ", " << value_as(applic.upper) << ") [MeV] for model '" + << label_ << "'"); + } + + data.electron = particles_.find(pdg::electron()); + CELER_VALIDATE(data.electron, + << "missing electron (required for model '" << label_ + << "')"); + + data.electron_mass = particles_.get(data.electron).mass(); + + CELER_ENSURE(data); + return data; +} + +//---------------------------------------------------------------------------// +} // namespace detail +} // namespace celeritas diff --git a/src/celeritas/em/process/MuIonizationProcess.cc b/src/celeritas/em/process/MuIonizationProcess.cc index b906dd2e49..c796dcbcd1 100644 --- a/src/celeritas/em/process/MuIonizationProcess.cc +++ b/src/celeritas/em/process/MuIonizationProcess.cc @@ -11,6 +11,8 @@ #include "corecel/Assert.hh" #include "corecel/cont/Range.hh" +#include "celeritas/em/interactor/detail/PhysicsConstants.hh" +#include "celeritas/em/model/BetheBlochModel.hh" #include "celeritas/em/model/BraggModel.hh" #include "celeritas/em/model/ICRU73QOModel.hh" #include "celeritas/em/model/MuBetheBlochModel.hh" @@ -39,12 +41,55 @@ MuIonizationProcess::MuIonizationProcess(SPConstParticles particles, //---------------------------------------------------------------------------// /*! * Construct the models associated with this process. + * + * \todo Possibly get model limits from imported data? */ auto MuIonizationProcess::build_models(ActionIdIter start_id) const -> VecModel { - return {std::make_shared(*start_id++, *particles_), - std::make_shared(*start_id++, *particles_), - std::make_shared(*start_id++, *particles_)}; + using IMC = ImportModelClass; + using SetApplicability = Model::SetApplicability; + + VecModel result; + + // Construct ICRU73QO mu- model + CELER_ASSERT(imported_.has_model(pdg::mu_minus(), IMC::icru_73_qo)); + Applicability mm; + mm.particle = particles_->find(pdg::mu_minus()); + mm.lower = zero_quantity(); + mm.upper = options_.bragg_icru73qo_upper_limit; + result.push_back(std::make_shared( + *start_id++, *particles_, SetApplicability{mm})); + + // Construct Bragg mu+ model + CELER_ASSERT(imported_.has_model(pdg::mu_plus(), IMC::bragg)); + Applicability mp = mm; + mp.particle = particles_->find(pdg::mu_plus()); + result.push_back(std::make_shared( + *start_id++, *particles_, SetApplicability{mp})); + + if (imported_.has_model(pdg::mu_minus(), IMC::bethe_bloch)) + { + // Older Geant4 versions use Bethe-Bloch at intermediate energies + CELER_ASSERT(imported_.has_model(pdg::mu_plus(), IMC::bethe_bloch)); + mm.lower = mm.upper; + mm.upper = options_.bethe_bloch_upper_limit; + mp.lower = mm.lower; + mp.upper = mm.upper; + result.push_back(std::make_shared( + *start_id++, *particles_, SetApplicability{mm, mp})); + } + + // Construct muon Bethe-Bloch model + CELER_ASSERT(imported_.has_model(pdg::mu_minus(), IMC::mu_bethe_bloch)); + CELER_ASSERT(imported_.has_model(pdg::mu_plus(), IMC::mu_bethe_bloch)); + mm.lower = mm.upper; + mm.upper = detail::high_energy_limit(); + mp.lower = mm.lower; + mp.upper = mm.upper; + result.push_back(std::make_shared( + *start_id++, *particles_, SetApplicability{mm, mp})); + + return result; } //---------------------------------------------------------------------------// diff --git a/src/celeritas/em/process/MuIonizationProcess.hh b/src/celeritas/em/process/MuIonizationProcess.hh index b648f78556..a6941e7268 100644 --- a/src/celeritas/em/process/MuIonizationProcess.hh +++ b/src/celeritas/em/process/MuIonizationProcess.hh @@ -25,6 +25,7 @@ class MuIonizationProcess : public Process public: //!@{ //! \name Type aliases + using Energy = units::MevEnergy; using SPConstParticles = std::shared_ptr; using SPConstImported = std::shared_ptr; //!@} @@ -32,8 +33,12 @@ class MuIonizationProcess : public Process // Options for electron and positron ionization struct Options { - bool use_integral_xs{true}; //!> Use integral method for sampling - //! discrete interaction length + //! Maximum energy for the Bragg and ICRU73QO models + Energy bragg_icru73qo_upper_limit{0.2}; //!< 200 keV + //! Maximum energy for the Bethe-Bloch model + Energy bethe_bloch_upper_limit{1e3}; //!< 1 GeV + //! Use integral method for sampling discrete interaction length + bool use_integral_xs{true}; }; public: diff --git a/src/celeritas/ext/detail/CelerEmStandardPhysics.cc b/src/celeritas/ext/detail/CelerEmStandardPhysics.cc index 08e08ee983..e5850e4cfa 100644 --- a/src/celeritas/ext/detail/CelerEmStandardPhysics.cc +++ b/src/celeritas/ext/detail/CelerEmStandardPhysics.cc @@ -467,7 +467,7 @@ void CelerEmStandardPhysics::add_e_processes(G4ParticleDefinition* p) * * \note Currently all muon processes are disabled by default * - * \todo Prior to version 11.1.0, Geant4 used the \c G4BetheBlochModel for muon + * \note Prior to version 11.1.0, Geant4 used the \c G4BetheBlochModel for muon * ionization between 200 keV and 1 GeV and the \c G4MuBetheBlochModel above 1 * GeV. Since version 11.1.0, the \c G4MuBetheBlochModel is used for all * energies above 200 keV. diff --git a/src/celeritas/phys/ImportedProcessAdapter.hh b/src/celeritas/phys/ImportedProcessAdapter.hh index 7334bd4ef9..8e1dbabd2c 100644 --- a/src/celeritas/phys/ImportedProcessAdapter.hh +++ b/src/celeritas/phys/ImportedProcessAdapter.hh @@ -7,6 +7,7 @@ //---------------------------------------------------------------------------// #pragma once +#include #include #include #include @@ -122,6 +123,9 @@ class ImportedProcessAdapter // Access the imported processes SPConstImported const& processes() const { return imported_; } + // Whether the given model is present in the process + inline bool has_model(PDGNumber, ImportModelClass) const; + private: using ImportTableId = OpaqueId; using ImportProcessId = ImportedProcesses::ImportProcessId; @@ -180,5 +184,19 @@ ImportedProcessAdapter::get_lambda(ParticleId id) const return imported_->get(iter->second.process).tables[tab.unchecked_get()]; } +//---------------------------------------------------------------------------// +/*! + * Whether the given model is present in the process. + */ +bool ImportedProcessAdapter::has_model(PDGNumber pdg, ImportModelClass imc) const +{ + auto const& models + = imported_->get(imported_->find({pdg, process_class_})).models; + return std::any_of( + models.begin(), models.end(), [&imc](ImportModel const& m) { + return m.model_class == imc; + }); +} + //---------------------------------------------------------------------------// } // namespace celeritas diff --git a/src/orange/detail/OrientedBoundingZone.hh b/src/orange/detail/OrientedBoundingZone.hh index c4fe9a3bac..92265c90db 100644 --- a/src/orange/detail/OrientedBoundingZone.hh +++ b/src/orange/detail/OrientedBoundingZone.hh @@ -31,10 +31,10 @@ namespace detail * the center of each bounding box is offset from the center of the OBZ * coordinate system. * - * Consequently, points in the unit's corrdinate system must be first - * transformed by into the OBZ coordinate system (resulting in - * "trans_pos"), then offset, i.e. translated, by or - * into the inner or outer bbox coordinate system (resulting + * Consequently, points in the unit's coordinate system must be first + * transformed by \c transform_id into the OBZ coordinate system (resulting in + * "trans_pos"), then offset, i.e. translated, by \c inner_offset_id or + * \c outer_offset_id into the inner or outer bbox coordinate system (resulting * in "offset_pos"). It is noted that these offset positions are always * automatically reflected into the first quadrant. * diff --git a/test/celeritas/CMakeLists.txt b/test/celeritas/CMakeLists.txt index 6247f592ff..61de25c977 100644 --- a/test/celeritas/CMakeLists.txt +++ b/test/celeritas/CMakeLists.txt @@ -122,6 +122,7 @@ celeritas_add_test(em/TsaiUrbanDistribution.test.cc ${_needs_double}) # Models celeritas_add_test(em/BetheHeitler.test.cc ${_needs_double}) +celeritas_add_test(em/BetheBloch.test.cc ${_needs_double}) celeritas_add_test(em/BraggICRU73QO.test.cc ${_needs_double}) celeritas_add_test(em/CombinedBrem.test.cc ${_needs_double}) celeritas_add_test(em/EPlusGG.test.cc ${_needs_double}) diff --git a/test/celeritas/em/BetheBloch.test.cc b/test/celeritas/em/BetheBloch.test.cc new file mode 100644 index 0000000000..16dd934ad7 --- /dev/null +++ b/test/celeritas/em/BetheBloch.test.cc @@ -0,0 +1,334 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2024 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/em/BetheBloch.test.cc +//---------------------------------------------------------------------------// +#include "corecel/cont/Range.hh" +#include "corecel/math/ArrayUtils.hh" +#include "celeritas/Quantities.hh" +#include "celeritas/em/distribution/BetheBlochEnergyDistribution.hh" +#include "celeritas/em/interactor/MuHadIonizationInteractor.hh" +#include "celeritas/em/model/BetheBlochModel.hh" +#include "celeritas/em/process/MuIonizationProcess.hh" +#include "celeritas/phys/CutoffView.hh" +#include "celeritas/phys/InteractionIO.hh" +#include "celeritas/phys/InteractorHostTestBase.hh" + +#include "celeritas_test.hh" + +namespace celeritas +{ +namespace test +{ +//---------------------------------------------------------------------------// +// TEST HARNESS +//---------------------------------------------------------------------------// + +class BetheBlochTest : public InteractorHostTestBase +{ + using Base = InteractorHostTestBase; + + protected: + void SetUp() override + { + using namespace units; + + // Set up shared material data + MaterialParams::Input mat_inp; + mat_inp.elements = {{AtomicNumber{29}, AmuMass{63.546}, {}, "Cu"}}; + mat_inp.materials = { + {native_value_from(MolCcDensity{0.141}), + 293.0, + MatterState::solid, + {{ElementId{0}, 1.0}}, + "Cu"}, + }; + this->set_material_params(mat_inp); + + // Set 1 keV electron cutoff + CutoffParams::Input cut_inp; + cut_inp.materials = this->material_params(); + cut_inp.particles = this->particle_params(); + cut_inp.cutoffs = {{pdg::electron(), {{MevEnergy{0.001}, 0.1234}}}}; + this->set_cutoff_params(cut_inp); + + // Set model data + auto const& particles = *this->particle_params(); + Applicability mu_minus; + mu_minus.particle = particles.find(pdg::mu_minus()); + mu_minus.lower + = MuIonizationProcess::Options{}.bragg_icru73qo_upper_limit; + mu_minus.upper = MuIonizationProcess::Options{}.bethe_bloch_upper_limit; + Applicability mu_plus = mu_minus; + mu_plus.particle = particles.find(pdg::mu_plus()); + model_ = std::make_shared( + ActionId{0}, particles, Model::SetApplicability{mu_minus, mu_plus}); + + // Set default particle to muon with energy of 10 MeV + this->set_inc_particle(pdg::mu_minus(), MevEnergy{10}); + this->set_inc_direction({0, 0, 1}); + this->set_material("Cu"); + } + + void sanity_check(Interaction const& interaction) const + { + // Check change to parent track + EXPECT_GT(this->particle_track().energy().value(), + interaction.energy.value()); + EXPECT_LT(0, interaction.energy.value()); + EXPECT_SOFT_EQ(1.0, norm(interaction.direction)); + EXPECT_EQ(Action::scattered, interaction.action); + + // Check secondaries + ASSERT_EQ(1, interaction.secondaries.size()); + + auto const& electron = interaction.secondaries.front(); + EXPECT_TRUE(electron); + EXPECT_EQ(model_->host_ref().electron, electron.particle_id); + EXPECT_GT(this->particle_track().energy().value(), + electron.energy.value()); + EXPECT_LT(0, electron.energy.value()); + EXPECT_SOFT_EQ(1.0, norm(electron.direction)); + + // Check conservation between primary and secondaries + this->check_conservation(interaction); + this->check_energy_conservation(interaction); + } + + protected: + std::shared_ptr model_; +}; + +//---------------------------------------------------------------------------// +// TESTS +//---------------------------------------------------------------------------// + +TEST_F(BetheBlochTest, distribution) +{ + int num_samples = 100000; + int num_bins = 12; + + MevEnergy cutoff{0.001}; + + std::vector counters; + std::vector min_energy; + std::vector max_energy; + for (real_type energy : {0.2, 1.0, 10.0, 1e2, 1e3, 1e4, 1e5, 1e7}) + { + this->set_inc_particle(pdg::mu_minus(), MevEnergy(energy)); + std::mt19937 rng; + + BetheBlochEnergyDistribution sample( + this->particle_track(), cutoff, model_->host_ref().electron_mass); + real_type min = value_as(sample.min_secondary_energy()); + real_type max = value_as(sample.max_secondary_energy()); + + std::vector count(num_bins); + for ([[maybe_unused]] int i : range(num_samples)) + { + auto r = value_as(sample(rng)); + ASSERT_GE(r, min); + ASSERT_LE(r, max); + int bin = int((1 / r - 1 / min) / (1 / max - 1 / min) * num_bins); + CELER_ASSERT(bin >= 0 && bin < num_bins); + ++count[bin]; + } + counters.insert(counters.end(), count.begin(), count.end()); + min_energy.push_back(min); + max_energy.push_back(max); + } + + static int const expected_counters[] = { + 8273, 8487, 8379, 8203, 8303, 8419, 8249, 8422, 8366, 8265, 8334, 8300, + 8281, 8499, 8383, 8211, 8309, 8428, 8256, 8435, 8371, 8263, 8320, 8244, + 8294, 8514, 8391, 8225, 8326, 8440, 8269, 8446, 8391, 8279, 8321, 8104, + 8283, 8499, 8389, 8211, 8312, 8427, 8257, 8440, 8380, 8278, 8340, 8184, + 8268, 8482, 8377, 8202, 8302, 8416, 8240, 8429, 8371, 8267, 8338, 8308, + 8267, 8480, 8377, 8202, 8301, 8415, 8239, 8429, 8371, 8266, 8338, 8315, + 8267, 8480, 8377, 8202, 8301, 8415, 8239, 8429, 8370, 8266, 8338, 8316, + 8267, 8480, 8377, 8202, 8301, 8415, 8239, 8429, 8370, 8266, 8338, 8316, + }; + static double const expected_min_energy[] + = {0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001}; + static double const expected_max_energy[] = { + 0.0038354680957569, + 0.019248476995285, + 0.20048052363148, + 2.7972680400033, + 100.69707462436, + 4855.7535710157, + 90256.629501068, + 9989193.9209199, + }; + EXPECT_VEC_EQ(expected_counters, counters); + EXPECT_VEC_SOFT_EQ(expected_min_energy, min_energy); + EXPECT_VEC_SOFT_EQ(expected_max_energy, max_energy); +} + +TEST_F(BetheBlochTest, basic) +{ + // Reserve 4 secondaries, one for each sample + int const num_samples = 4; + this->resize_secondaries(num_samples); + + // Create the interactor + MuHadIonizationInteractor interact( + model_->host_ref(), + this->particle_track(), + this->cutoff_params()->get(MaterialId{0}), + this->direction(), + this->secondary_allocator()); + RandomEngine& rng = this->rng(); + + std::vector energy; + std::vector costheta; + + // Produce four samples from the original incident energy + for (int i : range(num_samples)) + { + Interaction result = interact(rng); + SCOPED_TRACE(result); + this->sanity_check(result); + + EXPECT_EQ(result.secondaries.data(), + this->secondary_allocator().get().data() + i); + + energy.push_back(result.secondaries.front().energy.value()); + costheta.push_back(dot_product(result.direction, + result.secondaries.front().direction)); + } + + EXPECT_EQ(num_samples, this->secondary_allocator().get().size()); + + // Note: these are "gold" values based on the host RNG. + static double const expected_energy[] = { + 0.0071536254658967, + 0.0044460343975567, + 0.0051966857337682, + 0.0010332114718837, + }; + static double const expected_costheta[] = { + 0.20412950131105, + 0.16112014959461, + 0.17413342534288, + 0.07778872399166, + }; + + EXPECT_VEC_SOFT_EQ(expected_energy, energy); + EXPECT_VEC_SOFT_EQ(expected_costheta, costheta); + + // Next sample should fail because we're out of secondary buffer space + { + Interaction result = interact(rng); + EXPECT_EQ(0, result.secondaries.size()); + EXPECT_EQ(Action::failed, result.action); + } + + // No interaction when max secondary energy is below production cut + { + CutoffParams::Input cut_inp; + cut_inp.materials = this->material_params(); + cut_inp.particles = this->particle_params(); + cut_inp.cutoffs = {{pdg::electron(), {{MevEnergy{0.01}, 0.1234}}}}; + this->set_cutoff_params(cut_inp); + + this->set_inc_particle(pdg::mu_minus(), MevEnergy{0.2}); + MuHadIonizationInteractor interact( + model_->host_ref(), + this->particle_track(), + this->cutoff_params()->get(MaterialId{0}), + this->direction(), + this->secondary_allocator()); + + Interaction result = interact(rng); + EXPECT_EQ(0, result.secondaries.size()); + EXPECT_EQ(Action::unchanged, result.action); + } +} + +TEST_F(BetheBlochTest, stress_test) +{ + unsigned int const num_samples = 10000; + std::vector avg_engine_samples; + std::vector avg_energy; + std::vector avg_costheta; + + for (real_type inc_e : {0.2, 1.0, 10.0, 1e2, 1e3, 1e4, 1e6, 1e8}) + { + SCOPED_TRACE("Incident energy: " + std::to_string(inc_e)); + this->set_inc_particle(pdg::mu_minus(), MevEnergy{inc_e}); + + RandomEngine& rng = this->rng(); + RandomEngine::size_type num_particles_sampled = 0; + real_type energy = 0; + real_type costheta = 0; + + // Loop over several incident directions + for (Real3 const& inc_dir : + {Real3{0, 0, 1}, Real3{1, 0, 0}, Real3{1e-9, 0, 1}, Real3{1, 1, 1}}) + { + SCOPED_TRACE("Incident direction: " + to_string(inc_dir)); + this->set_inc_direction(inc_dir); + this->resize_secondaries(num_samples); + + // Create interactor + MuHadIonizationInteractor interact( + model_->host_ref(), + this->particle_track(), + this->cutoff_params()->get(MaterialId{0}), + this->direction(), + this->secondary_allocator()); + + // Loop over many particles + for (unsigned int i = 0; i < num_samples; ++i) + { + Interaction result = interact(rng); + this->sanity_check(result); + + energy += result.secondaries.front().energy.value(); + costheta += dot_product(result.direction, + result.secondaries.front().direction); + } + EXPECT_EQ(num_samples, this->secondary_allocator().get().size()); + num_particles_sampled += num_samples; + } + avg_engine_samples.push_back(real_type(rng.count()) + / num_particles_sampled); + avg_energy.push_back(energy / num_particles_sampled); + avg_costheta.push_back(costheta / num_particles_sampled); + } + + // Gold values for average number of calls to RNG + static double const expected_avg_engine_samples[] + = {6.0069, 6.011, 6.0185, 6.0071, 6.0002, 6, 6, 6}; + static double const expected_avg_energy[] = { + 0.001820244315187, + 0.0030955371350616, + 0.0051011191515049, + 0.0071137840944271, + 0.010051462815079, + 0.0091324190520607, + 0.010625686784878, + 0.012563698826237, + }; + static double const expected_avg_costheta[] = { + 0.67374005035636, + 0.37023194384465, + 0.14030216439644, + 0.06933001323056, + 0.060759849187212, + 0.060385882216971, + 0.060522416476511, + 0.060172252455851, + }; + + EXPECT_VEC_SOFT_EQ(expected_avg_engine_samples, avg_engine_samples); + EXPECT_VEC_SOFT_EQ(expected_avg_energy, avg_energy); + EXPECT_VEC_SOFT_EQ(expected_avg_costheta, avg_costheta); +} + +//---------------------------------------------------------------------------// +} // namespace test +} // namespace celeritas diff --git a/test/celeritas/em/BraggICRU73QO.test.cc b/test/celeritas/em/BraggICRU73QO.test.cc index 252f0612ce..36bf62a1a8 100644 --- a/test/celeritas/em/BraggICRU73QO.test.cc +++ b/test/celeritas/em/BraggICRU73QO.test.cc @@ -8,8 +8,12 @@ #include "corecel/cont/Range.hh" #include "corecel/math/ArrayUtils.hh" #include "celeritas/Quantities.hh" -#include "celeritas/em/interactor/BraggICRU73QOInteractor.hh" +#include "celeritas/em/distribution/BraggICRU73QOEnergyDistribution.hh" +#include "celeritas/em/interactor/MuHadIonizationInteractor.hh" #include "celeritas/em/interactor/detail/PhysicsConstants.hh" +#include "celeritas/em/model/BraggModel.hh" +#include "celeritas/em/model/ICRU73QOModel.hh" +#include "celeritas/em/process/MuIonizationProcess.hh" #include "celeritas/phys/CutoffView.hh" #include "celeritas/phys/InteractionIO.hh" #include "celeritas/phys/InteractorHostTestBase.hh" @@ -52,32 +56,30 @@ class BraggICRU73QOTest : public InteractorHostTestBase cut_inp.cutoffs = {{pdg::electron(), {{MevEnergy{0.001}, 0.1234}}}}; this->set_cutoff_params(cut_inp); - // Set model data: default to ICRU73QO auto const& particles = *this->particle_params(); - data_.electron = particles.find(pdg::electron()); - data_.electron_mass = particles.get(data_.electron).mass(); - this->set_icru73qo(); + + // Set ICRU73QO model data + Applicability mu_minus; + mu_minus.particle = particles.find(pdg::mu_minus()); + mu_minus.lower = zero_quantity(); + mu_minus.upper + = MuIonizationProcess::Options{}.bragg_icru73qo_upper_limit; + icru73qo_model_ = std::make_shared( + ActionId{0}, particles, Model::SetApplicability{mu_minus}); + + // Set Bragg model data + Applicability mu_plus = mu_minus; + mu_plus.particle = particles.find(pdg::mu_plus()); + bragg_model_ = std::make_shared( + ActionId{0}, particles, Model::SetApplicability{mu_plus}); // Set default particle to muon with energy of 100 keV + inc_particle_ = pdg::mu_minus(); this->set_inc_particle(inc_particle_, MevEnergy{0.1}); this->set_inc_direction({0, 0, 1}); this->set_material("Cu"); } - void set_bragg() - { - inc_particle_ = pdg::mu_plus(); - data_.inc_particle = this->particle_params()->find(inc_particle_); - data_.lowest_kin_energy = detail::bragg_lowest_kin_energy(); - } - - void set_icru73qo() - { - inc_particle_ = pdg::mu_minus(); - data_.inc_particle = this->particle_params()->find(inc_particle_); - data_.lowest_kin_energy = detail::icru73qo_lowest_kin_energy(); - } - void sanity_check(Interaction const& interaction) const { // Check change to parent track @@ -92,7 +94,7 @@ class BraggICRU73QOTest : public InteractorHostTestBase auto const& electron = interaction.secondaries.front(); EXPECT_TRUE(electron); - EXPECT_EQ(data_.electron, electron.particle_id); + EXPECT_EQ(bragg_model_->host_ref().electron, electron.particle_id); EXPECT_GT(this->particle_track().energy().value(), electron.energy.value()); EXPECT_LT(0, electron.energy.value()); @@ -104,7 +106,8 @@ class BraggICRU73QOTest : public InteractorHostTestBase } protected: - BraggICRU73QOData data_; + std::shared_ptr bragg_model_; + std::shared_ptr icru73qo_model_; PDGNumber inc_particle_; }; @@ -112,12 +115,74 @@ class BraggICRU73QOTest : public InteractorHostTestBase // TESTS //---------------------------------------------------------------------------// +TEST_F(BraggICRU73QOTest, distribution) +{ + int num_samples = 100000; + int num_bins = 12; + + MevEnergy cutoff{1e-6}; + + std::vector counters; + std::vector min_energy; + std::vector max_energy; + for (real_type energy : {1e-4, 1e-3, 1e-2, 0.1, 0.2, 0.5, 1.0}) + { + this->set_inc_particle(pdg::mu_minus(), MevEnergy(energy)); + std::mt19937 rng; + + BraggICRU73QOEnergyDistribution sample( + this->particle_track(), + cutoff, + bragg_model_->host_ref().electron_mass); + real_type min = value_as(sample.min_secondary_energy()); + real_type max = value_as(sample.max_secondary_energy()); + + std::vector count(num_bins); + for ([[maybe_unused]] int i : range(num_samples)) + { + auto r = value_as(sample(rng)); + ASSERT_GE(r, min); + ASSERT_LE(r, max); + int bin = int((1 / r - 1 / min) / (1 / max - 1 / min) * num_bins); + CELER_ASSERT(bin >= 0 && bin < num_bins); + ++count[bin]; + } + counters.insert(counters.end(), count.begin(), count.end()); + min_energy.push_back(min); + max_energy.push_back(max); + } + + static int const expected_counters[] = { + 8267, 8480, 8377, 8202, 8301, 8415, 8239, 8429, 8370, 8266, 8338, 8316, + 8267, 8480, 8377, 8202, 8301, 8415, 8239, 8429, 8370, 8266, 8338, 8316, + 8267, 8480, 8377, 8202, 8301, 8415, 8239, 8429, 8370, 8266, 8338, 8316, + 8267, 8480, 8377, 8202, 8301, 8415, 8239, 8429, 8371, 8266, 8338, 8315, + 8267, 8480, 8377, 8202, 8301, 8415, 8239, 8429, 8371, 8266, 8338, 8315, + 8267, 8480, 8377, 8202, 8301, 8415, 8239, 8429, 8371, 8266, 8338, 8315, + 8267, 8480, 8377, 8202, 8301, 8415, 8239, 8429, 8371, 8266, 8338, 8315, + }; + static double const expected_min_energy[] + = {1e-06, 1e-06, 1e-06, 1e-06, 1e-06, 1e-06, 1e-06}; + static double const expected_max_energy[] = { + 1.9159563630249e-06, + 1.915964366753e-05, + 0.00019160444039615, + 0.001916844768863, + 0.0038354680957569, + 0.0096020089408745, + 0.019248476995285, + }; + EXPECT_VEC_EQ(expected_counters, counters); + EXPECT_VEC_SOFT_EQ(expected_min_energy, min_energy); + EXPECT_VEC_SOFT_EQ(expected_max_energy, max_energy); +} + TEST_F(BraggICRU73QOTest, basic) { std::vector energy; std::vector costheta; - auto sample = [&]() { + auto sample = [&](MuHadIonizationData const& data) { // Reserve 4 secondaries, one for each sample int const num_samples = 4; this->resize_secondaries(num_samples); @@ -128,8 +193,8 @@ TEST_F(BraggICRU73QOTest, basic) this->set_inc_particle(inc_particle_, MevEnergy{0.1}); // Create the interactor - BraggICRU73QOInteractor interact( - data_, + MuHadIonizationInteractor interact( + data, this->particle_track(), this->cutoff_params()->get(MaterialId{0}), this->direction(), @@ -163,8 +228,8 @@ TEST_F(BraggICRU73QOTest, basic) // No interaction when max secondary energy is below production cut { this->set_inc_particle(inc_particle_, MevEnergy{0.0011}); - BraggICRU73QOInteractor interact( - data_, + MuHadIonizationInteractor interact( + data, this->particle_track(), this->cutoff_params()->get(MaterialId{0}), this->direction(), @@ -178,8 +243,8 @@ TEST_F(BraggICRU73QOTest, basic) // Sample ICRU73QO model with incident mu- { - this->set_icru73qo(); - sample(); + inc_particle_ = pdg::mu_minus(); + sample(icru73qo_model_->host_ref()); static double const expected_energy[] = {0.0014458653777536, 0.001251648293082, @@ -195,8 +260,8 @@ TEST_F(BraggICRU73QOTest, basic) } // Sample Bragg model with incident mu+ { - this->set_bragg(); - sample(); + inc_particle_ = pdg::mu_plus(); + sample(bragg_model_->host_ref()); static double const expected_energy[] = {0.00022900204776481, 0.0014511488605566, @@ -217,7 +282,7 @@ TEST_F(BraggICRU73QOTest, stress_test) std::vector avg_energy; std::vector avg_costheta; - auto sample = [&]() { + auto sample = [&](MuHadIonizationData const& data) { unsigned int const num_samples = 10000; avg_engine_samples.clear(); @@ -245,12 +310,12 @@ TEST_F(BraggICRU73QOTest, stress_test) this->resize_secondaries(num_samples); // Create interactor - BraggICRU73QOInteractor interact( - data_, - this->particle_track(), - this->cutoff_params()->get(MaterialId{0}), - this->direction(), - this->secondary_allocator()); + MuHadIonizationInteractor + interact(data, + this->particle_track(), + this->cutoff_params()->get(MaterialId{0}), + this->direction(), + this->secondary_allocator()); // Loop over many particles for (unsigned int i = 0; i < num_samples; ++i) @@ -275,8 +340,8 @@ TEST_F(BraggICRU73QOTest, stress_test) // Sample ICRU73QO model with incident mu- { - this->set_icru73qo(); - sample(); + inc_particle_ = pdg::mu_minus(); + sample(icru73qo_model_->host_ref()); static double const expected_avg_engine_samples[] = {6.0027, 6.0021, 6.003, 6.0034, 6.0047}; @@ -296,8 +361,8 @@ TEST_F(BraggICRU73QOTest, stress_test) } // Sample Bragg model with incident mu+ { - this->set_bragg(); - sample(); + inc_particle_ = pdg::mu_plus(); + sample(bragg_model_->host_ref()); static double const expected_avg_engine_samples[] = {6.0004, 6.0004, 6.0006, 6.0003, 6.0005}; diff --git a/test/celeritas/em/MuBetheBloch.test.cc b/test/celeritas/em/MuBetheBloch.test.cc index f8d9cecbed..f75d090252 100644 --- a/test/celeritas/em/MuBetheBloch.test.cc +++ b/test/celeritas/em/MuBetheBloch.test.cc @@ -8,7 +8,10 @@ #include "corecel/cont/Range.hh" #include "corecel/math/ArrayUtils.hh" #include "celeritas/Quantities.hh" -#include "celeritas/em/interactor/MuBetheBlochInteractor.hh" +#include "celeritas/em/distribution/MuBBEnergyDistribution.hh" +#include "celeritas/em/interactor/MuHadIonizationInteractor.hh" +#include "celeritas/em/model/MuBetheBlochModel.hh" +#include "celeritas/em/process/MuIonizationProcess.hh" #include "celeritas/phys/CutoffView.hh" #include "celeritas/phys/InteractionIO.hh" #include "celeritas/phys/InteractorHostTestBase.hh" @@ -53,10 +56,15 @@ class MuBetheBlochTest : public InteractorHostTestBase // Set model data auto const& particles = *this->particle_params(); - data_.electron = particles.find(pdg::electron()); - data_.mu_minus = particles.find(pdg::mu_minus()); - data_.mu_plus = particles.find(pdg::mu_plus()); - data_.electron_mass = particles.get(data_.electron).mass(); + Applicability mu_minus; + mu_minus.particle = particles.find(pdg::mu_minus()); + mu_minus.lower + = MuIonizationProcess::Options{}.bragg_icru73qo_upper_limit; + mu_minus.upper = detail::high_energy_limit(); + Applicability mu_plus = mu_minus; + mu_plus.particle = particles.find(pdg::mu_plus()); + model_ = std::make_shared( + ActionId{0}, particles, Model::SetApplicability{mu_minus, mu_plus}); // Set default particle to muon with energy of 1 GeV this->set_inc_particle(pdg::mu_minus(), MevEnergy{1e3}); @@ -78,7 +86,7 @@ class MuBetheBlochTest : public InteractorHostTestBase auto const& electron = interaction.secondaries.front(); EXPECT_TRUE(electron); - EXPECT_EQ(data_.electron, electron.particle_id); + EXPECT_EQ(model_->host_ref().electron, electron.particle_id); EXPECT_GT(this->particle_track().energy().value(), electron.energy.value()); EXPECT_LT(0, electron.energy.value()); @@ -90,13 +98,75 @@ class MuBetheBlochTest : public InteractorHostTestBase } protected: - MuBetheBlochData data_; + std::shared_ptr model_; }; //---------------------------------------------------------------------------// // TESTS //---------------------------------------------------------------------------// +TEST_F(MuBetheBlochTest, distribution) +{ + int num_samples = 100000; + int num_bins = 12; + + MevEnergy cutoff{0.001}; + + std::vector counters; + std::vector min_energy; + std::vector max_energy; + for (real_type energy : {0.2, 1.0, 10.0, 1e2, 1e3, 1e4, 1e5, 1e7}) + { + this->set_inc_particle(pdg::mu_minus(), MevEnergy(energy)); + std::mt19937 rng; + + MuBBEnergyDistribution sample( + this->particle_track(), cutoff, model_->host_ref().electron_mass); + real_type min = value_as(sample.min_secondary_energy()); + real_type max = value_as(sample.max_secondary_energy()); + + std::vector count(num_bins); + for ([[maybe_unused]] int i : range(num_samples)) + { + auto r = value_as(sample(rng)); + ASSERT_GE(r, min); + ASSERT_LE(r, max); + int bin = int((1 / r - 1 / min) / (1 / max - 1 / min) * num_bins); + CELER_ASSERT(bin >= 0 && bin < num_bins); + ++count[bin]; + } + counters.insert(counters.end(), count.begin(), count.end()); + min_energy.push_back(min); + max_energy.push_back(max); + } + + static int const expected_counters[] = { + 8273, 8487, 8379, 8203, 8303, 8419, 8249, 8422, 8366, 8265, 8334, 8300, + 8281, 8499, 8383, 8211, 8309, 8428, 8256, 8435, 8371, 8263, 8320, 8244, + 8294, 8514, 8391, 8225, 8326, 8440, 8269, 8446, 8391, 8279, 8321, 8104, + 8283, 8499, 8389, 8211, 8312, 8427, 8257, 8440, 8380, 8278, 8340, 8184, + 8272, 8488, 8374, 8191, 8308, 8410, 8232, 8419, 8370, 8281, 8338, 8317, + 8278, 8490, 8370, 8191, 8311, 8404, 8246, 8413, 8371, 8255, 8335, 8336, + 8308, 8467, 8361, 8164, 8309, 8404, 8226, 8442, 8359, 8296, 8330, 8334, + 8353, 8428, 8388, 8264, 8315, 8312, 8227, 8499, 8323, 8280, 8266, 8345, + }; + static double const expected_min_energy[] + = {0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001}; + static double const expected_max_energy[] = { + 0.0038354680957569, + 0.019248476995285, + 0.20048052363148, + 2.7972680400033, + 100.69707462436, + 4855.7535710157, + 90256.629501068, + 9989193.9209199, + }; + EXPECT_VEC_EQ(expected_counters, counters); + EXPECT_VEC_SOFT_EQ(expected_min_energy, min_energy); + EXPECT_VEC_SOFT_EQ(expected_max_energy, max_energy); +} + TEST_F(MuBetheBlochTest, basic) { // Reserve 4 secondaries, one for each sample @@ -104,11 +174,12 @@ TEST_F(MuBetheBlochTest, basic) this->resize_secondaries(num_samples); // Create the interactor - MuBetheBlochInteractor interact(data_, - this->particle_track(), - this->cutoff_params()->get(MaterialId{0}), - this->direction(), - this->secondary_allocator()); + MuHadIonizationInteractor interact( + model_->host_ref(), + this->particle_track(), + this->cutoff_params()->get(MaterialId{0}), + this->direction(), + this->secondary_allocator()); RandomEngine& rng = this->rng(); std::vector energy; @@ -160,8 +231,8 @@ TEST_F(MuBetheBlochTest, basic) this->set_cutoff_params(cut_inp); this->set_inc_particle(pdg::mu_minus(), MevEnergy{0.2}); - MuBetheBlochInteractor interact( - data_, + MuHadIonizationInteractor interact( + model_->host_ref(), this->particle_track(), this->cutoff_params()->get(MaterialId{0}), this->direction(), @@ -199,8 +270,8 @@ TEST_F(MuBetheBlochTest, stress_test) this->resize_secondaries(num_samples); // Create interactor - MuBetheBlochInteractor interact( - data_, + MuHadIonizationInteractor interact( + model_->host_ref(), this->particle_track(), this->cutoff_params()->get(MaterialId{0}), this->direction(), diff --git a/test/celeritas/phys/ProcessBuilder.test.cc b/test/celeritas/phys/ProcessBuilder.test.cc index 107f644ab2..9a1418e94f 100644 --- a/test/celeritas/phys/ProcessBuilder.test.cc +++ b/test/celeritas/phys/ProcessBuilder.test.cc @@ -582,7 +582,9 @@ TEST_F(ProcessBuilderTest, mu_ionization) // Test model auto models = process->build_models(ActionIdIter{}); - ASSERT_EQ(3, models.size()); + // Note that newer versions of Geant4 use the \c G4MuBetheBloch model for + // all energies above 200 so will only have three models + ASSERT_EQ(4, models.size()); ASSERT_TRUE(models[0]); EXPECT_EQ("ioni-icru73qo", models[0]->label()); EXPECT_EQ(1, models[0]->applicability().size()); @@ -590,8 +592,11 @@ TEST_F(ProcessBuilderTest, mu_ionization) EXPECT_EQ("ioni-bragg", models[1]->label()); EXPECT_EQ(1, models[1]->applicability().size()); ASSERT_TRUE(models[2]); - EXPECT_EQ("ioni-mu-bethe-bloch", models[2]->label()); - auto all_applic = models[2]->applicability(); + EXPECT_EQ("ioni-bethe-bloch", models[2]->label()); + EXPECT_EQ(2, models[2]->applicability().size()); + ASSERT_TRUE(models[3]); + EXPECT_EQ("ioni-mu-bethe-bloch", models[3]->label()); + auto all_applic = models[3]->applicability(); ASSERT_EQ(2, all_applic.size()); for (auto mat_id : range(MaterialId{this->material()->num_materials()}))