em/interactor/MuBremsstrahlungInteractor.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/MuBremsstrahlungInteractor.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/ArrayUtils.hh"
0013 #include "corecel/random/distribution/ReciprocalDistribution.hh"
0014 #include "corecel/random/distribution/RejectionSampler.hh"
0015 #include "celeritas/Constants.hh"
0016 #include "celeritas/Quantities.hh"
0017 #include "celeritas/em/data/MuBremsstrahlungData.hh"
0018 #include "celeritas/em/distribution/MuAngularDistribution.hh"
0019 #include "celeritas/em/xs/MuBremsDiffXsCalculator.hh"
0020 #include "celeritas/mat/ElementView.hh"
0021 #include "celeritas/mat/MaterialView.hh"
0022 #include "celeritas/phys/CutoffView.hh"
0023 #include "celeritas/phys/Interaction.hh"
0024 #include "celeritas/phys/InteractionUtils.hh"
0025 #include "celeritas/phys/ParticleTrackView.hh"
0026 #include "celeritas/phys/Secondary.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
0089 //---------------------------------------------------------------------------//
0090 // INLINE DEFINITIONS
0091 //---------------------------------------------------------------------------//
0092 /*!
0093  * Construct with shared and state data.
0094  */
0095 CELER_FUNCTION MuBremsstrahlungInteractor::MuBremsstrahlungInteractor(
0096     MuBremsstrahlungData const& shared,
0097     ParticleTrackView const& particle,
0098     Real3 const& inc_direction,
0099     CutoffView const& cutoffs,
0100     StackAllocator<Secondary>& allocate,
0101     MaterialView const& material,
0102     ElementComponentId elcomp_id)
0103     : shared_(shared)
0104     , inc_direction_(inc_direction)
0105     , allocate_(allocate)
0106     , element_(material.element_record(elcomp_id))
0107     , particle_(particle)
0108     , calc_dcs_(
0109           element_, particle.energy(), particle.mass(), shared.electron_mass)
0110     , sample_energy_{value_as<Energy>(cutoffs.energy(shared.gamma)),
0111                      value_as<Energy>(particle_.energy())}
0112 {
0113     CELER_EXPECT(particle.particle_id() == shared_.mu_minus
0114                  || particle.particle_id() == shared_.mu_plus);
0115     CELER_EXPECT(particle_.energy() > cutoffs.energy(shared.gamma));
0116
0117     // Calculate rejection envelope: *assume* the highest cross section
0118     // is at its lowest value
0119     real_type gamma_cutoff = value_as<Energy>(cutoffs.energy(shared.gamma));
0120     envelope_ = gamma_cutoff * calc_dcs_(Energy{gamma_cutoff});
0121 }
0122
0123 //---------------------------------------------------------------------------//
0124 /*!
0125  * Sample using the muon bremsstrahlung model.
0126  */
0127 template<class Engine>
0128 CELER_FUNCTION Interaction MuBremsstrahlungInteractor::operator()(Engine& rng)
0129 {
0130     // Allocate space for gamma
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     // Sample the energy transfer
0139     real_type gamma_energy;
0140     do
0141     {
0142         gamma_energy = sample_energy_(rng);
0143     } while (RejectionSampler{gamma_energy * calc_dcs_(Energy{gamma_energy}),
0144                               envelope_}(rng));
0145
0146     MuAngularDistribution sample_costheta(
0147         particle_.energy(), particle_.mass(), Energy{gamma_energy});
0148
0149     // Update kinematics of the final state and return this interaction
0150     return detail::BremFinalStateHelper(particle_.energy(),
0151                                         inc_direction_,
0152                                         particle_.momentum(),
0153                                         shared_.gamma,
0154                                         Energy{gamma_energy},
0155                                         sample_costheta(rng),
0156                                         secondary)(rng);
0157 }
0158
0159 //---------------------------------------------------------------------------//
0160 }  // namespace celeritas