mucf/interactor/DDMucfInteractor.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/mucf/interactor/DDMucfInteractor.hh
0006 //---------------------------------------------------------------------------//
0007 #pragma once
0008
0009 #include "corecel/Macros.hh"
0010 #include "corecel/data/StackAllocator.hh"
0011 #include "celeritas/mucf/data/DTMixMucfData.hh"
0012 #include "celeritas/phys/Interaction.hh"
0013 #include "celeritas/phys/Secondary.hh"
0014
0015 #include "detail/MucfInteractorUtils.hh"
0016
0017 namespace celeritas
0018 {
0019 //---------------------------------------------------------------------------//
0020 /*!
0021  * Muon-catalyzed fusion of \f$ (dd)_\mu \f$ molecules.
0022  *
0023  * Fusion channels:
0024  * - \f$ ^3\text{He} + \mu + n \f$
0025  * - \f$ (^3\text{He})_\mu + n \f$
0026  * - \f$ t + \mu + p \f$
0027  *
0028  * The mass ratios between \f$ ^3\text{He} \f$ and the neutron, and between
0029  * tritium and the proton, are both about 3:1. This leads to the neutron and
0030  * proton kinetic energies being 3/4 of the total kinetic energy available in
0031  * their respective channels.
0032  *
0033  * \note This interactor has a similar implementation as \c DTMucfInteractor ,
0034  * where energy is not fully conserved. See its documentation for details.
0035  */
0036 class DDMucfInteractor
0037 {
0038   public:
0039     //! \todo Implement muonictriton_proton (\f$ (t)_\mu + p \f$)
0040     enum class Channel
0041     {
0042         helium3_muon_neutron,  //!< \f$ ^3\text{He} + \mu + n \f$
0043         muonichelium3_neutron,  //!< \f$ (^3\text{He})_\mu + n \f$
0044         tritium_muon_proton,  //!< \f$ t + \mu + p \f$
0045         size_
0046     };
0047
0048     // Construct from shared and state data
0049     inline CELER_FUNCTION DDMucfInteractor(
0050         NativeCRef<DTMixMucfData> const& data,
0051         Channel channel,
0052         StackAllocator<Secondary>& allocate);
0053
0054     // Sample an interaction with the given RNG
0055     template<class Engine>
0056     inline CELER_FUNCTION Interaction operator()(Engine& rng);
0057
0058   private:
0059     // Shared constant physics properties
0060     NativeCRef<DTMixMucfData> const& data_;
0061     // Selected fusion channel
0062     Channel channel_;
0063     // Allocate space for secondary particles
0064     StackAllocator<Secondary>& allocate_;
0065     // Number of secondaries per channel
0066     EnumArray<Channel, size_type> num_secondaries_{
0067         3,  // helium3_muon_neutron
0068         2,  // muonichelium3_neutron
0069         3  // tritium_muon_proton
0070     };
0071
0072     // Neutron channels total kinetic energy
0073     inline CELER_FUNCTION units::MevEnergy total_energy_neutron_channels() const
0074     {
0075         return units::MevEnergy{3.3};
0076     }
0077
0078     // Proton channel total kinetic energy
0079     inline CELER_FUNCTION units::MevEnergy total_energy_proton_channel() const
0080     {
0081         return units::MevEnergy{4.03};
0082     }
0083
0084     // Outgoing neutron kinetic energy
0085     inline CELER_FUNCTION units::MevEnergy neutron_kinetic_energy() const
0086     {
0087         return 0.75 * this->total_energy_neutron_channels();
0088     }
0089
0090     // Outgoing proton kinetic energy
0091     inline CELER_FUNCTION units::MevEnergy proton_kinetic_energy() const
0092     {
0093         return 0.75 * this->total_energy_proton_channel();
0094     }
0095 };
0096
0097 //---------------------------------------------------------------------------//
0098 // INLINE DEFINITIONS
0099 //---------------------------------------------------------------------------//
0100 /*!
0101  * Construct with shared and state data.
0102  */
0103 CELER_FUNCTION DDMucfInteractor::DDMucfInteractor(
0104     NativeCRef<DTMixMucfData> const& data,
0105     Channel channel,
0106     StackAllocator<Secondary>& allocate)
0107     : data_(data), channel_(channel), allocate_(allocate)
0108 {
0109     CELER_EXPECT(data_);
0110     CELER_EXPECT(channel < Channel::size_);
0111 }
0112
0113 //---------------------------------------------------------------------------//
0114 /*!
0115  * Sample a dd muonic molecule fusion.
0116  */
0117 template<class Engine>
0118 CELER_FUNCTION Interaction DDMucfInteractor::operator()(Engine& rng)
0119 {
0120     // Allocate space for the final fusion channel
0121     Secondary* sec = allocate_(num_secondaries_[channel_]);
0122     if (sec == nullptr)
0123     {
0124         // Failed to allocate space for secondaries
0125         return Interaction::from_failure();
0126     }
0127
0128     // Short aliases for secondary indices
0129     size_type const neutron_idx{0}, proton_idx{0};
0130     size_type const muon_idx{1}, muonic_he3{1};
0131     size_type const he3_idx{2}, tritium_idx{2};
0132
0133     IsotropicDistribution sample_isotropic;
0134
0135     switch (channel_)
0136     {
0137         case Channel::helium3_muon_neutron: {
0138             // Neutron: random direction with known energy
0139             sec[neutron_idx]
0140                 = detail::sample_mucf_secondary(data_.particle_ids.neutron,
0141                                                 this->neutron_kinetic_energy(),
0142                                                 rng);
0143
0144             // Muon: random direction with energy sampled from its CDF
0145             sec[muon_idx] = detail::sample_mucf_muon(
0146                 data_.particle_ids.mu_minus,
0147                 NonuniformGridCalculator{data_.muon_energy_cdf, data_.reals},
0148                 rng);
0149
0150             // Helium-3: momentum conservation
0151             sec[he3_idx]
0152                 = detail::calc_third_secondary(sec[neutron_idx],
0153                                                data_.particle_masses.neutron,
0154                                                sec[muon_idx],
0155                                                data_.particle_masses.mu_minus,
0156                                                data_.particle_ids.he3,
0157                                                data_.particle_masses.he3);
0158             break;
0159         }
0160
0161         case Channel::muonichelium3_neutron: {
0162             // Neutron: random direction with known energy
0163             sec[neutron_idx]
0164                 = detail::sample_mucf_secondary(data_.particle_ids.neutron,
0165                                                 this->neutron_kinetic_energy(),
0166                                                 rng);
0167
0168             // Muonic helium-3: momentum conservation
0169             sec[muonic_he3].particle_id = data_.particle_ids.muonic_he3;
0170             sec[muonic_he3].energy = this->total_energy_neutron_channels()
0171                                      - this->neutron_kinetic_energy();
0172             sec[muonic_he3].direction
0173                 = detail::opposite(sec[neutron_idx].direction);
0174             break;
0175         }
0176
0177         case Channel::tritium_muon_proton: {
0178             // Proton: random direction with known energy
0179             sec[proton_idx] = detail::sample_mucf_secondary(
0180                 data_.particle_ids.proton, this->proton_kinetic_energy(), rng);
0181
0182             // Muon: random direction with energy sampled from its CDF
0183             sec[muon_idx] = detail::sample_mucf_muon(
0184                 data_.particle_ids.mu_minus,
0185                 NonuniformGridCalculator{data_.muon_energy_cdf, data_.reals},
0186                 rng);
0187
0188             // tritium: momentum conservation
0189             sec[tritium_idx]
0190                 = detail::calc_third_secondary(sec[proton_idx],
0191                                                data_.particle_masses.proton,
0192                                                sec[muon_idx],
0193                                                data_.particle_masses.mu_minus,
0194                                                data_.particle_ids.triton,
0195                                                data_.particle_masses.triton);
0196             break;
0197         }
0198
0199         default:
0200             CELER_ASSERT_UNREACHABLE();
0201     }
0202
0203     // Kill primary and generate secondaries
0204     Interaction result = Interaction::from_absorption();
0205     result.secondaries = {sec, num_secondaries_[channel_]};
0206     return result;
0207 }
0208
0209 //---------------------------------------------------------------------------//
0210 }  // namespace celeritas