em/interactor/MuPairProductionInteractor.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/MuPairProductionInteractor.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/ArrayOperators.hh"
0014 #include "corecel/math/ArrayUtils.hh"
0015 #include "corecel/random/distribution/UniformRealDistribution.hh"
0016 #include "celeritas/Constants.hh"
0017 #include "celeritas/Quantities.hh"
0018 #include "celeritas/em/data/MuPairProductionData.hh"
0019 #include "celeritas/em/distribution/MuAngularDistribution.hh"
0020 #include "celeritas/em/distribution/MuPPEnergyDistribution.hh"
0021 #include "celeritas/mat/ElementView.hh"
0022 #include "celeritas/phys/CutoffView.hh"
0023 #include "celeritas/phys/Interaction.hh"
0024 #include "celeritas/phys/ParticleTrackView.hh"
0025 #include "celeritas/phys/Secondary.hh"
0026
0027 namespace celeritas
0028 {
0029 //---------------------------------------------------------------------------//
0030 /*!
0031  * Perform electron-positron pair production by muons.
0032  *
0033  * \note This performs the same sampling routine as in Geant4's
0034  * G4MuPairProductionModel and as documented in the Geant4 Physics Reference
0035  * Manual (Release 11.1) section 11.3.
0036  */
0037 class MuPairProductionInteractor
0038 {
0039   public:
0040     // Construct with shared and state data
0041     inline CELER_FUNCTION
0042     MuPairProductionInteractor(NativeCRef<MuPairProductionData> const& shared,
0043                                ParticleTrackView const& particle,
0044                                CutoffView const& cutoffs,
0045                                ElementView const& element,
0046                                Real3 const& inc_direction,
0047                                StackAllocator<Secondary>& allocate);
0048
0049     // Sample an interaction with the given RNG
0050     template<class Engine>
0051     inline CELER_FUNCTION Interaction operator()(Engine& rng);
0052
0053   private:
0054     //// TYPES ////
0055
0056     using Energy = units::MevEnergy;
0057     using Mass = units::MevMass;
0058     using Momentum = units::MevMomentum;
0059     using UniformRealDist = UniformRealDistribution<real_type>;
0060
0061     //// DATA ////
0062
0063     // Shared model data
0064     NativeCRef<MuPairProductionData> const& shared_;
0065     // Allocate space for the secondary particle
0066     StackAllocator<Secondary>& allocate_;
0067     // Incident direction
0068     Real3 const& inc_direction_;
0069     // Incident particle energy [MeV]
0070     Energy inc_energy_;
0071     // Incident particle mass
0072     Mass inc_mass_;
0073     // Incident particle momentum [MeV / c]
0074     real_type inc_momentum_;
0075     // Sample the azimuthal angle
0076     UniformRealDist sample_phi_;
0077     // Sampler for the electron-positron pair energy
0078     MuPPEnergyDistribution sample_energy_;
0079
0080     //// HELPER FUNCTIONS ////
0081
0082     // Calculate the secondary particle momentum from the sampled energy
0083     inline CELER_FUNCTION real_type calc_momentum(Energy) const;
0084 };
0085
0086 //---------------------------------------------------------------------------//
0087 // INLINE DEFINITIONS
0088 //---------------------------------------------------------------------------//
0089 /*!
0090  * Construct with shared and state data.
0091  */
0092 CELER_FUNCTION MuPairProductionInteractor::MuPairProductionInteractor(
0093     NativeCRef<MuPairProductionData> const& shared,
0094     ParticleTrackView const& particle,
0095     CutoffView const& cutoffs,
0096     ElementView const& element,
0097     Real3 const& inc_direction,
0098     StackAllocator<Secondary>& allocate)
0099     : shared_(shared)
0100     , allocate_(allocate)
0101     , inc_direction_(inc_direction)
0102     , inc_energy_(particle.energy())
0103     , inc_mass_(particle.mass())
0104     , inc_momentum_(value_as<Momentum>(particle.momentum()))
0105     , sample_phi_(0, real_type(2 * constants::pi))
0106     , sample_energy_(shared, particle, cutoffs, element)
0107 {
0108     CELER_EXPECT(particle.particle_id() == shared.ids.mu_minus
0109                  || particle.particle_id() == shared.ids.mu_plus);
0110 }
0111
0112 //---------------------------------------------------------------------------//
0113 /*!
0114  * Simulate electron-positron pair production by muons.
0115  */
0116 template<class Engine>
0117 CELER_FUNCTION Interaction MuPairProductionInteractor::operator()(Engine& rng)
0118 {
0119     // Allocate secondary electron and positron
0120     Secondary* secondaries = allocate_(2);
0121     if (secondaries == nullptr)
0122     {
0123         // Failed to allocate space for a secondary
0124         return Interaction::from_failure();
0125     }
0126
0127     // Sample the electron and positron energies
0128     auto energy = sample_energy_(rng);
0129     Energy pair_energy = energy.electron + energy.positron;
0130
0131     // Sample the secondary directions
0132     MuAngularDistribution sample_costheta(inc_energy_, inc_mass_, pair_energy);
0133     real_type phi = sample_phi_(rng);
0134
0135     // Create the secondary electron
0136     Secondary& electron = secondaries[0];
0137     electron.particle_id = shared_.ids.electron;
0138     electron.energy = energy.electron;
0139     electron.direction
0140         = rotate(from_spherical(sample_costheta(rng), phi), inc_direction_);
0141
0142     // Create the secondary positron
0143     Secondary& positron = secondaries[1];
0144     positron.particle_id = shared_.ids.positron;
0145     positron.energy = energy.positron;
0146     positron.direction
0147         = rotate(from_spherical(sample_costheta(rng), phi + constants::pi),
0148                  inc_direction_);
0149
0150     // Construct interaction for change to the incident muon
0151     Interaction result;
0152     result.secondaries = {secondaries, 2};
0153     result.energy = inc_energy_ - pair_energy;
0154     result.direction = make_unit_vector(
0155         inc_momentum_ * inc_direction_
0156         - this->calc_momentum(electron.energy) * electron.direction
0157         - this->calc_momentum(positron.energy) * positron.direction);
0158
0159     return result;
0160 }
0161
0162 //---------------------------------------------------------------------------//
0163 /*!
0164  * Calculate the secondary particle momentum from the sampled energy.
0165  */
0166 CELER_FUNCTION real_type
0167 MuPairProductionInteractor::calc_momentum(Energy energy) const
0168 {
0169     return std::sqrt(ipow<2>(value_as<Energy>(energy))
0170                      + 2 * value_as<Mass>(shared_.electron_mass)
0171                            * value_as<Energy>(energy));
0172 }
0173
0174 //---------------------------------------------------------------------------//
0175 }  // namespace celeritas