em/interactor/MuHadIonizationInteractor.hh

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