em/interactor/MuBremsstrahlungInteractor.hh

0001 //----------------------------------*-C++-*----------------------------------//
0002 // Copyright 2021-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/MuBremsstrahlungInteractor.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/ArrayUtils.hh"
0014 #include "celeritas/Constants.hh"
0015 #include "celeritas/Quantities.hh"
0016 #include "celeritas/em/data/MuBremsstrahlungData.hh"
0017 #include "celeritas/em/xs/MuBremsDiffXsCalculator.hh"
0018 #include "celeritas/mat/ElementView.hh"
0019 #include "celeritas/mat/MaterialView.hh"
0020 #include "celeritas/phys/CutoffView.hh"
0021 #include "celeritas/phys/Interaction.hh"
0022 #include "celeritas/phys/InteractionUtils.hh"
0023 #include "celeritas/phys/ParticleTrackView.hh"
0024 #include "celeritas/phys/Secondary.hh"
0025 #include "celeritas/random/distribution/ReciprocalDistribution.hh"
0026 #include "celeritas/random/distribution/RejectionSampler.hh"
0027
0028 #include "detail/BremFinalStateHelper.hh"
0029
0030 namespace celeritas
0031 {
0032 //---------------------------------------------------------------------------//
0033 /*!
0034  * Perform muon bremsstrahlung interaction.
0035  *
0036  * This is a model for the Bremsstrahlung process for muons. Given an incident
0037  * muon, the class computes the change to the incident muon direction and
0038  * energy, and it adds a single secondary gamma to the secondary stack.
0039  *
0040  * \note This performs the same sampling routine as in Geant4's
0041  * G4MuBremsstrahlungModel class, as documented in section 11.2
0042  * of the Geant4 Physics Reference (release 10.6).
0043  */
0044 class MuBremsstrahlungInteractor
0045 {
0046     //!@{
0047     //! \name Type aliases
0048     using Energy = units::MevEnergy;
0049     using Mass = units::MevMass;
0050     using Momentum = units::MevMomentum;
0051     //!@}
0052
0053   public:
0054     // Construct with shared and state data
0055     inline CELER_FUNCTION
0056     MuBremsstrahlungInteractor(MuBremsstrahlungData const& shared,
0057                                ParticleTrackView const& particle,
0058                                Real3 const& inc_direction,
0059                                CutoffView const& cutoffs,
0060                                StackAllocator<Secondary>& allocate,
0061                                MaterialView const& material,
0062                                ElementComponentId elcomp_id);
0063
0064     // Sample an interaction with the given RNG
0065     template<class Engine>
0066     inline CELER_FUNCTION Interaction operator()(Engine& rng);
0067
0068   private:
0069     //// DATA ////
0070
0071     // Shared constant physics properties
0072     MuBremsstrahlungData const& shared_;
0073     // Incident direction
0074     Real3 const& inc_direction_;
0075     // Allocate space for one or more secondary particles
0076     StackAllocator<Secondary>& allocate_;
0077     // Element properties
0078     ElementView const element_;
0079     // Incident particle
0080     ParticleTrackView const& particle_;
0081     // Differential cross section calculator
0082     MuBremsDiffXsCalculator calc_dcs_;
0083     // Distribution to sample energy
0084     ReciprocalDistribution<> sample_energy_;
0085     // Envelope distribution for rejection sampling of gamma energy
0086     real_type envelope_;
0087
0088     //// CONSTANTS ////
0089
0090     //// HELPER FUNCTIONS ////
0091
0092     template<class Engine>
0093     inline CELER_FUNCTION real_type sample_cos_theta(real_type gamma_energy,
0094                                                      Engine& rng) const;
0095 };
0096
0097 //---------------------------------------------------------------------------//
0098 // INLINE DEFINITIONS
0099 //---------------------------------------------------------------------------//
0100 /*!
0101  * Construct with shared and state data.
0102  */
0103 CELER_FUNCTION MuBremsstrahlungInteractor::MuBremsstrahlungInteractor(
0104     MuBremsstrahlungData const& shared,
0105     ParticleTrackView const& particle,
0106     Real3 const& inc_direction,
0107     CutoffView const& cutoffs,
0108     StackAllocator<Secondary>& allocate,
0109     MaterialView const& material,
0110     ElementComponentId elcomp_id)
0111     : shared_(shared)
0112     , inc_direction_(inc_direction)
0113     , allocate_(allocate)
0114     , element_(material.make_element_view(elcomp_id))
0115     , particle_(particle)
0116     , calc_dcs_(
0117           element_, particle.energy(), particle.mass(), shared.electron_mass)
0118     , sample_energy_{value_as<Energy>(cutoffs.energy(shared.gamma)),
0119                      value_as<Energy>(particle_.energy())}
0120 {
0121     CELER_EXPECT(particle.particle_id() == shared_.mu_minus
0122                  || particle.particle_id() == shared_.mu_plus);
0123     CELER_EXPECT(particle_.energy() > cutoffs.energy(shared.gamma));
0124
0125     // Calculate rejection envelope: *assume* the highest cross section
0126     // is at its lowerst value
0127     real_type gamma_cutoff = value_as<Energy>(cutoffs.energy(shared.gamma));
0128     envelope_ = gamma_cutoff * calc_dcs_(Energy{gamma_cutoff});
0129 }
0130
0131 //---------------------------------------------------------------------------//
0132 /*!
0133  * Sample using the muon bremsstrahlung model.
0134  */
0135 template<class Engine>
0136 CELER_FUNCTION Interaction MuBremsstrahlungInteractor::operator()(Engine& rng)
0137 {
0138     // Allocate space for gamma
0139     Secondary* secondary = allocate_(1);
0140     if (secondary == nullptr)
0141     {
0142         // Failed to allocate space for a secondary
0143         return Interaction::from_failure();
0144     }
0145
0146     // Sample the energy transfer
0147     real_type gamma_energy;
0148     do
0149     {
0150         gamma_energy = sample_energy_(rng);
0151     } while (RejectionSampler{gamma_energy * calc_dcs_(Energy{gamma_energy}),
0152                               envelope_}(rng));
0153
0154     // Update kinematics of the final state and return this interaction
0155     return detail::BremFinalStateHelper(
0156         particle_.energy(),
0157         inc_direction_,
0158         particle_.momentum(),
0159         shared_.gamma,
0160         Energy{gamma_energy},
0161         this->sample_cos_theta(gamma_energy, rng),
0162         secondary)(rng);
0163 }
0164
0165 //---------------------------------------------------------------------------//
0166 /*!
0167  * Sample cosine of the angle between incident and secondary particles.
0168  */
0169 template<class Engine>
0170 CELER_FUNCTION real_type MuBremsstrahlungInteractor::sample_cos_theta(
0171     real_type gamma_energy, Engine& rng) const
0172 {
0173     real_type gamma = particle_.lorentz_factor();
0174     real_type r_max_sq = ipow<2>(
0175         gamma * constants::pi * real_type(0.5)
0176         * min(real_type(1.0),
0177               gamma * value_as<Mass>(particle_.mass()) / gamma_energy - 1));
0178     real_type a = generate_canonical(rng) * r_max_sq / (1 + r_max_sq);
0179
0180     return std::cos(std::sqrt(a / (1 - a)) / gamma);
0181 }
0182
0183 //---------------------------------------------------------------------------//
0184 }  // namespace celeritas