mucf/interactor/DTMucfInteractor.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/DTMucfInteractor.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$ (dt)_\mu \f$ molecules.
0022  *
0023  * Fusion channels:
0024  * - \f$ \alpha + \mu + n \f$
0025  * - \f$ (\alpha)_\mu + n \f$
0026  *
0027  * \warning This implementation has an incorrect energy and momentum
0028  * conservation implementation. Acceleron assumes an isotropic direction for
0029  * both neutron and muon in the \f$ \alpha + \mu + n \f$ channel, which leads
0030  * to the alpha particle either conserving energy or momentum but not both
0031  * simultaneously. The current implementation results in a roughly correct
0032  * total energy within \f$ K_\text{total} = [17.5, 17.9] \f$ MeV, instead of
0033  * the expected 17.6 MeV.
0034  */
0035 class DTMucfInteractor
0036 {
0037   public:
0038     enum class Channel
0039     {
0040         alpha_muon_neutron,  //!< \f$ \alpha + \mu + n \f$
0041         muonicalpha_neutron,  //!< \f$ (\alpha)_\mu + n \f$
0042         size_
0043     };
0044
0045     // Construct from shared and state data
0046     inline CELER_FUNCTION DTMucfInteractor(
0047         NativeCRef<DTMixMucfData> const& data,
0048         Channel channel,
0049         StackAllocator<Secondary>& allocate);
0050
0051     // Sample an interaction with the given RNG
0052     template<class Engine>
0053     inline CELER_FUNCTION Interaction operator()(Engine& rng);
0054
0055   private:
0056     // Shared constant physics properties
0057     NativeCRef<DTMixMucfData> const& data_;
0058     // Selected fusion channel
0059     Channel channel_{Channel::size_};
0060     // Allocate space for secondary particles
0061     StackAllocator<Secondary>& allocate_;
0062     // Number of secondaries per channel
0063     EnumArray<Channel, size_type> num_secondaries_{
0064         3,  // alpha_muon_neutron
0065         2  // muonicalpha_neutron
0066     };
0067
0068     // Outgoing neutron kinetic energy
0069     inline CELER_FUNCTION units::MevEnergy neutron_kinetic_energy() const
0070     {
0071         return units::MevEnergy{14.1};
0072     }
0073
0074     // Total fusion kinetic energy
0075     inline CELER_FUNCTION units::MevEnergy total_kinetic_energy() const
0076     {
0077         return units::MevEnergy{17.6};
0078     }
0079 };
0080
0081 //---------------------------------------------------------------------------//
0082 // INLINE DEFINITIONS
0083 //---------------------------------------------------------------------------//
0084 /*!
0085  * Construct with shared data and channel selection.
0086  */
0087 CELER_FUNCTION DTMucfInteractor::DTMucfInteractor(
0088     NativeCRef<DTMixMucfData> const& data,
0089     Channel const channel,
0090     StackAllocator<Secondary>& allocate)
0091     : data_(data), channel_(channel), allocate_(allocate)
0092 {
0093     CELER_EXPECT(data_);
0094     CELER_EXPECT(channel_ < Channel::size_);
0095 }
0096
0097 //---------------------------------------------------------------------------//
0098 /*!
0099  * Sample a dt muonic molecule fusion.
0100  */
0101 template<class Engine>
0102 CELER_FUNCTION Interaction DTMucfInteractor::operator()(Engine& rng)
0103 {
0104     // Allocate space for the final fusion channel
0105     Secondary* sec = allocate_(num_secondaries_[channel_]);
0106     if (sec == nullptr)
0107     {
0108         // Failed to allocate space for secondaries
0109         return Interaction::from_failure();
0110     }
0111
0112     size_type const neutron_idx{0};  // Both channels
0113     size_type const muon_idx{1}, alpha_idx{2};  // Channel::alpha_muon_neutron
0114     size_type const muonicalpha_idx{1};  // Channel::muonicalpha_neutron
0115
0116     IsotropicDistribution sample_isotropic;
0117
0118     // Neutron is the same on both cases: 14.1 MeV with random direction
0119     sec[neutron_idx] = detail::sample_mucf_secondary(
0120         data_.particle_ids.neutron, this->neutron_kinetic_energy(), rng);
0121
0122     switch (channel_)
0123     {
0124         case Channel::alpha_muon_neutron: {
0125             // Muon: random direction with energy sampled from its CDF
0126             sec[muon_idx] = detail::sample_mucf_muon(
0127                 data_.particle_ids.mu_minus,
0128                 NonuniformGridCalculator{data_.muon_energy_cdf, data_.reals},
0129                 rng);
0130
0131             // Alpha: Final state calculated via momentum conservation
0132             sec[alpha_idx]
0133                 = detail::calc_third_secondary(sec[neutron_idx],
0134                                                data_.particle_masses.neutron,
0135                                                sec[muon_idx],
0136                                                data_.particle_masses.mu_minus,
0137                                                data_.particle_ids.alpha,
0138                                                data_.particle_masses.alpha);
0139             break;
0140         }
0141
0142         case Channel::muonicalpha_neutron: {
0143             // Muonic alpha: Equal and opposite momentum to neutron
0144             sec[muonicalpha_idx].particle_id = data_.particle_ids.muonic_alpha;
0145             sec[muonicalpha_idx].energy = this->total_kinetic_energy()
0146                                           - this->neutron_kinetic_energy();
0147             sec[muonicalpha_idx].direction
0148                 = detail::opposite(sec[neutron_idx].direction);
0149             break;
0150         }
0151
0152         default:
0153             CELER_ASSERT_UNREACHABLE();
0154     }
0155
0156     // Kill muon primary and generate fusion secondaries
0157     Interaction result = Interaction::from_absorption();
0158     result.secondaries = {sec, num_secondaries_[channel_]};
0159     return result;
0160 }
0161
0162 //---------------------------------------------------------------------------//
0163 }  // namespace celeritas