em/interactor/MuHadIonizationInteractor.hh

0001 //----------------------------------*-C++-*----------------------------------//
0002 // Copyright 2024 UT-Battelle, LLC, and other Celeritas developers.
0003 // See the top-level COPYRIGHT file for details.
0004 // SPDX-License-Identifier: (Apache-2.0 OR MIT)
0005 //---------------------------------------------------------------------------//
0006 //! \file celeritas/em/interactor/MuHadIonizationInteractor.hh
0007 //---------------------------------------------------------------------------//
0008 #pragma once
0009
0010 #include "corecel/Macros.hh"
0011 #include "corecel/Types.hh"
0012 #include "corecel/data/StackAllocator.hh"
0013 #include "corecel/math/Algorithms.hh"
0014 #include "corecel/math/ArrayUtils.hh"
0015 #include "celeritas/Constants.hh"
0016 #include "celeritas/Quantities.hh"
0017 #include "celeritas/em/data/MuHadIonizationData.hh"
0018 #include "celeritas/phys/CutoffView.hh"
0019 #include "celeritas/phys/Interaction.hh"
0020 #include "celeritas/phys/ParticleTrackView.hh"
0021 #include "celeritas/phys/Secondary.hh"
0022 #include "celeritas/random/distribution/InverseSquareDistribution.hh"
0023 #include "celeritas/random/distribution/RejectionSampler.hh"
0024
0025 #include "detail/IoniFinalStateHelper.hh"
0026 #include "detail/PhysicsConstants.hh"
0027
0028 namespace celeritas
0029 {
0030 //---------------------------------------------------------------------------//
0031 /*!
0032  * Perform the discrete part of the muon or hadron ionization process.
0033  *
0034  * This simulates the production of delta rays by incident muons or hadrons.
0035  * The same basic sampling routine is used by multiple models, but the energy
0036  * of the secondary is sampled from a distribution unique to the model.
0037  *
0038  * \note This performs the same sampling routine as in Geant4's \c
0039  * G4BetheBlochModel, \c G4MuBetheBlochModel, \c G4BraggModel, and \c
0040  * G4ICRU73QOModel, as documented in the Geant4 Physics Reference Manual
0041  * release 11.2 sections 11.1 and 12.1.5.
0042  */
0043 template<class EnergySampler>
0044 class MuHadIonizationInteractor
0045 {
0046   public:
0047     //!@{
0048     //! \name Type aliases
0049     using Mass = units::MevMass;
0050     using Energy = units::MevEnergy;
0051     using Momentum = units::MevMomentum;
0052     //!@}
0053
0054   public:
0055     //! Construct with shared and state data
0056     inline CELER_FUNCTION
0057     MuHadIonizationInteractor(MuHadIonizationData const& shared,
0058                               ParticleTrackView const& particle,
0059                               CutoffView const& cutoffs,
0060                               Real3 const& inc_direction,
0061                               StackAllocator<Secondary>& allocate);
0062
0063     // Sample an interaction with the given RNG
0064     template<class Engine>
0065     inline CELER_FUNCTION Interaction operator()(Engine& rng);
0066
0067   private:
0068     //// DATA ////
0069
0070     // Allocate space for the secondary particle
0071     StackAllocator<Secondary>& allocate_;
0072     // Incident direction
0073     Real3 const& inc_direction_;
0074     // Incident particle energy [MeV]
0075     Energy inc_energy_;
0076     // Incident particle momentum [MeV / c]
0077     Momentum inc_momentum_;
0078     // Muon mass [MeV / c^2]
0079     Mass inc_mass_;
0080     // Electron mass [MeV / c^2]
0081     Mass electron_mass_;
0082     // Secondary electron ID
0083     ParticleId electron_id_;
0084     // Maximum energy of the secondary electron [MeV]
0085     real_type max_secondary_energy_;
0086     // Secondary electron energy distribution
0087     EnergySampler sample_energy_;
0088 };
0089
0090 //---------------------------------------------------------------------------//
0091 // INLINE DEFINITIONS
0092 //---------------------------------------------------------------------------//
0093 /*!
0094  * Construct with shared and state data.
0095  */
0096 template<class ES>
0097 CELER_FUNCTION MuHadIonizationInteractor<ES>::MuHadIonizationInteractor(
0098     MuHadIonizationData const& shared,
0099     ParticleTrackView const& particle,
0100     CutoffView const& cutoffs,
0101     Real3 const& inc_direction,
0102     StackAllocator<Secondary>& allocate)
0103     : allocate_(allocate)
0104     , inc_direction_(inc_direction)
0105     , inc_energy_(particle.energy())
0106     , inc_momentum_(particle.momentum())
0107     , inc_mass_(particle.mass())
0108     , electron_mass_(shared.electron_mass)
0109     , electron_id_(shared.electron)
0110     , sample_energy_(particle, cutoffs.energy(electron_id_), electron_mass_)
0111 {
0112     CELER_EXPECT(inc_energy_ > sample_energy_.min_secondary_energy());
0113 }
0114
0115 //---------------------------------------------------------------------------//
0116 /*!
0117  * Simulate discrete muon ionization.
0118  */
0119 template<class ES>
0120 template<class Engine>
0121 CELER_FUNCTION Interaction MuHadIonizationInteractor<ES>::operator()(Engine& rng)
0122 {
0123     if (sample_energy_.min_secondary_energy()
0124         >= sample_energy_.max_secondary_energy())
0125     {
0126         // No interaction if the maximum secondary energy is below the limit
0127         return Interaction::from_unchanged();
0128     }
0129
0130     // Allocate secondary electron
0131     Secondary* secondary = allocate_(1);
0132     if (secondary == nullptr)
0133     {
0134         // Failed to allocate space for a secondary
0135         return Interaction::from_failure();
0136     }
0137
0138     // Update kinematics of the final state and return the interaction
0139     return detail::IoniFinalStateHelper(inc_energy_,
0140                                         inc_direction_,
0141                                         inc_momentum_,
0142                                         inc_mass_,
0143                                         sample_energy_(rng),
0144                                         electron_mass_,
0145                                         electron_id_,
0146                                         secondary)(rng);
0147 }
0148
0149 //---------------------------------------------------------------------------//
0150 }  // namespace celeritas