diff --git a/src/celeritas/mucf/data/DTMixMucfData.hh b/src/celeritas/mucf/data/DTMixMucfData.hh index 56d2426450..3586060459 100644 --- a/src/celeritas/mucf/data/DTMixMucfData.hh +++ b/src/celeritas/mucf/data/DTMixMucfData.hh @@ -30,8 +30,9 @@ struct MucfParticleIds //!@{ //! Elementary particles and nuclei - ParticleId neutron; ParticleId proton; + ParticleId triton; + ParticleId neutron; ParticleId alpha; ParticleId he3; //!@} @@ -48,8 +49,9 @@ struct MucfParticleIds //! Check whether all particles are assigned CELER_FUNCTION explicit operator bool() const { - return mu_minus && neutron && proton && alpha && he3 && muonic_hydrogen - && muonic_alpha && muonic_triton && muonic_he3; + return mu_minus && proton && triton && neutron && alpha && he3 + && muonic_hydrogen && muonic_alpha && muonic_triton + && muonic_he3; } }; @@ -64,8 +66,9 @@ struct MucfParticleMasses //!@{ //! Elementary particles and nuclei - units::MevMass neutron; units::MevMass proton; + units::MevMass triton; + units::MevMass neutron; units::MevMass alpha; units::MevMass he3; //!@} @@ -82,9 +85,10 @@ struct MucfParticleMasses //! Check whether all data are assigned CELER_FUNCTION explicit operator bool() const { - return mu_minus > zero_quantity() && neutron > zero_quantity() - && proton > zero_quantity() && alpha > zero_quantity() - && he3 > zero_quantity() && muonic_hydrogen > zero_quantity() + return mu_minus > zero_quantity() && proton > zero_quantity() + && triton > zero_quantity() && neutron > zero_quantity() + && alpha > zero_quantity() && he3 > zero_quantity() + && muonic_hydrogen > zero_quantity() && muonic_alpha > zero_quantity() && muonic_triton > zero_quantity() && muonic_he3 > zero_quantity(); diff --git a/src/celeritas/mucf/interactor/DDMucfInteractor.hh b/src/celeritas/mucf/interactor/DDMucfInteractor.hh index b75fa71c02..c76e71b5b6 100644 --- a/src/celeritas/mucf/interactor/DDMucfInteractor.hh +++ b/src/celeritas/mucf/interactor/DDMucfInteractor.hh @@ -8,13 +8,12 @@ #include "corecel/Macros.hh" #include "corecel/data/StackAllocator.hh" -#include "celeritas/mat/ElementView.hh" -#include "celeritas/mat/MaterialView.hh" #include "celeritas/mucf/data/DTMixMucfData.hh" #include "celeritas/phys/Interaction.hh" -#include "celeritas/phys/ParticleTrackView.hh" #include "celeritas/phys/Secondary.hh" +#include "detail/MucfInteractorUtils.hh" + namespace celeritas { //---------------------------------------------------------------------------// @@ -24,17 +23,25 @@ namespace celeritas * Fusion channels: * - \f$ ^3\text{He} + \mu + n \f$ * - \f$ (^3\text{He})_\mu + n \f$ - * - \f$ ^3\text{H} + \mu + p \f$ + * - \f$ t + \mu + p \f$ + * + * The mass ratios between \f$ ^3\text{He} \f$ and the neutron, and between + * tritium and the proton, are both about 3:1. This leads to the neutron and + * proton kinetic energies being 3/4 of the total kinetic energy available in + * their respective channels. + * + * \note This interactor has a similar implementation as \c DTMucfInteractor , + * where energy is not fully conserved. See its documentation for details. */ class DDMucfInteractor { public: - //! \todo Implement muonichydrogen3_proton (\f$ (^3\text{H})_\mu + p \f$) + //! \todo Implement muonictriton_proton (\f$ (t)_\mu + p \f$) enum class Channel { helium3_muon_neutron, //!< \f$ ^3\text{He} + \mu + n \f$ muonichelium3_neutron, //!< \f$ (^3\text{He})_\mu + n \f$ - hydrogen3_muon_proton, //!< \f$ ^3\text{H} + \mu + p \f$ + tritium_muon_proton, //!< \f$ t + \mu + p \f$ size_ }; @@ -52,20 +59,39 @@ class DDMucfInteractor // Shared constant physics properties NativeCRef const& data_; // Selected fusion channel - Channel channel_{Channel::size_}; + Channel channel_; // Allocate space for secondary particles StackAllocator& allocate_; // Number of secondaries per channel EnumArray num_secondaries_{ 3, // helium3_muon_neutron 2, // muonichelium3_neutron - 3 // hydrogen3_muon_proton + 3 // tritium_muon_proton }; - // Sample Interaction secondaries - template - inline CELER_FUNCTION Span - sample_secondaries(Secondary* secondaries /*, other args */, Engine&); + // Neutron channels total kinetic energy + inline CELER_FUNCTION units::MevEnergy total_energy_neutron_channels() const + { + return units::MevEnergy{3.3}; + } + + // Proton channel total kinetic energy + inline CELER_FUNCTION units::MevEnergy total_energy_proton_channel() const + { + return units::MevEnergy{4.03}; + } + + // Outgoing neutron kinetic energy + inline CELER_FUNCTION units::MevEnergy neutron_kinetic_energy() const + { + return 0.75 * this->total_energy_neutron_channels(); + } + + // Outgoing proton kinetic energy + inline CELER_FUNCTION units::MevEnergy proton_kinetic_energy() const + { + return 0.75 * this->total_energy_proton_channel(); + } }; //---------------------------------------------------------------------------// @@ -76,7 +102,7 @@ class DDMucfInteractor */ CELER_FUNCTION DDMucfInteractor::DDMucfInteractor(NativeCRef const& data, - Channel const channel, + Channel channel, StackAllocator& allocate) : data_(data), channel_(channel), allocate_(allocate) { @@ -92,36 +118,88 @@ template CELER_FUNCTION Interaction DDMucfInteractor::operator()(Engine& rng) { // Allocate space for the final fusion channel - Secondary* secondaries = allocate_(num_secondaries_[channel_]); - if (secondaries == nullptr) + Secondary* sec = allocate_(num_secondaries_[channel_]); + if (sec == nullptr) { // Failed to allocate space for secondaries return Interaction::from_failure(); } - // Kill primary and generate secondaries - Interaction result = Interaction::from_absorption(); - result.secondaries = this->sample_secondaries(secondaries, rng); + // Short aliases for secondary indices + size_type const neutron_idx{0}, proton_idx{0}; + size_type const muon_idx{1}, muonic_he3{1}; + size_type const he3_idx{2}, tritium_idx{2}; - return result; -} + IsotropicDistribution sample_isotropic; -//---------------------------------------------------------------------------// -/*! - * Sample the secondaries of the selected channel. - * - * Since secondaries come from an at rest interaction, their final state is - * a simple combination of random direction + momentum conservation - */ -template -CELER_FUNCTION Span -DDMucfInteractor::sample_secondaries(Secondary* secondaries /*, other args */, - Engine&) -{ - // TODO: switch on channel_ - CELER_ASSERT_UNREACHABLE(); + switch (channel_) + { + case Channel::helium3_muon_neutron: { + // Neutron: random direction with known energy + sec[neutron_idx] = detail::sample_mucf_secondary( + data_.particle_ids.neutron, this->neutron_kinetic_energy(), rng); + + // Muon: random direction with energy sampled from its CDF + sec[muon_idx] = detail::sample_mucf_muon( + data_.particle_ids.mu_minus, + NonuniformGridCalculator{data_.muon_energy_cdf, data_.reals}, + rng); + + // Helium-3: momentum conservation + sec[he3_idx] + = detail::calc_third_secondary(sec[neutron_idx], + data_.particle_masses.neutron, + sec[muon_idx], + data_.particle_masses.mu_minus, + data_.particle_ids.he3, + data_.particle_masses.he3); + break; + } + + case Channel::muonichelium3_neutron: { + // Neutron: random direction with known energy + sec[neutron_idx] = detail::sample_mucf_secondary( + data_.particle_ids.neutron, this->neutron_kinetic_energy(), rng); - return Span{secondaries, num_secondaries_[channel_]}; + // Muonic helium-3: momentum conservation + sec[muonic_he3].particle_id = data_.particle_ids.muonic_he3; + sec[muonic_he3].energy = this->total_energy_neutron_channels() + - this->neutron_kinetic_energy(); + sec[muonic_he3].direction + = detail::opposite(sec[neutron_idx].direction); + break; + } + + case Channel::tritium_muon_proton: { + // Proton: random direction with known energy + sec[proton_idx] = detail::sample_mucf_secondary( + data_.particle_ids.proton, this->proton_kinetic_energy(), rng); + + // Muon: random direction with energy sampled from its CDF + sec[muon_idx] = detail::sample_mucf_muon( + data_.particle_ids.mu_minus, + NonuniformGridCalculator{data_.muon_energy_cdf, data_.reals}, + rng); + + // tritium: momentum conservation + sec[tritium_idx] + = detail::calc_third_secondary(sec[proton_idx], + data_.particle_masses.proton, + sec[muon_idx], + data_.particle_masses.mu_minus, + data_.particle_ids.triton, + data_.particle_masses.triton); + break; + } + + default: + CELER_ASSERT_UNREACHABLE(); + } + + // Kill primary and generate secondaries + Interaction result = Interaction::from_absorption(); + result.secondaries = {sec, num_secondaries_[channel_]}; + return result; } //---------------------------------------------------------------------------// diff --git a/src/celeritas/mucf/interactor/DTMucfInteractor.hh b/src/celeritas/mucf/interactor/DTMucfInteractor.hh index fb7b131a91..ac6c5dec9f 100644 --- a/src/celeritas/mucf/interactor/DTMucfInteractor.hh +++ b/src/celeritas/mucf/interactor/DTMucfInteractor.hh @@ -8,14 +8,12 @@ #include "corecel/Macros.hh" #include "corecel/data/StackAllocator.hh" -#include "corecel/math/ArrayUtils.hh" -#include "corecel/random/distribution/GenerateCanonical.hh" -#include "corecel/random/distribution/IsotropicDistribution.hh" -#include "celeritas/grid/NonuniformGridCalculator.hh" #include "celeritas/mucf/data/DTMixMucfData.hh" #include "celeritas/phys/Interaction.hh" #include "celeritas/phys/Secondary.hh" +#include "detail/MucfInteractorUtils.hh" + namespace celeritas { //---------------------------------------------------------------------------// @@ -74,14 +72,11 @@ class DTMucfInteractor return units::MevEnergy{14.1}; } - // Calculate momentum magnitude from energy and mass - inline CELER_FUNCTION real_type calc_momentum( - units::MevEnergy const energy, units::MevMass const mass) const; - - // Calculate kinetic energy from momentum and mass - inline CELER_FUNCTION units::MevEnergy - calc_kinetic_energy(Real3 const& momentum_vec, - units::MevMass const mass) const; + // Total fusion kinetic energy + inline CELER_FUNCTION units::MevEnergy total_kinetic_energy() const + { + return units::MevEnergy{17.6}; + } }; //---------------------------------------------------------------------------// @@ -120,57 +115,38 @@ CELER_FUNCTION Interaction DTMucfInteractor::operator()(Engine& rng) size_type const muonicalpha_idx{1}; // Channel::muonicalpha_neutron IsotropicDistribution sample_isotropic; - // Grid range is [0,1), with its domain being the muon energy in MeV - NonuniformGridCalculator sample_muon_energy(data_.muon_energy_cdf, - data_.reals); // Neutron is the same on both cases: 14.1 MeV with random direction - sec[neutron_idx].particle_id = data_.particle_ids.neutron; - sec[neutron_idx].energy = this->neutron_kinetic_energy(); - sec[neutron_idx].direction = sample_isotropic(rng); + sec[neutron_idx] = detail::sample_mucf_secondary( + data_.particle_ids.neutron, this->neutron_kinetic_energy(), rng); switch (channel_) { case Channel::alpha_muon_neutron: { // Muon: random direction with energy sampled from its CDF - sec[muon_idx].particle_id = data_.particle_ids.mu_minus; - sec[muon_idx].direction = sample_isotropic(rng); - sec[muon_idx].energy = units::MevEnergy{ - sample_muon_energy(generate_canonical(rng))}; + sec[muon_idx] = detail::sample_mucf_muon( + data_.particle_ids.mu_minus, + NonuniformGridCalculator{data_.muon_energy_cdf, data_.reals}, + rng); // Alpha: Final state calculated via momentum conservation - sec[alpha_idx].particle_id = data_.particle_ids.alpha; - - // Momentum: p_alpha = - (p_neutron + p_muon) - auto const neutron_momentum = this->calc_momentum( - sec[neutron_idx].energy, data_.particle_masses.neutron); - auto const muon_momentum = this->calc_momentum( - sec[muon_idx].energy, data_.particle_masses.mu_minus); - - Real3 alpha_momentum_vec; - for (auto i : range(3)) - { - alpha_momentum_vec[i] - = -(sec[neutron_idx].direction[i] * neutron_momentum - + sec[muon_idx].direction[i] * muon_momentum); - } - sec[alpha_idx].direction = make_unit_vector(alpha_momentum_vec); - - // Energy: K = sqrt(p^2 + m^2) - m - sec[alpha_idx].energy = this->calc_kinetic_energy( - alpha_momentum_vec, data_.particle_masses.alpha); + sec[alpha_idx] + = detail::calc_third_secondary(sec[neutron_idx], + data_.particle_masses.neutron, + sec[muon_idx], + data_.particle_masses.mu_minus, + data_.particle_ids.alpha, + data_.particle_masses.alpha); break; } case Channel::muonicalpha_neutron: { // Muonic alpha: Equal and opposite momentum to neutron sec[muonicalpha_idx].particle_id = data_.particle_ids.muonic_alpha; - sec[muonicalpha_idx].energy = sec[neutron_idx].energy; - for (auto i : range(3)) - { - sec[muonicalpha_idx].direction[i] - = -sec[neutron_idx].direction[i]; - } + sec[muonicalpha_idx].energy = this->total_kinetic_energy() + - this->neutron_kinetic_energy(); + sec[muonicalpha_idx].direction + = detail::opposite(sec[neutron_idx].direction); break; } @@ -184,34 +160,5 @@ CELER_FUNCTION Interaction DTMucfInteractor::operator()(Engine& rng) return result; } -//---------------------------------------------------------------------------// -/*! - * Calculate momentum magnitude from particle energy and mass via - * \f$ p = \sqrt{K^2 + 2mK} \f$ . - */ -CELER_FUNCTION real_type DTMucfInteractor::calc_momentum( - units::MevEnergy const energy, units::MevMass const mass) const - -{ - return std::sqrt(ipow<2>(value_as(energy)) - + 2 * value_as(mass) - * value_as(energy)); -} - -//---------------------------------------------------------------------------// -/*! - * Calculate kinetic energy given a particle's momentum and mass via - * \f$ K = \sqrt{p^2 - m^2} - m\f$ . - */ -CELER_FUNCTION units::MevEnergy -DTMucfInteractor::calc_kinetic_energy(Real3 const& momentum_vec, - units::MevMass const mass) const -{ - real_type momentum_mag = norm(momentum_vec); - return units::MevEnergy{std::sqrt(ipow<2>(momentum_mag) - + ipow<2>(value_as(mass))) - - value_as(mass)}; -} - //---------------------------------------------------------------------------// } // namespace celeritas diff --git a/src/celeritas/mucf/interactor/detail/MucfInteractorUtils.hh b/src/celeritas/mucf/interactor/detail/MucfInteractorUtils.hh new file mode 100644 index 0000000000..213b3b0aaf --- /dev/null +++ b/src/celeritas/mucf/interactor/detail/MucfInteractorUtils.hh @@ -0,0 +1,136 @@ +//------------------------------- -*- C++ -*- -------------------------------// +// Copyright Celeritas contributors: see top-level COPYRIGHT file for details +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/mucf/interactor/detail/MucfInteractorUtils.hh +//---------------------------------------------------------------------------// +#pragma once + +#include "corecel/Macros.hh" +#include "corecel/Types.hh" +#include "corecel/math/Algorithms.hh" +#include "corecel/math/ArrayUtils.hh" +#include "corecel/random/distribution/GenerateCanonical.hh" +#include "corecel/random/distribution/IsotropicDistribution.hh" +#include "celeritas/Quantities.hh" +#include "celeritas/Units.hh" +#include "celeritas/grid/NonuniformGridCalculator.hh" +#include "celeritas/phys/Secondary.hh" + +namespace celeritas +{ +namespace detail +{ +//---------------------------------------------------------------------------// +/*! + * Calculate momentum magnitude from particle energy and mass via + * \f$ p = \sqrt{K^2 + 2mK} \f$ . + */ +inline CELER_FUNCTION real_type calc_momentum(units::MevEnergy const energy, + units::MevMass const mass) +{ + return std::sqrt(ipow<2>(value_as(energy)) + + 2 * value_as(mass) + * value_as(energy)); +} + +//---------------------------------------------------------------------------// +/*! + * Calculate kinetic energy given a particle's momentum and mass via + * \f$ K = \sqrt{p^2 - m^2} - m\f$ . + */ +inline CELER_FUNCTION units::MevEnergy +calc_kinetic_energy(Real3 const& momentum_vec, units::MevMass const mass) +{ + real_type momentum_mag = norm(momentum_vec); + return units::MevEnergy{std::sqrt(ipow<2>(momentum_mag) + + ipow<2>(value_as(mass))) + - value_as(mass)}; +} + +//---------------------------------------------------------------------------// +/*! + * Sample muCF secondary with isotropic random direction and known outgoing + * energy. + */ +template +inline CELER_FUNCTION Secondary sample_mucf_secondary(ParticleId pid, + units::MevEnergy energy, + Engine& rng) +{ + Secondary result; + result.particle_id = pid; + result.energy = energy; + result.direction = IsotropicDistribution()(rng); + return result; +} + +//---------------------------------------------------------------------------// +/*! + * Sample a muon secondary with isotropic random direction and energy from the + * provided CDF. + * + * \note The muon grid range is [0,1) and its domain is the energy in MeV. + */ +template +inline CELER_FUNCTION Secondary sample_mucf_muon( + ParticleId pid, NonuniformGridCalculator sample_energy, Engine& rng) +{ + Secondary result; + result.particle_id = pid; + result.energy = units::MevEnergy{sample_energy(generate_canonical(rng))}; + result.direction = IsotropicDistribution()(rng); + return result; +} + +//---------------------------------------------------------------------------// +/*! + * Return opposite direction. + */ +inline CELER_FUNCTION Real3 opposite(Real3 const& vec) +{ + Real3 result; + for (auto i : range(3)) + { + result[i] = -vec[i]; + } + return result; +} + +//---------------------------------------------------------------------------// +/*! + * Return the momentum of the third secondary from an at-rest three-body + * sampling, once the other two are known. + * + * \todo This may be expanded to do a full three-body energy + momentum + * conservation. + */ +inline CELER_FUNCTION Secondary +calc_third_secondary(Secondary sec_a, + units::MevMass const mass_a, + Secondary sec_b, + units::MevMass const mass_b, + ParticleId pid_c, + units::MevMass const mass_c) +{ + auto const momentum_a_mag = calc_momentum(sec_a.energy, mass_a); + auto const momentum_b_mag = calc_momentum(sec_b.energy, mass_b); + + // Calculate direction from momentum conservation: p3 = - (p1 + p2) + Real3 momentum_vec; + for (auto i : range(3)) + { + momentum_vec[i] = -(sec_a.direction[i] * momentum_a_mag + + sec_b.direction[i] * momentum_b_mag); + } + + Secondary result; + result.particle_id = pid_c; + result.direction = make_unit_vector(momentum_vec); + result.energy = calc_kinetic_energy(momentum_vec, mass_c); + return result; +} + +//---------------------------------------------------------------------------// +} // namespace detail +} // namespace celeritas diff --git a/src/celeritas/phys/PDGNumber.hh b/src/celeritas/phys/PDGNumber.hh index f97f1fdf13..e1f1a2c3ec 100644 --- a/src/celeritas/phys/PDGNumber.hh +++ b/src/celeritas/phys/PDGNumber.hh @@ -137,6 +137,7 @@ CELER_DEFINE_PDGNUMBER(proton, 2212) CELER_DEFINE_PDGNUMBER(anti_proton, -2212) // Ions +CELER_DEFINE_PDGNUMBER(hydrogen, 1000010010) CELER_DEFINE_PDGNUMBER(deuteron, 1000010020) CELER_DEFINE_PDGNUMBER(anti_deuteron, -1000010020) CELER_DEFINE_PDGNUMBER(triton, 1000010030) @@ -147,9 +148,11 @@ CELER_DEFINE_PDGNUMBER(alpha, 1000020040) CELER_DEFINE_PDGNUMBER(anti_alpha, -1000020040) // Muonic nuclei +CELER_DEFINE_PDGNUMBER(muonic_hydrogen, 2000010010) CELER_DEFINE_PDGNUMBER(muonic_deuteron, 2000010020) CELER_DEFINE_PDGNUMBER(muonic_triton, 2000010030) CELER_DEFINE_PDGNUMBER(muonic_alpha, 2000020040) +CELER_DEFINE_PDGNUMBER(muonic_he3, 2000020030) //!@} diff --git a/test/celeritas/CMakeLists.txt b/test/celeritas/CMakeLists.txt index dccf4a402b..5a21748cbb 100644 --- a/test/celeritas/CMakeLists.txt +++ b/test/celeritas/CMakeLists.txt @@ -151,8 +151,9 @@ celeritas_add_test(decay/MuDecay.test.cc ${_needs_double}) #-----------------------------------------------------------------------------# # MuCF -celeritas_add_test(mucf/MuonicAtomSelector.test.cc) +celeritas_add_test(mucf/DDMucfInteractor.test.cc) celeritas_add_test(mucf/DTMucfInteractor.test.cc) +celeritas_add_test(mucf/MuonicAtomSelector.test.cc) #-----------------------------------------------------------------------------# # EM diff --git a/test/celeritas/mucf/DDMucfInteractor.test.cc b/test/celeritas/mucf/DDMucfInteractor.test.cc new file mode 100644 index 0000000000..283534290e --- /dev/null +++ b/test/celeritas/mucf/DDMucfInteractor.test.cc @@ -0,0 +1,309 @@ +//------------------------------- -*- C++ -*- -------------------------------// +// Copyright Celeritas contributors: see top-level COPYRIGHT file for details +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/mucf/DDMucfInteractor.test.cc +//---------------------------------------------------------------------------// +#include "celeritas/mucf/interactor/DDMucfInteractor.hh" + +#include "corecel/cont/Range.hh" +#include "corecel/math/ArrayUtils.hh" +#include "celeritas/Quantities.hh" +#include "celeritas/grid/NonuniformGridBuilder.hh" +#include "celeritas/inp/MucfPhysics.hh" + +#include "MucfInteractorHostTestBase.hh" +#include "celeritas_test.hh" + +namespace celeritas +{ +namespace test +{ +//---------------------------------------------------------------------------// +// TEST HARNESS +//---------------------------------------------------------------------------// + +class DDMucfInteractorTest : public MucfInteractorHostTestBase +{ + using Base = MucfInteractorHostTestBase; + using Channel = DDMucfInteractor::Channel; + using MevEnergy = units::MevEnergy; + using MevMass = units::MevMass; + + protected: + void SetUp() override + { + auto const& params = *this->particle_params(); + this->set_material("hdt_fuel"); + + HostVal host_data; + + // Set up particle IDs + host_data.particle_ids.mu_minus = params.find(pdg::mu_minus()); + host_data.particle_ids.neutron = params.find(pdg::neutron()); + host_data.particle_ids.proton = params.find(pdg::proton()); + host_data.particle_ids.he3 = params.find(pdg::he3()); + host_data.particle_ids.triton = params.find(pdg::triton()); + host_data.particle_ids.muonic_he3 = params.find(pdg::muonic_he3()); + + // Set up particle masses + host_data.particle_masses.mu_minus + = params.get(host_data.particle_ids.mu_minus).mass(); + host_data.particle_masses.neutron + = params.get(host_data.particle_ids.neutron).mass(); + host_data.particle_masses.proton + = params.get(host_data.particle_ids.proton).mass(); + host_data.particle_masses.he3 + = params.get(host_data.particle_ids.he3).mass(); + host_data.particle_masses.triton + = params.get(host_data.particle_ids.triton).mass(); + host_data.particle_masses.muonic_he3 + = params.get(host_data.particle_ids.muonic_he3).mass(); + + // Set up muon energy CDF + auto const inp_data = inp::MucfPhysics::from_default(); + NonuniformGridBuilder build_grid_record{&host_data.reals}; + host_data.muon_energy_cdf = build_grid_record(inp_data.muon_energy_cdf); + + // Construct collection + data_ = ParamsDataStore{std::move(host_data)}; + + // At-rest muon primary + this->set_inc_particle(pdg::mu_minus(), MevEnergy{0.0}); + this->set_inc_direction({1, 0, 0}); + } + + // Detailed validation of the interaction result + void + validate_interaction(Interaction const& interaction, Channel channel) const + { + EXPECT_LT(channel, Channel::size_); + + // Primary muon should be killed + EXPECT_EQ(Action::absorbed, interaction.action); + + auto const& host_data = data_.host_ref(); + auto const& sec = interaction.secondaries; + + // Verify channel-specific data + if (channel == Channel::helium3_muon_neutron) + { + ASSERT_EQ(num_secondaries_[channel], sec.size()); + + // First particle is the outgoing neutron with 0.75 * 3.3 MeV + EXPECT_EQ(host_data.particle_ids.neutron, sec[0].particle_id); + EXPECT_SOFT_EQ(0.75 * 3.3, sec[0].energy.value()); + + // Check particles + EXPECT_EQ(host_data.particle_ids.mu_minus, sec[1].particle_id); + EXPECT_EQ(host_data.particle_ids.he3, sec[2].particle_id); + + // Check approximate energy conservation + real_type total_kinetic_energy = 0; + for (auto const& s : interaction.secondaries) + { + total_kinetic_energy += s.energy.value(); + } + // The total kinetic energy range is within [3.2, 3.5] MeV + EXPECT_SOFT_NEAR(3.3, total_kinetic_energy, 0.3); + + // Check momentum conservation + auto const neutron_p_mag = this->calc_momentum( + sec[0].energy, host_data.particle_masses.neutron); + auto const muon_p_mag = this->calc_momentum( + sec[1].energy, host_data.particle_masses.mu_minus); + + Real3 he3_momentum, total_momentum; + for (auto i : range(3)) + { + real_type neutron_momentum_i = sec[0].direction[i] + * neutron_p_mag; + real_type muon_momentum_i = sec[1].direction[i] * muon_p_mag; + he3_momentum[i] = -(neutron_momentum_i + muon_momentum_i); + total_momentum[i] = neutron_momentum_i + muon_momentum_i + + he3_momentum[i]; + } + + EXPECT_VEC_SOFT_EQ(sec[2].direction, + make_unit_vector(he3_momentum)); + EXPECT_VEC_SOFT_EQ(Real3{}, total_momentum); + } + + else if (channel == Channel::muonichelium3_neutron) + { + ASSERT_EQ(num_secondaries_[channel], + interaction.secondaries.size()); + + // Check particle types + EXPECT_EQ(host_data.particle_ids.neutron, + interaction.secondaries[0].particle_id); + EXPECT_EQ(host_data.particle_ids.muonic_he3, + interaction.secondaries[1].particle_id); + + // First particle is the outgoing neutron with 0.75 * 3.3 MeV + EXPECT_SOFT_EQ(0.75 * 3.3, + interaction.secondaries[0].energy.value()); + + // Check directions are opposite + EXPECT_SOFT_EQ(-1.0, + dot_product(interaction.secondaries[0].direction, + interaction.secondaries[1].direction)); + } + + else if (channel == Channel::tritium_muon_proton) + { + ASSERT_EQ(num_secondaries_[channel], sec.size()); + + // First particle is the outgoing proton with 0.75 * 4.03 MeV + EXPECT_EQ(host_data.particle_ids.proton, sec[0].particle_id); + EXPECT_SOFT_EQ(0.75 * 4.03, sec[0].energy.value()); + + // Check particles + EXPECT_EQ(host_data.particle_ids.mu_minus, sec[1].particle_id); + EXPECT_EQ(host_data.particle_ids.triton, sec[2].particle_id); + + // Check approximate energy conservation + // The total kinetic energy is only very roughly 4.03 MeV due to + // simplistic sampling. + real_type total_kinetic_energy = 0; + for (auto const& s : interaction.secondaries) + { + total_kinetic_energy += s.energy.value(); + } + EXPECT_SOFT_NEAR(4.03, total_kinetic_energy, 0.5); + + // Check momentum conservation + auto const proton_p_mag = this->calc_momentum( + sec[0].energy, host_data.particle_masses.proton); + auto const muon_p_mag = this->calc_momentum( + sec[1].energy, host_data.particle_masses.mu_minus); + + Real3 triton_momentum, total_momentum; + for (auto i : range(3)) + { + real_type proton_momentum_i = sec[0].direction[i] + * proton_p_mag; + real_type muon_momentum_i = sec[1].direction[i] * muon_p_mag; + triton_momentum[i] = -(proton_momentum_i + muon_momentum_i); + total_momentum[i] = proton_momentum_i + muon_momentum_i + + triton_momentum[i]; + } + + EXPECT_VEC_SOFT_EQ(sec[2].direction, + make_unit_vector(triton_momentum)); + EXPECT_VEC_SOFT_EQ(Real3{}, total_momentum); + } + } + + // Momentum magnitude (p = \sqrt{K^2 + 2mK}) + real_type calc_momentum(MevEnergy energy, MevMass mass) const + { + return std::sqrt(ipow<2>(value_as(energy)) + + 2 * value_as(mass) + * value_as(energy)); + } + + protected: + ParamsDataStore data_; + EnumArray num_secondaries_{ + 3, // helium3_muon_neutron + 2, // muonichelium3_neutron + 3 // tritium_muon_proton + }; +}; + +//---------------------------------------------------------------------------// +// TESTS +//---------------------------------------------------------------------------// + +TEST_F(DDMucfInteractorTest, helium3_muon_neutron) +{ + auto const channel = DDMucfInteractor::Channel::helium3_muon_neutron; + + size_type const num_samples = 4; + this->resize_secondaries(num_samples * num_secondaries_[channel]); + + // Run interactor + DDMucfInteractor interact( + data_.host_ref(), channel, this->secondary_allocator()); + + auto& rng = this->rng(); + for ([[maybe_unused]] auto i : range(num_samples)) + { + auto result = interact(rng); + this->validate_interaction(result, channel); + } +} + +//---------------------------------------------------------------------------// +TEST_F(DDMucfInteractorTest, muonichelium3_neutron) +{ + auto const channel = DDMucfInteractor::Channel::muonichelium3_neutron; + + size_type const num_samples = 4; + this->resize_secondaries(num_samples * num_secondaries_[channel]); + + // Run interactor + DDMucfInteractor interact( + data_.host_ref(), channel, this->secondary_allocator()); + + auto& rng = this->rng(); + for ([[maybe_unused]] auto i : range(num_samples)) + { + auto result = interact(rng); + this->validate_interaction(result, channel); + } +} + +//---------------------------------------------------------------------------// +TEST_F(DDMucfInteractorTest, tritium_muon_proton) +{ + auto const channel = DDMucfInteractor::Channel::tritium_muon_proton; + + size_type const num_samples = 4; + this->resize_secondaries(num_samples * num_secondaries_[channel]); + + // Run interactor + DDMucfInteractor interact( + data_.host_ref(), channel, this->secondary_allocator()); + + auto& rng = this->rng(); + for ([[maybe_unused]] auto i : range(num_samples)) + { + auto result = interact(rng); + this->validate_interaction(result, channel); + } +} + +//---------------------------------------------------------------------------// +TEST_F(DDMucfInteractorTest, stress_test) +{ + size_type const num_samples = 10000; + real_type total_avg_secondaries{0}; + + for (auto channel : {DDMucfInteractor::Channel::helium3_muon_neutron, + DDMucfInteractor::Channel::muonichelium3_neutron, + DDMucfInteractor::Channel::tritium_muon_proton}) + { + this->resize_secondaries(num_samples * num_secondaries_[channel]); + + DDMucfInteractor interact( + data_.host_ref(), channel, this->secondary_allocator()); + + auto& rng = this->rng(); + for ([[maybe_unused]] auto i : range(num_samples)) + { + auto result = interact(rng); + total_avg_secondaries += result.secondaries.size(); + } + } + total_avg_secondaries /= 3 * num_samples; // Average over all channels + + // (3 + 2 + 3) / 3 + static real_type const expected_total_avg_secondaries{8. / 3.}; + EXPECT_SOFT_EQ(expected_total_avg_secondaries, total_avg_secondaries); +} + +//---------------------------------------------------------------------------// +} // namespace test +} // namespace celeritas diff --git a/test/celeritas/mucf/DTMucfInteractor.test.cc b/test/celeritas/mucf/DTMucfInteractor.test.cc index fc63f1b118..1d2e47fd71 100644 --- a/test/celeritas/mucf/DTMucfInteractor.test.cc +++ b/test/celeritas/mucf/DTMucfInteractor.test.cc @@ -11,7 +11,6 @@ #include "celeritas/Quantities.hh" #include "celeritas/grid/NonuniformGridBuilder.hh" #include "celeritas/inp/MucfPhysics.hh" -#include "celeritas/phys/InteractionIO.hh" #include "MucfInteractorHostTestBase.hh" #include "celeritas_test.hh" @@ -145,10 +144,6 @@ class DTMucfInteractorTest : public MucfInteractorHostTestBase EXPECT_EQ(host_data.particle_ids.muonic_alpha, interaction.secondaries[1].particle_id); - // Check kinetic energies are equal - EXPECT_SOFT_EQ(interaction.secondaries[0].energy.value(), - interaction.secondaries[1].energy.value()); - // Check directions are opposite EXPECT_SOFT_EQ(-1.0, dot_product(interaction.secondaries[0].direction, @@ -220,7 +215,7 @@ TEST_F(DTMucfInteractorTest, muonicalpha_neutron) //---------------------------------------------------------------------------// TEST_F(DTMucfInteractorTest, stress_test) { - size_type const num_samples = 10000; + size_type const num_samples = 100; real_type total_avg_secondaries{0}; for (auto channel : {DTMucfInteractor::Channel::alpha_muon_neutron, diff --git a/test/celeritas/mucf/MucfInteractorHostTestBase.cc b/test/celeritas/mucf/MucfInteractorHostTestBase.cc index e39de6f2bb..1889f4b3be 100644 --- a/test/celeritas/mucf/MucfInteractorHostTestBase.cc +++ b/test/celeritas/mucf/MucfInteractorHostTestBase.cc @@ -62,7 +62,8 @@ MucfInteractorHostBase::MucfInteractorHostBase() muon_mass, ElementaryCharge{1}, native_value_from(muon_decay_constant)}, - // Nuclei and ions + + // Ions {"proton", pdg::proton(), protium_mass, @@ -89,11 +90,33 @@ MucfInteractorHostBase::MucfInteractorHostBase() ElementaryCharge{2}, stable_decay_constant}, {"he3", pdg::he3(), he3_mass, ElementaryCharge{2}, stable_decay_constant}, + + // Muonic atoms + {"muonic_hydrogen", + pdg::muonic_hydrogen(), + protium_mass + muon_mass, + ElementaryCharge{1}, + stable_decay_constant}, + {"muonic_deuteron", + pdg::muonic_deuteron(), + deuterium_mass + muon_mass, + ElementaryCharge{1}, + stable_decay_constant}, + {"muonic_triton", + pdg::muonic_triton(), + tritium_mass + muon_mass, + ElementaryCharge{1}, + stable_decay_constant}, {"muonic_alpha", pdg::muonic_alpha(), alpha_mass + muon_mass, ElementaryCharge{1}, native_value_from(muon_decay_constant)}, + {"muonic_he3", + pdg::muonic_he3(), + he3_mass + muon_mass, + ElementaryCharge{1}, + native_value_from(muon_decay_constant)}, }; this->set_particle_params(std::move(par_inp));