decay/interactor/MuDecayInteractor.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/decay/interactor/MuDecayInteractor.hh
0006 //---------------------------------------------------------------------------//
0007 #pragma once
0008
0009 #include "corecel/Constants.hh"
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/random/distribution/RejectionSampler.hh"
0015 #include "geocel/random/IsotropicDistribution.hh"
0016 #include "celeritas/Quantities.hh"
0017 #include "celeritas/decay/data/MuDecayData.hh"
0018 #include "celeritas/phys/FourVector.hh"
0019 #include "celeritas/phys/Interaction.hh"
0020 #include "celeritas/phys/ParticleTrackView.hh"
0021 #include "celeritas/phys/Secondary.hh"
0022
0023 namespace celeritas
0024 {
0025 //---------------------------------------------------------------------------//
0026 /*!
0027  * Perform muon decay.
0028  *
0029  * Only one decay channel is implemented, with muons decaying to
0030  * \f[
0031  * \mu^- \longrightarrow e^- \bar{\nu}_e \nu_\mu
0032  * \f]
0033  * or
0034  * \f[
0035  * \mu^+ \longrightarrow e^+ \nu_e \bar{\nu}_\mu
0036  * \f].
0037  *
0038  * This interactor follows \c G4MuonDecayChannel::DecayIt and the Physics
0039  * Reference Manual, Release 11.2, section 4.2.3. The sampling happens at the
0040  * muon's rest frame, with the result being uniformly rotated and finally
0041  * boosted to the lab frame.
0042  *
0043  * As it is a three-body decay, the energy sampling happens for \f$ e^\pm \f$
0044  * and \f$ \nu_e(\bar{\nu}_e) \f$, with the \f$ \nu_\mu(\bar{\nu}_\mu) \f$
0045  * final energy directly calculated from energy conservation. The sampling loop
0046  * selects fractional energies \f$[0,1)\f$ for the first two particles. A
0047  * fraction of 1 yields the maximum possible kinetic energy for said particle,
0048  * defined as \f$ E_\text{max} = \frac{m_\mu}{2} - m_e \f$. For the electron
0049  * neutrino, its energy fraction \f$ f_{E_{\nu_e}} \f$ is sampled from the PDF
0050  * \f$ f(x) = 6x(1-x), \ x \in [0,1) \f$, using the rejection method, with
0051  * proposal distribution PDF \f$ g(x) = U(0,1) \f$ and bounding constant
0052  * \f$ M = 1.5 \f$. The charged lepton's fractional energy \f$ f_{E_e} \f$ is
0053  * then selected uniformly obeying \f$ f_{E_{\nu_e}} + f_{E_e} >= 1 \f$ . The
0054  * remaining fractional energy for the muon neutrino is
0055  * \f$ f_{E_{\nu_\mu}} = 2 - f_{E_e} - f_{E_{\nu_e}} \f$.
0056  *
0057  * \note Neutrinos are currently not returned by this interactor as they are
0058  * not tracked down or transported and would significantly increase secondary
0059  * memory allocation usage. See discussion and commit history at
0060  * https://github.com/celeritas-project/celeritas/pull/1456
0061  */
0062 class MuDecayInteractor
0063 {
0064   public:
0065     // Construct with shared and state data
0066     inline CELER_FUNCTION
0067     MuDecayInteractor(MuDecayData const& shared,
0068                       ParticleTrackView const& particle,
0069                       Real3 const& inc_direction,
0070                       StackAllocator<Secondary>& allocate);
0071
0072     // Sample an interaction with the given RNG
0073     template<class Engine>
0074     inline CELER_FUNCTION Interaction operator()(Engine& rng);
0075
0076   private:
0077     //// TYPES ////
0078
0079     using Energy = units::MevEnergy;
0080     using Momentum = units::MevMomentum;
0081     using Mass = units::MevMass;
0082
0083     //// DATA ////
0084
0085     // Constant data
0086     MuDecayData const& shared_;
0087     // Incident muon energy
0088     Energy const inc_energy_;
0089     // Allocate space for secondary particles (electron only)
0090     StackAllocator<Secondary>& allocate_;
0091     // Define decay channel based on muon or anti-muon primary
0092     ParticleId sec_id_;
0093     // Incident muon four vector
0094     FourVector inc_fourvec_;
0095     // Maximum electron energy [MeV]
0096     real_type max_energy_;
0097
0098     //// HELPER FUNCTIONS ////
0099
0100     // Boost four vector from the rest frame to the lab frame
0101     inline CELER_FUNCTION FourVector to_lab_frame(Real3 const& dir,
0102                                                   Momentum momentum,
0103                                                   Mass mass) const;
0104
0105     // Calculate particle momentum (or kinetic energy) in the center of mass
0106     inline CELER_FUNCTION Momentum calc_momentum(real_type energy_frac,
0107                                                  Mass mass) const;
0108 };
0109
0110 //---------------------------------------------------------------------------//
0111 // INLINE DEFINITIONS
0112 //---------------------------------------------------------------------------//
0113 /*!
0114  * Construct with shared and state data.
0115  *
0116  * \note Geant4 physics manual defines \f$ E_{max} = m_\mu / 2 \f$,
0117  * while the source code (since v10.2.0 at least) defines
0118  * \f$ E_{max} = m_\mu / 2 - m_e \f$ .
0119  * The source code implementation leads to a total CM energy of ~104.6
0120  * MeV instead of the expected 105.7 MeV (muon mass), which is achieved by
0121  * using the physics manual definition.
0122  */
0123 CELER_FUNCTION
0124 MuDecayInteractor::MuDecayInteractor(MuDecayData const& shared,
0125                                      ParticleTrackView const& particle,
0126                                      Real3 const& inc_direction,
0127                                      StackAllocator<Secondary>& allocate)
0128     : shared_(shared)
0129     , inc_energy_(particle.energy())
0130     , allocate_(allocate)
0131     , sec_id_((particle.particle_id() == shared_.mu_minus_id)
0132                   ? shared_.electron_id
0133                   : shared_.positron_id)
0134     , inc_fourvec_{FourVector::from_particle(particle, inc_direction)}
0135     , max_energy_(real_type{0.5} * shared_.muon_mass.value()
0136                   - shared_.electron_mass.value())
0137 {
0138     CELER_EXPECT(shared_);
0139     CELER_EXPECT(particle.particle_id() == shared_.mu_minus_id
0140                  || particle.particle_id() == shared_.mu_plus_id);
0141 }
0142
0143 //---------------------------------------------------------------------------//
0144 /*!
0145  * Sample the muon decay.
0146  */
0147 template<class Engine>
0148 CELER_FUNCTION Interaction MuDecayInteractor::operator()(Engine& rng)
0149 {
0150     // Allocate secondaries
0151     Secondary* secondaries = allocate_(1);
0152     if (secondaries == nullptr)
0153     {
0154         // Failed to allocate secondaries
0155         return Interaction::from_failure();
0156     }
0157
0158     real_type electron_energy_frac{};
0159     real_type electron_nu_energy_frac{};
0160     do
0161     {
0162         do
0163         {
0164             electron_nu_energy_frac = generate_canonical(rng);
0165         } while (RejectionSampler(
0166             electron_nu_energy_frac * (real_type{1} - electron_nu_energy_frac),
0167             real_type{0.25})(rng));
0168
0169         electron_energy_frac = generate_canonical(rng);
0170     } while (electron_nu_energy_frac + electron_energy_frac < real_type{1});
0171
0172     // Decay isotropically in rest frame and boost secondaries to the lab frame
0173     auto charged_lep_fv = this->to_lab_frame(
0174         IsotropicDistribution{}(rng),
0175         this->calc_momentum(electron_energy_frac, shared_.electron_mass),
0176         shared_.electron_mass);
0177
0178     // Return charged lepton only
0179     Interaction result = Interaction::from_absorption();
0180     result.secondaries = {secondaries, 1};
0181     result.secondaries[0].particle_id = sec_id_;
0182     // Interaction stores kinetic energy; FourVector stores total energy
0183     result.secondaries[0].energy
0184         = Energy{charged_lep_fv.energy - shared_.electron_mass.value()};
0185     result.secondaries[0].direction = make_unit_vector(charged_lep_fv.mom);
0186
0187     return result;
0188 }
0189
0190 //---------------------------------------------------------------------------//
0191 /*!
0192  * Boost secondary to the lab frame.
0193  *
0194  * \note This assumes the primary to be at rest and, thus, there is no need
0195  * to perform an inverse boost of the primary at the CM frame.
0196  */
0197 CELER_FUNCTION FourVector MuDecayInteractor::to_lab_frame(Real3 const& dir,
0198                                                           Momentum momentum,
0199                                                           Mass mass) const
0200 {
0201     CELER_EXPECT(is_soft_unit_vector(dir));
0202     CELER_EXPECT(momentum > zero_quantity());
0203     CELER_EXPECT(mass >= zero_quantity());
0204
0205     auto lepton_fv = FourVector::from_mass_momentum(mass, momentum, dir);
0206     boost(boost_vector(inc_fourvec_), &lepton_fv);
0207
0208     return lepton_fv;
0209 }
0210
0211 //---------------------------------------------------------------------------//
0212 /*!
0213  * Calculate final particle momentum (or kinetic energy) from its sampled
0214  * fractional energy.
0215  */
0216 CELER_FUNCTION auto
0217 MuDecayInteractor::calc_momentum(real_type energy_frac,
0218                                  Mass mass) const -> Momentum
0219 {
0220     return Momentum{std::sqrt(ipow<2>(energy_frac * max_energy_)
0221                               + 2 * energy_frac * max_energy_ * mass.value())};
0222 }
0223
0224 //---------------------------------------------------------------------------//
0225 }  // namespace celeritas