From 4f9b7880d4984ac8b49c8dce8eff2d632589dfdd Mon Sep 17 00:00:00 2001 From: Stefano Tognini Date: Tue, 20 Jan 2026 16:31:52 -0500 Subject: [PATCH 01/42] Draft implementation of the `DTMucfInteractor` --- src/celeritas/mucf/data/DTMixMucfData.hh | 50 +++++++- .../mucf/interactor/DTMucfInteractor.hh | 107 ++++++++++++++---- src/celeritas/mucf/model/DTMixMucfModel.cc | 26 +++-- src/celeritas/phys/PDGNumber.hh | 1 + 4 files changed, 150 insertions(+), 34 deletions(-) diff --git a/src/celeritas/mucf/data/DTMixMucfData.hh b/src/celeritas/mucf/data/DTMixMucfData.hh index 63affe3867..63109034fa 100644 --- a/src/celeritas/mucf/data/DTMixMucfData.hh +++ b/src/celeritas/mucf/data/DTMixMucfData.hh @@ -53,6 +53,44 @@ struct MucfParticleIds } }; +//---------------------------------------------------------------------------// +/*! + * Particle masses used by the \c DTMixMucfModel . + */ +struct MucfParticleMasses +{ + //! Primary + units::MevMass mu_minus; + + //!@{ + //! Elementary particles and nuclei + units::MevMass neutron; + units::MevMass proton; + units::MevMass alpha; + units::MevMass he3; + //!@} + + //!@{ + //! Muonic atoms + units::MevMass muonic_hydrogen; + units::MevMass muonic_deuteron; + units::MevMass muonic_triton; + units::MevMass muonic_alpha; + units::MevMass muonic_he3; + //!@} + + //! 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() + && muonic_alpha > zero_quantity() + && muonic_triton > zero_quantity() + && muonic_he3 > zero_quantity(); + } +}; + //---------------------------------------------------------------------------// /*! * Data for for the \c DTMixMucfModel . @@ -67,8 +105,9 @@ struct DTMixMucfData using GridRecord = NonuniformGridRecord; using CycleTimesArray = EnumArray>; - //! Particle IDs - MucfParticleIds particles; + //! Particles + MucfParticleIds particle_ids; + MucfParticleMasses particle_masses; //! Muon CDF energy grid for sampling outgoing muCF muons GridRecord muon_energy_cdf; //! \todo Verify energy unit @@ -87,8 +126,8 @@ struct DTMixMucfData //! Check whether the data are assigned explicit CELER_FUNCTION operator bool() const { - return particles && muon_energy_cdf && !mucfmatid_to_matid.empty() - && !cycle_times.empty() + return particle_ids && particle_masses && muon_energy_cdf + && !mucfmatid_to_matid.empty() && !cycle_times.empty() && (mucfmatid_to_matid.size() == cycle_times.size()); } @@ -99,7 +138,8 @@ struct DTMixMucfData CELER_EXPECT(other); //! \todo Finish implementation - this->particles = other.particles; + this->particle_ids = other.particle_ids; + this->particle_masses = other.particle_masses; this->reals = other.reals; this->muon_energy_cdf = other.muon_energy_cdf; this->mucfmatid_to_matid = other.mucfmatid_to_matid; diff --git a/src/celeritas/mucf/interactor/DTMucfInteractor.hh b/src/celeritas/mucf/interactor/DTMucfInteractor.hh index 568bf3ba0e..9a68b57103 100644 --- a/src/celeritas/mucf/interactor/DTMucfInteractor.hh +++ b/src/celeritas/mucf/interactor/DTMucfInteractor.hh @@ -8,11 +8,12 @@ #include "corecel/Macros.hh" #include "corecel/data/StackAllocator.hh" -#include "celeritas/mat/ElementView.hh" -#include "celeritas/mat/MaterialView.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/ParticleTrackView.hh" #include "celeritas/phys/Secondary.hh" namespace celeritas @@ -52,16 +53,28 @@ class DTMucfInteractor Channel channel_{Channel::size_}; // Allocate space for secondary particles StackAllocator& allocate_; + // Number of secondaries per channel EnumArray num_secondaries_{ 3, // alpha_muon_neutron 2 // muonicalpha_neutron }; - // Sample Interaction secondaries - template - inline CELER_FUNCTION Span - sample_secondaries(Secondary* secondaries /*, other args */, Engine&); + // Total energy generated by the dt fusion reaction + inline CELER_FUNCTION units::MevEnergy total_fusion_energy() const + { + return units::MevEnergy{17.6}; + } + + // Neutron kinetic energy + inline CELER_FUNCTION units::MevEnergy neutron_kinetic_energy() const + { + return units::MevEnergy{14.1}; + } + + // Calculate momentum magnitude from energy and mass + inline CELER_FUNCTION real_type calc_momentum(units::MevEnergy energy, + units::MevMass mass) const; }; //---------------------------------------------------------------------------// @@ -83,6 +96,9 @@ DTMucfInteractor::DTMucfInteractor(NativeCRef const& data, //---------------------------------------------------------------------------// /*! * Sample a dt muonic molecule fusion. + * + * \note Since secondaries come from an at rest interaction, their final state + * is a simple combination of random direction + momentum conservation. */ template CELER_FUNCTION Interaction DTMucfInteractor::operator()(Engine& rng) @@ -95,28 +111,79 @@ CELER_FUNCTION Interaction DTMucfInteractor::operator()(Engine& rng) return Interaction::from_failure(); } + IsotropicDistribution sample_isotropic; + NonuniformGridCalculator sample_muon_energy(data_.muon_energy_cdf, + data_.reals); + + // Neutron is the same on both cases: 14.1 MeV with random direction + secondaries[0].particle_id = data_.particle_ids.neutron; + secondaries[0].energy = this->neutron_kinetic_energy(); + secondaries[0].direction = sample_isotropic(rng); + + switch (channel_) + { + case Channel::alpha_muon_neutron: { + // Muon: Final state sampled from CDF + secondaries[1].particle_id = data_.particle_ids.mu_minus; + secondaries[1].direction = sample_isotropic(rng); + secondaries[1].energy = units::MevEnergy{ + sample_muon_energy(generate_canonical(rng))}; + + // Alpha: Final state calculated via energy/momentum conservation: + // E_alpha = E_total - E_neutron - E_muon + // p_alpha = - (p_neutron + p_muon) + auto const neutron_momentum_mag = this->calc_momentum( + secondaries[0].energy, data_.particle_masses.neutron); + auto const muon_momentum_mag = this->calc_momentum( + secondaries[1].energy, data_.particle_masses.mu_minus); + Real3 alpha_momentum_vec; + for (auto i : range(3)) + { + alpha_momentum_vec[i] + = -(secondaries[0].direction[i] * neutron_momentum_mag + + secondaries[1].direction[i] * muon_momentum_mag); + } + + secondaries[2].particle_id = data_.particle_ids.alpha; + secondaries[2].direction = make_unit_vector(alpha_momentum_vec); + secondaries[2].energy = this->total_fusion_energy() + - secondaries[0].energy + - secondaries[1].energy; + break; + } + + case Channel::muonicalpha_neutron: { + // Muonic alpha: Equal and opposite momentum to neutron + secondaries[1].particle_id = data_.particle_ids.muonic_alpha; + secondaries[1].energy = secondaries[0].energy; + for (auto i : range(3)) + { + secondaries[1].direction[i] = -secondaries[0].direction[i]; + } + break; + } + + default: + CELER_ASSERT_UNREACHABLE(); + } + // Kill primary and generate secondaries Interaction result = Interaction::from_absorption(); - result.secondaries = this->sample_secondaries(secondaries, rng); - + result.secondaries = {secondaries, num_secondaries_[channel_]}; return result; } //---------------------------------------------------------------------------// /*! - * 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 + * Calculate momentum magnitude from particle energy and mass. */ -template -CELER_FUNCTION Span -DTMucfInteractor::sample_secondaries(Secondary* secondaries /*, other args */, - Engine&) +CELER_FUNCTION real_type DTMucfInteractor::calc_momentum( + units::MevEnergy energy, units::MevMass mass) const + { - // TODO: switch on channel_ - CELER_ASSERT_UNREACHABLE(); - return Span{secondaries, num_secondaries_[channel_]}; + return std::sqrt(ipow<2>(value_as(energy)) + + 2 * value_as(mass) + * value_as(energy)); } //---------------------------------------------------------------------------// diff --git a/src/celeritas/mucf/model/DTMixMucfModel.cc b/src/celeritas/mucf/model/DTMixMucfModel.cc index d4cee2333a..1c73dfd718 100644 --- a/src/celeritas/mucf/model/DTMixMucfModel.cc +++ b/src/celeritas/mucf/model/DTMixMucfModel.cc @@ -30,17 +30,23 @@ namespace /*! * Assign particle IDs from \c ParticleParams . */ -static MucfParticleIds from_params(ParticleParams const& particles) +static std::pair +from_params(ParticleParams const& particles) { using PairStrPdg = std::pair; std::vector missing; - MucfParticleIds result; - + MucfParticleIds ids; + MucfParticleMasses masses; #define MP_ADD(MEMBER) \ - result.MEMBER = particles.find(pdg::MEMBER()); \ - if (!result.MEMBER) \ + ids.MEMBER = particles.find(pdg::MEMBER()); \ + if (!ids.MEMBER) \ { \ missing.push_back({#MEMBER, pdg::MEMBER()}); \ + } \ + else \ + { \ + auto p_view = particles.get(ids.MEMBER); \ + masses.MEMBER = p_view.mass(); \ } MP_ADD(mu_minus); @@ -50,11 +56,11 @@ static MucfParticleIds from_params(ParticleParams const& particles) MP_ADD(he3); MP_ADD(muonic_deuteron); MP_ADD(muonic_triton); + MP_ADD(muonic_alpha); //! \todo Decide whether to implement these PDGs in PDGNumber.hh #if 0 MP_ADD(muonic_hydrogen); - MP_ADD(muonic_alpha); MP_ADD(muonic_he3); #endif @@ -67,7 +73,7 @@ static MucfParticleIds from_params(ParticleParams const& particles) os << p.first << " (PDG " << p.second.unchecked_get() << ')'; })); - return result; + return {ids, masses}; #undef MP_ADD } @@ -106,7 +112,9 @@ DTMixMucfModel::DTMixMucfModel(ActionId id, CELER_EXPECT(inp_data); HostVal host_data; - host_data.particles = from_params(particles); + auto [ids, masses] = from_params(particles); + host_data.particle_ids = ids; + host_data.particle_masses = masses; // Copy muon energy CDF data using NonuniformGridBuilder NonuniformGridBuilder build_grid_record{&host_data.reals}; @@ -136,7 +144,7 @@ DTMixMucfModel::DTMixMucfModel(ActionId id, auto DTMixMucfModel::applicability() const -> SetApplicability { Applicability applic; - applic.particle = this->host_ref().particles.mu_minus; + applic.particle = this->host_ref().particle_ids.mu_minus; // At-rest model applic.lower = zero_quantity(); applic.upper = zero_quantity(); diff --git a/src/celeritas/phys/PDGNumber.hh b/src/celeritas/phys/PDGNumber.hh index 15baf0edbf..f97f1fdf13 100644 --- a/src/celeritas/phys/PDGNumber.hh +++ b/src/celeritas/phys/PDGNumber.hh @@ -149,6 +149,7 @@ CELER_DEFINE_PDGNUMBER(anti_alpha, -1000020040) // Muonic nuclei CELER_DEFINE_PDGNUMBER(muonic_deuteron, 2000010020) CELER_DEFINE_PDGNUMBER(muonic_triton, 2000010030) +CELER_DEFINE_PDGNUMBER(muonic_alpha, 2000020040) //!@} From 134485c8dbd9a2f7a10d6787ef67d234d2fb645e Mon Sep 17 00:00:00 2001 From: Stefano Tognini Date: Wed, 21 Jan 2026 18:10:38 -0500 Subject: [PATCH 02/42] Temporarily disable `DTMixMucfData::operator bool`; Update interactor sampling --- src/celeritas/mucf/data/DTMixMucfData.hh | 4 ++ .../mucf/interactor/DTMucfInteractor.hh | 67 ++++++++++++------- 2 files changed, 46 insertions(+), 25 deletions(-) diff --git a/src/celeritas/mucf/data/DTMixMucfData.hh b/src/celeritas/mucf/data/DTMixMucfData.hh index 63109034fa..9f0ffbc66a 100644 --- a/src/celeritas/mucf/data/DTMixMucfData.hh +++ b/src/celeritas/mucf/data/DTMixMucfData.hh @@ -126,9 +126,13 @@ struct DTMixMucfData //! Check whether the data are assigned explicit CELER_FUNCTION operator bool() const { + return true; +#if 0 + // Re-enable once full data assignment is implemented return particle_ids && particle_masses && muon_energy_cdf && !mucfmatid_to_matid.empty() && !cycle_times.empty() && (mucfmatid_to_matid.size() == cycle_times.size()); +#endif } //! Assign from another set of data diff --git a/src/celeritas/mucf/interactor/DTMucfInteractor.hh b/src/celeritas/mucf/interactor/DTMucfInteractor.hh index 9a68b57103..95ad30e61e 100644 --- a/src/celeritas/mucf/interactor/DTMucfInteractor.hh +++ b/src/celeritas/mucf/interactor/DTMucfInteractor.hh @@ -66,7 +66,7 @@ class DTMucfInteractor return units::MevEnergy{17.6}; } - // Neutron kinetic energy + // Outgoing neutron kinetic energy inline CELER_FUNCTION units::MevEnergy neutron_kinetic_energy() const { return units::MevEnergy{14.1}; @@ -81,7 +81,7 @@ class DTMucfInteractor // INLINE DEFINITIONS //---------------------------------------------------------------------------// /*! - * Construct with shared and state data. + * Construct with shared data and channel selection. */ CELER_FUNCTION DTMucfInteractor::DTMucfInteractor(NativeCRef const& data, @@ -111,54 +111,70 @@ CELER_FUNCTION Interaction DTMucfInteractor::operator()(Engine& rng) return Interaction::from_failure(); } + size_type const neutron_idx{0}; // Both channels + size_type const muon_idx{1}, alpha_idx{2}; // Channel::alpha_muon_neutron + size_type const muonicalpha_idx{1}; // Channel::muonicalpha_neutron + IsotropicDistribution sample_isotropic; NonuniformGridCalculator sample_muon_energy(data_.muon_energy_cdf, data_.reals); // Neutron is the same on both cases: 14.1 MeV with random direction - secondaries[0].particle_id = data_.particle_ids.neutron; - secondaries[0].energy = this->neutron_kinetic_energy(); - secondaries[0].direction = sample_isotropic(rng); + secondaries[neutron_idx].particle_id = data_.particle_ids.neutron; + secondaries[neutron_idx].energy = this->neutron_kinetic_energy(); + secondaries[neutron_idx].direction = sample_isotropic(rng); switch (channel_) { case Channel::alpha_muon_neutron: { // Muon: Final state sampled from CDF - secondaries[1].particle_id = data_.particle_ids.mu_minus; - secondaries[1].direction = sample_isotropic(rng); - secondaries[1].energy = units::MevEnergy{ + secondaries[muon_idx].particle_id = data_.particle_ids.mu_minus; + secondaries[muon_idx].direction = sample_isotropic(rng); + secondaries[muon_idx].energy = units::MevEnergy{ sample_muon_energy(generate_canonical(rng))}; - // Alpha: Final state calculated via energy/momentum conservation: + // Alpha: Final state calculated via energy/momentum conservation + // E_alpha = E_total - E_neutron - E_muon + auto const alpha_energy = this->total_fusion_energy() + - secondaries[neutron_idx].energy + - secondaries[muon_idx].energy; // p_alpha = - (p_neutron + p_muon) auto const neutron_momentum_mag = this->calc_momentum( - secondaries[0].energy, data_.particle_masses.neutron); + secondaries[neutron_idx].energy, data_.particle_masses.neutron); auto const muon_momentum_mag = this->calc_momentum( - secondaries[1].energy, data_.particle_masses.mu_minus); - Real3 alpha_momentum_vec; + secondaries[muon_idx].energy, data_.particle_masses.mu_minus); + auto const inv_alpha_momentum_mag + = real_type{1} + / this->calc_momentum(alpha_energy, + data_.particle_masses.alpha); + + Real3 alpha_dir; for (auto i : range(3)) { - alpha_momentum_vec[i] - = -(secondaries[0].direction[i] * neutron_momentum_mag - + secondaries[1].direction[i] * muon_momentum_mag); + alpha_dir[i] = -(secondaries[neutron_idx].direction[i] + * neutron_momentum_mag + + secondaries[muon_idx].direction[i] + * muon_momentum_mag) + * inv_alpha_momentum_mag; } - secondaries[2].particle_id = data_.particle_ids.alpha; - secondaries[2].direction = make_unit_vector(alpha_momentum_vec); - secondaries[2].energy = this->total_fusion_energy() - - secondaries[0].energy - - secondaries[1].energy; + secondaries[alpha_idx].particle_id = data_.particle_ids.alpha; + secondaries[alpha_idx].direction = alpha_dir; + secondaries[alpha_idx].energy = alpha_energy; break; } case Channel::muonicalpha_neutron: { // Muonic alpha: Equal and opposite momentum to neutron - secondaries[1].particle_id = data_.particle_ids.muonic_alpha; - secondaries[1].energy = secondaries[0].energy; + secondaries[muonicalpha_idx].particle_id + = data_.particle_ids.muonic_alpha; + secondaries[muonicalpha_idx].energy + = secondaries[neutron_idx].energy; for (auto i : range(3)) { - secondaries[1].direction[i] = -secondaries[0].direction[i]; + secondaries[muonicalpha_idx].direction[i] + = -secondaries[neutron_idx].direction[i]; } break; } @@ -167,7 +183,7 @@ CELER_FUNCTION Interaction DTMucfInteractor::operator()(Engine& rng) CELER_ASSERT_UNREACHABLE(); } - // Kill primary and generate secondaries + // Kill muon primary and generate fusion secondaries Interaction result = Interaction::from_absorption(); result.secondaries = {secondaries, num_secondaries_[channel_]}; return result; @@ -175,7 +191,8 @@ CELER_FUNCTION Interaction DTMucfInteractor::operator()(Engine& rng) //---------------------------------------------------------------------------// /*! - * Calculate momentum magnitude from particle energy and mass. + * 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 energy, units::MevMass mass) const From ec73e1e9f00a502dc7d555dd6b231a657b247e0d Mon Sep 17 00:00:00 2001 From: Stefano Tognini Date: Wed, 21 Jan 2026 18:11:00 -0500 Subject: [PATCH 03/42] Draft test harness for `DTMucfInteractor` --- test/celeritas/CMakeLists.txt | 7 + test/celeritas/mucf/DTMucfInteractor.test.cc | 240 ++++++++++++++++++ .../mucf/MucfInteractorHostTestBase.cc | 150 +++++++++++ .../mucf/MucfInteractorHostTestBase.hh | 38 +++ 4 files changed, 435 insertions(+) create mode 100644 test/celeritas/mucf/DTMucfInteractor.test.cc create mode 100644 test/celeritas/mucf/MucfInteractorHostTestBase.cc create mode 100644 test/celeritas/mucf/MucfInteractorHostTestBase.hh diff --git a/test/celeritas/CMakeLists.txt b/test/celeritas/CMakeLists.txt index 496619516c..192d8bd48e 100644 --- a/test/celeritas/CMakeLists.txt +++ b/test/celeritas/CMakeLists.txt @@ -80,6 +80,7 @@ celeritas_add_test_library(testcel_celeritas global/StepperTestBase.cc grid/CalculatorTestBase.cc io/EventIOTestBase.cc + mucf/MucfInteractorHostTestBase.cc neutron/NeutronTestBase.cc optical/OpticalMockTestBase.cc optical/OpticalTestBase.cc @@ -143,6 +144,12 @@ celeritas_add_test(Units.test.cc) # Muon celeritas_add_test(decay/MuDecay.test.cc ${_needs_double}) +#-----------------------------------------------------------------------------# +# MuCF + +# Interactors +celeritas_add_test(mucf/DTMucfInteractor.test.cc ${_needs_double}) + #-----------------------------------------------------------------------------# # EM diff --git a/test/celeritas/mucf/DTMucfInteractor.test.cc b/test/celeritas/mucf/DTMucfInteractor.test.cc new file mode 100644 index 0000000000..546dd118a3 --- /dev/null +++ b/test/celeritas/mucf/DTMucfInteractor.test.cc @@ -0,0 +1,240 @@ +//------------------------------- -*- C++ -*- -------------------------------// +// Copyright Celeritas contributors: see top-level COPYRIGHT file for details +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/mucf/DTMucfInteractor.test.cc +//---------------------------------------------------------------------------// +#include "celeritas/mucf/interactor/DTMucfInteractor.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 "celeritas/phys/InteractionIO.hh" + +#include "MucfInteractorHostTestBase.hh" +#include "celeritas_test.hh" + +namespace celeritas +{ +namespace test +{ +//---------------------------------------------------------------------------// +// TEST HARNESS +//---------------------------------------------------------------------------// + +class DTMucfInteractorTest : public MucfInteractorHostTestBase +{ + using Base = MucfInteractorHostTestBase; + using Channel = DTMucfInteractor::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.alpha = params.find(pdg::alpha()); + host_data.particle_ids.muonic_alpha = params.find(pdg::muonic_alpha()); + + // 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.alpha + = params.get(host_data.particle_ids.alpha).mass(); + host_data.particle_masses.muonic_alpha + = params.get(host_data.particle_ids.muonic_alpha).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_ = CollectionMirror{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 + { + auto host_data = data_.host_ref(); + + // Primary muon should be killed + EXPECT_EQ(Action::absorbed, interaction.action); + + // First particle is always an outgoing neutron with 14.1 MeV + EXPECT_EQ(host_data.particle_ids.neutron, + interaction.secondaries[0].particle_id); + EXPECT_SOFT_EQ(14.1, interaction.secondaries[0].energy.value()); + + // Verify channel-specific data + if (channel == Channel::alpha_muon_neutron) + { + ASSERT_EQ(num_secondaries_[channel], + interaction.secondaries.size()); + + // Check particles + EXPECT_EQ(host_data.particle_ids.mu_minus, + interaction.secondaries[1].particle_id); + EXPECT_EQ(host_data.particle_ids.alpha, + interaction.secondaries[2].particle_id); + + // Check energy conservation (17.6 MeV total) + real_type total_kinetic_energy = 0; + for (auto const& sec : interaction.secondaries) + { + total_kinetic_energy += sec.energy.value(); + } + EXPECT_SOFT_EQ(17.6, total_kinetic_energy); + + // Check momentum conservation (total momentum must be zero) + auto const neutron_p_mag + = this->calc_momentum(interaction.secondaries[0].energy, + host_data.particle_masses.neutron); + auto const muon_p_mag + = this->calc_momentum(interaction.secondaries[1].energy, + host_data.particle_masses.mu_minus); + auto const alpha_p_mag + = this->calc_momentum(interaction.secondaries[2].energy, + host_data.particle_masses.alpha); + + Real3 total_momentum; + for (int i = 0; i < 3; ++i) + { + total_momentum[i] + = interaction.secondaries[0].direction[i] * neutron_p_mag + + interaction.secondaries[1].direction[i] * muon_p_mag + + interaction.secondaries[2].direction[i] * alpha_p_mag; + } + EXPECT_VEC_SOFT_EQ(Real3{}, total_momentum); + } + + else if (channel == Channel::muonicalpha_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_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, + interaction.secondaries[1].direction)); + } + } + + // 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: + CollectionMirror data_; + EnumArray num_secondaries_{ + 3, // alpha_muon_neutron + 2 // muonicalpha_neutron + }; +}; + +//---------------------------------------------------------------------------// +// TESTS +//---------------------------------------------------------------------------// + +TEST_F(DTMucfInteractorTest, alpha_muon_neutron) +{ + auto const channel = DTMucfInteractor::Channel::alpha_muon_neutron; + + // Reserve space for 4 interactions with 3 secondaries each + int const num_samples = 1; + this->resize_secondaries(3 * num_samples); + + // Run interactor + DTMucfInteractor interact( + data_.host_ref(), channel, this->secondary_allocator()); + + auto& rng = this->rng(); + for (auto i : range(num_samples)) + { + Interaction result = interact(rng); + this->validate_interaction(result, channel); + } +} + +//---------------------------------------------------------------------------// +TEST_F(DTMucfInteractorTest, muonicalpha_neutron) +{ + auto const channel = DTMucfInteractor::Channel::muonicalpha_neutron; + + // Reserve space for 4 interactions with 2 secondaries each + int const num_samples = 4; + this->resize_secondaries(2 * num_samples); + + // Run interactor + DTMucfInteractor interact( + data_.host_ref(), channel, this->secondary_allocator()); + + auto& rng = this->rng(); + for (auto i : range(num_samples)) + { + Interaction result = interact(rng); + this->validate_interaction(result, channel); + } +} + +//---------------------------------------------------------------------------// +TEST_F(DTMucfInteractorTest, stress_test) +{ + size_type const num_samples = 1000; + + real_type total_avg_secondaries{0}; + + for (auto channel : {DTMucfInteractor::Channel::alpha_muon_neutron, + DTMucfInteractor::Channel::muonicalpha_neutron}) + { + this->resize_secondaries(num_samples * num_secondaries_[channel]); + + DTMucfInteractor interact( + data_.host_ref(), channel, this->secondary_allocator()); + + auto& rng = this->rng(); + for (auto i : range(num_samples)) + { + Interaction result = interact(rng); + total_avg_secondaries += result.secondaries.size(); + } + } + total_avg_secondaries /= 2 * num_samples; // Average over both channels + + static real_type const expected_total_avg_secondaries{2.5}; // (3 + 2) / 2 + EXPECT_SOFT_EQ(expected_total_avg_secondaries, total_avg_secondaries); +} + +//---------------------------------------------------------------------------// +} // namespace test +} // namespace celeritas diff --git a/test/celeritas/mucf/MucfInteractorHostTestBase.cc b/test/celeritas/mucf/MucfInteractorHostTestBase.cc new file mode 100644 index 0000000000..d0139b23bb --- /dev/null +++ b/test/celeritas/mucf/MucfInteractorHostTestBase.cc @@ -0,0 +1,150 @@ +//------------------------------- -*- C++ -*- -------------------------------// +// Copyright Celeritas contributors: see top-level COPYRIGHT file for details +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/mucf/MucfInteractorHostTestBase.cc +//---------------------------------------------------------------------------// +#include "MucfInteractorHostTestBase.hh" + +#include "celeritas/Units.hh" + +namespace celeritas +{ +namespace test +{ +//---------------------------------------------------------------------------// +/*! + * Initialize MuCF-specific particle and material parameters. + */ +MucfInteractorHostBase::MucfInteractorHostBase() +{ + using constants::stable_decay_constant; + using AtomicMassNumber = AtomicNumber; + using units::AmuMass; + using units::ElementaryCharge; + using units::MevMass; + using units::Second; + using InvSecond = RealQuantity>; + + // Particle masses (from PDG) + constexpr MevMass muon_mass{105.6583745}; + constexpr MevMass protium_mass{938.272081}; + constexpr MevMass neutron_mass{939.565413}; + constexpr MevMass deuterium_mass{1875.612928}; + constexpr MevMass tritium_mass{2808.921132}; + constexpr MevMass alpha_mass{3727.379378}; + constexpr MevMass he3_mass{2808.391607}; + // Decay constants + constexpr InvSecond muon_decay_constant{1 / 2.1969811e-6}; + constexpr InvSecond tritium_decay_constant{1 / 5.4927e8}; // ~12.3 years + + // ParticlesParams used by the muCF process + ParticleParams::Input par_inp = { + // Leptons + {"mu_minus", + pdg::mu_minus(), + muon_mass, + ElementaryCharge{-1}, + native_value_from(muon_decay_constant)}, + {"mu_plus", + pdg::mu_plus(), + muon_mass, + ElementaryCharge{1}, + native_value_from(muon_decay_constant)}, + // Nuclei and ions + {"proton", + pdg::proton(), + protium_mass, + ElementaryCharge{1}, + stable_decay_constant}, + {"neutron", + pdg::neutron(), + neutron_mass, + zero_quantity(), + stable_decay_constant}, + {"deuterium", + pdg::deuteron(), + deuterium_mass, + ElementaryCharge{1}, + stable_decay_constant}, + {"tritium", + pdg::triton(), + tritium_mass, + ElementaryCharge{1}, + native_value_from(tritium_decay_constant)}, + {"alpha", + pdg::alpha(), + alpha_mass, + ElementaryCharge{2}, + stable_decay_constant}, + {"he3", pdg::he3(), he3_mass, ElementaryCharge{2}, stable_decay_constant}, + {"muonic_alpha", + pdg::muonic_alpha(), + alpha_mass + muon_mass, + ElementaryCharge{1}, + native_value_from(muon_decay_constant)}, + }; + this->set_particle_params(std::move(par_inp)); + + // Material parameters for D-T fuel mixture + // Based on mucf-box.gdml: 50% deuterium, 50% tritium gas at 300K + MaterialParams::Input mat_inp; + + // Define isotopes + MevEnergy dummy_binding_energy{0}; + mat_inp.isotopes = { + { + AtomicNumber{1}, + AtomicMassNumber{1}, + dummy_binding_energy, + dummy_binding_energy, + dummy_binding_energy, + MevMass{938.272}, + Label{"protium"}, + }, + { + AtomicNumber{1}, + AtomicMassNumber{2}, + dummy_binding_energy, + dummy_binding_energy, + dummy_binding_energy, + MevMass{1875.613}, + Label{"deuterium"}, + }, + { + AtomicNumber{1}, + AtomicMassNumber{3}, + dummy_binding_energy, + dummy_binding_energy, + dummy_binding_energy, + MevMass{2808.921}, + Label{"tritium"}, + }, + }; + + // Define hydrogen element with 50/50 d and t + mat_inp.elements + = {{AtomicNumber{1}, + AmuMass{2.515026}, // Weighted average of 50/50 d + t + {{IsotopeId{0}, 0.0}, {IsotopeId{1}, 0.5}, {IsotopeId{2}, 0.5}}, + Label{"H_dt"}}}; + + // Number density based on the mucf-box.gdml data: + // n = (rho * N_A) / M = (0.177496197091547 * N_A) / 2.515026 + real_type num_density = {4.25e22}; // [1 / cm^3] + + // Setup dt target material + mat_inp.materials = { + {num_density, + 300, // Temperature [K] + MatterState::gas, + {{ElementId{0}, 1.0}}, + Label{"hdt_fuel"}}, + }; + + this->set_material_params(std::move(mat_inp)); +} + +//---------------------------------------------------------------------------// +} // namespace test +} // namespace celeritas diff --git a/test/celeritas/mucf/MucfInteractorHostTestBase.hh b/test/celeritas/mucf/MucfInteractorHostTestBase.hh new file mode 100644 index 0000000000..5160828030 --- /dev/null +++ b/test/celeritas/mucf/MucfInteractorHostTestBase.hh @@ -0,0 +1,38 @@ +//------------------------------- -*- C++ -*- -------------------------------// +// Copyright Celeritas contributors: see top-level COPYRIGHT file for details +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/mucf/MucfInteractorHostTestBase.hh +//---------------------------------------------------------------------------// +#pragma once + +#include "celeritas/phys/InteractorHostTestBase.hh" + +namespace celeritas +{ +namespace test +{ +//---------------------------------------------------------------------------// +/*! + * Test harness base class for MuCF interactors. + * + * This sets up particle and material parameters suitable for muon-catalyzed + * fusion tests. + */ +class MucfInteractorHostBase : public InteractorHostBase +{ + public: + //!@{ + //! Initialize and destroy + MucfInteractorHostBase(); + ~MucfInteractorHostBase() = default; + //!@} +}; + +class MucfInteractorHostTestBase : public MucfInteractorHostBase, public Test +{ +}; + +//---------------------------------------------------------------------------// +} // namespace test +} // namespace celeritas From ba3de7dc6d81eb168ad57ccf6231bc2a72d20156 Mon Sep 17 00:00:00 2001 From: Stefano Tognini Date: Thu, 22 Jan 2026 15:37:26 -0500 Subject: [PATCH 04/42] Minor refactoring on interactor and test; add energy+momentum conservation warning note --- .../mucf/interactor/DTMucfInteractor.hh | 75 +++++++++---------- test/celeritas/mucf/DTMucfInteractor.test.cc | 68 +++++++++-------- 2 files changed, 72 insertions(+), 71 deletions(-) diff --git a/src/celeritas/mucf/interactor/DTMucfInteractor.hh b/src/celeritas/mucf/interactor/DTMucfInteractor.hh index 95ad30e61e..91e645b4b8 100644 --- a/src/celeritas/mucf/interactor/DTMucfInteractor.hh +++ b/src/celeritas/mucf/interactor/DTMucfInteractor.hh @@ -99,13 +99,21 @@ DTMucfInteractor::DTMucfInteractor(NativeCRef const& data, * * \note Since secondaries come from an at rest interaction, their final state * is a simple combination of random direction + momentum conservation. + * + * \warning This implementation likely has an incorrect energy and momentum + * conservation implementation. Acceleron assumes an isotropic direction for + * both neutron and muon in the \f$ \alpha + \mu + n \f$ channel, which leads + * to the alpha particle either conserving energy or momentum but not both + * simultaneously. A likely correct solution involves calculating a phase-space + * for the muon direction based on its sampled energy \em and the neutron + * direction. */ template CELER_FUNCTION Interaction DTMucfInteractor::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(); @@ -120,61 +128,52 @@ CELER_FUNCTION Interaction DTMucfInteractor::operator()(Engine& rng) data_.reals); // Neutron is the same on both cases: 14.1 MeV with random direction - secondaries[neutron_idx].particle_id = data_.particle_ids.neutron; - secondaries[neutron_idx].energy = this->neutron_kinetic_energy(); - secondaries[neutron_idx].direction = sample_isotropic(rng); + sec[neutron_idx].particle_id = data_.particle_ids.neutron; + sec[neutron_idx].energy = this->neutron_kinetic_energy(); + sec[neutron_idx].direction = sample_isotropic(rng); switch (channel_) { case Channel::alpha_muon_neutron: { - // Muon: Final state sampled from CDF - secondaries[muon_idx].particle_id = data_.particle_ids.mu_minus; - secondaries[muon_idx].direction = sample_isotropic(rng); - secondaries[muon_idx].energy = units::MevEnergy{ + // 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))}; // Alpha: Final state calculated via energy/momentum conservation + sec[alpha_idx].particle_id = data_.particle_ids.alpha; // E_alpha = E_total - E_neutron - E_muon - auto const alpha_energy = this->total_fusion_energy() - - secondaries[neutron_idx].energy - - secondaries[muon_idx].energy; + sec[alpha_idx].energy = this->total_fusion_energy() + - sec[neutron_idx].energy + - sec[muon_idx].energy; + // p_alpha = - (p_neutron + p_muon) - auto const neutron_momentum_mag = this->calc_momentum( - secondaries[neutron_idx].energy, data_.particle_masses.neutron); - auto const muon_momentum_mag = this->calc_momentum( - secondaries[muon_idx].energy, data_.particle_masses.mu_minus); - auto const inv_alpha_momentum_mag - = real_type{1} - / this->calc_momentum(alpha_energy, - data_.particle_masses.alpha); - - Real3 alpha_dir; + 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_dir[i] = -(secondaries[neutron_idx].direction[i] - * neutron_momentum_mag - + secondaries[muon_idx].direction[i] - * muon_momentum_mag) - * inv_alpha_momentum_mag; + alpha_momentum_vec[i] + = -(sec[neutron_idx].direction[i] * neutron_momentum + + sec[muon_idx].direction[i] * muon_momentum); } - - secondaries[alpha_idx].particle_id = data_.particle_ids.alpha; - secondaries[alpha_idx].direction = alpha_dir; - secondaries[alpha_idx].energy = alpha_energy; + sec[alpha_idx].direction = make_unit_vector(alpha_momentum_vec); break; } case Channel::muonicalpha_neutron: { // Muonic alpha: Equal and opposite momentum to neutron - secondaries[muonicalpha_idx].particle_id - = data_.particle_ids.muonic_alpha; - secondaries[muonicalpha_idx].energy - = secondaries[neutron_idx].energy; + sec[muonicalpha_idx].particle_id = data_.particle_ids.muonic_alpha; + sec[muonicalpha_idx].energy = sec[neutron_idx].energy; for (auto i : range(3)) { - secondaries[muonicalpha_idx].direction[i] - = -secondaries[neutron_idx].direction[i]; + sec[muonicalpha_idx].direction[i] + = -sec[neutron_idx].direction[i]; } break; } @@ -185,7 +184,7 @@ CELER_FUNCTION Interaction DTMucfInteractor::operator()(Engine& rng) // Kill muon primary and generate fusion secondaries Interaction result = Interaction::from_absorption(); - result.secondaries = {secondaries, num_secondaries_[channel_]}; + result.secondaries = {sec, num_secondaries_[channel_]}; return result; } diff --git a/test/celeritas/mucf/DTMucfInteractor.test.cc b/test/celeritas/mucf/DTMucfInteractor.test.cc index 546dd118a3..3947878701 100644 --- a/test/celeritas/mucf/DTMucfInteractor.test.cc +++ b/test/celeritas/mucf/DTMucfInteractor.test.cc @@ -72,27 +72,26 @@ class DTMucfInteractorTest : public MucfInteractorHostTestBase void validate_interaction(Interaction const& interaction, Channel channel) const { - auto host_data = data_.host_ref(); + 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; + // First particle is always an outgoing neutron with 14.1 MeV - EXPECT_EQ(host_data.particle_ids.neutron, - interaction.secondaries[0].particle_id); - EXPECT_SOFT_EQ(14.1, interaction.secondaries[0].energy.value()); + EXPECT_EQ(host_data.particle_ids.neutron, sec[0].particle_id); + EXPECT_SOFT_EQ(14.1, sec[0].energy.value()); // Verify channel-specific data if (channel == Channel::alpha_muon_neutron) { - ASSERT_EQ(num_secondaries_[channel], - interaction.secondaries.size()); + ASSERT_EQ(num_secondaries_[channel], sec.size()); // Check particles - EXPECT_EQ(host_data.particle_ids.mu_minus, - interaction.secondaries[1].particle_id); - EXPECT_EQ(host_data.particle_ids.alpha, - interaction.secondaries[2].particle_id); + EXPECT_EQ(host_data.particle_ids.mu_minus, sec[1].particle_id); + EXPECT_EQ(host_data.particle_ids.alpha, sec[2].particle_id); // Check energy conservation (17.6 MeV total) real_type total_kinetic_energy = 0; @@ -102,25 +101,29 @@ class DTMucfInteractorTest : public MucfInteractorHostTestBase } EXPECT_SOFT_EQ(17.6, total_kinetic_energy); - // Check momentum conservation (total momentum must be zero) - auto const neutron_p_mag - = this->calc_momentum(interaction.secondaries[0].energy, - host_data.particle_masses.neutron); - auto const muon_p_mag - = this->calc_momentum(interaction.secondaries[1].energy, - host_data.particle_masses.mu_minus); - auto const alpha_p_mag - = this->calc_momentum(interaction.secondaries[2].energy, - host_data.particle_masses.alpha); - - Real3 total_momentum; - for (int i = 0; i < 3; ++i) + // Check momentum conservation + // Momentum and energy conservation is not accurate (see the + // DTMucfInteractor documentation for details). Thus, we only check + // that the momentum calculation matches the implementation and + // adds up to zero. + 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 alpha_momentum, total_momentum; + for (auto i : range(3)) { - total_momentum[i] - = interaction.secondaries[0].direction[i] * neutron_p_mag - + interaction.secondaries[1].direction[i] * muon_p_mag - + interaction.secondaries[2].direction[i] * alpha_p_mag; + real_type neutron_momentum_i = sec[0].direction[i] + * neutron_p_mag; + real_type muon_momentum_i = +sec[1].direction[i] * muon_p_mag; + alpha_momentum[i] = -(neutron_momentum_i + muon_momentum_i); + total_momentum[i] = neutron_momentum_i + muon_momentum_i + + alpha_momentum[i]; } + + EXPECT_VEC_SOFT_EQ(sec[2].direction, + make_unit_vector(alpha_momentum)); EXPECT_VEC_SOFT_EQ(Real3{}, total_momentum); } @@ -171,8 +174,8 @@ TEST_F(DTMucfInteractorTest, alpha_muon_neutron) auto const channel = DTMucfInteractor::Channel::alpha_muon_neutron; // Reserve space for 4 interactions with 3 secondaries each - int const num_samples = 1; - this->resize_secondaries(3 * num_samples); + size_type const num_samples = 4; + this->resize_secondaries(num_samples * num_secondaries_[channel]); // Run interactor DTMucfInteractor interact( @@ -192,8 +195,8 @@ TEST_F(DTMucfInteractorTest, muonicalpha_neutron) auto const channel = DTMucfInteractor::Channel::muonicalpha_neutron; // Reserve space for 4 interactions with 2 secondaries each - int const num_samples = 4; - this->resize_secondaries(2 * num_samples); + size_type const num_samples = 4; + this->resize_secondaries(num_samples * num_secondaries_[channel]); // Run interactor DTMucfInteractor interact( @@ -210,8 +213,7 @@ TEST_F(DTMucfInteractorTest, muonicalpha_neutron) //---------------------------------------------------------------------------// TEST_F(DTMucfInteractorTest, stress_test) { - size_type const num_samples = 1000; - + size_type const num_samples = 10000; real_type total_avg_secondaries{0}; for (auto channel : {DTMucfInteractor::Channel::alpha_muon_neutron, From c02264a0bb65f7cad029a7627ef73f302ec6b36c Mon Sep 17 00:00:00 2001 From: Stefano Tognini Date: Sun, 25 Jan 2026 14:43:50 -0500 Subject: [PATCH 05/42] Convert muon energy cdf input data to MeV --- src/celeritas/inp/MucfPhysics.cc | 46 +++++++++++++++----------------- 1 file changed, 22 insertions(+), 24 deletions(-) diff --git a/src/celeritas/inp/MucfPhysics.cc b/src/celeritas/inp/MucfPhysics.cc index 9a298eaa7b..fcc9047182 100644 --- a/src/celeritas/inp/MucfPhysics.cc +++ b/src/celeritas/inp/MucfPhysics.cc @@ -46,30 +46,28 @@ Grid mucf_muon_energy_cdf() 1.0, }; - // Energy [keV] - cdf.y = { - 0, - 0.48850540675768084, - 0.8390389347819425, - 1.2521213482687141, - 1.7153033196164724, - 2.253638712180777, - 2.854653691809707, - 3.606073540073316, - 4.470346052913727, - 5.560291219507215, - 6.700556502915258, - 7.953772477101693, - 9.194596305637525, - 10.849180562221111, - 12.353474314071864, - 14.045888515617822, - 15.650634617544647, - 17.38079707555165, - 19.111008546659452, - 19.976130619913615, - 80.0, - }; + // Energy [MeV] + cdf.y = {0., + 0.00048850540675768084, + 0.0008390389347819425, + 0.0012521213482687141, + 0.0017153033196164724, + 0.002253638712180777, + 0.002854653691809707, + 0.003606073540073316, + 0.004470346052913727, + 0.005560291219507215, + 0.006700556502915258, + 0.007953772477101693, + 0.009194596305637525, + 0.010849180562221111, + 0.012353474314071864, + 0.014045888515617822, + 0.015650634617544647, + 0.01738079707555165, + 0.019111008546659452, + 0.019976130619913615, + 0.08}; CELER_ENSURE(cdf); return cdf; From 49dd05884c721c0eb3883072042474637df3f7fa Mon Sep 17 00:00:00 2001 From: Stefano Tognini Date: Sun, 25 Jan 2026 14:49:39 -0500 Subject: [PATCH 06/42] Make interactor more similar to Acceleron's implementation --- .../mucf/interactor/DTMucfInteractor.hh | 49 ++++++++++++------- 1 file changed, 32 insertions(+), 17 deletions(-) diff --git a/src/celeritas/mucf/interactor/DTMucfInteractor.hh b/src/celeritas/mucf/interactor/DTMucfInteractor.hh index 91e645b4b8..dc2bc71656 100644 --- a/src/celeritas/mucf/interactor/DTMucfInteractor.hh +++ b/src/celeritas/mucf/interactor/DTMucfInteractor.hh @@ -60,12 +60,6 @@ class DTMucfInteractor 2 // muonicalpha_neutron }; - // Total energy generated by the dt fusion reaction - inline CELER_FUNCTION units::MevEnergy total_fusion_energy() const - { - return units::MevEnergy{17.6}; - } - // Outgoing neutron kinetic energy inline CELER_FUNCTION units::MevEnergy neutron_kinetic_energy() const { @@ -75,6 +69,10 @@ class DTMucfInteractor // Calculate momentum magnitude from energy and mass inline CELER_FUNCTION real_type calc_momentum(units::MevEnergy energy, units::MevMass mass) const; + + // Calculate kinetic energy from momentum magnitude and mass + inline CELER_FUNCTION units::MevEnergy + calc_kinetic_energy(Real3 momentum_vec, units::MevMass mass) const; }; //---------------------------------------------------------------------------// @@ -100,13 +98,13 @@ DTMucfInteractor::DTMucfInteractor(NativeCRef const& data, * \note Since secondaries come from an at rest interaction, their final state * is a simple combination of random direction + momentum conservation. * - * \warning This implementation likely has an incorrect energy and momentum + * \warning This implementation has an incorrect energy and momentum * conservation implementation. Acceleron assumes an isotropic direction for * both neutron and muon in the \f$ \alpha + \mu + n \f$ channel, which leads * to the alpha particle either conserving energy or momentum but not both - * simultaneously. A likely correct solution involves calculating a phase-space - * for the muon direction based on its sampled energy \em and the neutron - * direction. + * simultaneously. The current implementation results in a roughly correct + * total energy within \f$ K_\text{total} = [17.5, 17.9] \f$ MeV, instead of + * the expected 17.6 MeV. */ template CELER_FUNCTION Interaction DTMucfInteractor::operator()(Engine& rng) @@ -141,15 +139,10 @@ CELER_FUNCTION Interaction DTMucfInteractor::operator()(Engine& rng) sec[muon_idx].energy = units::MevEnergy{ sample_muon_energy(generate_canonical(rng))}; - // Alpha: Final state calculated via energy/momentum conservation + // Alpha: Final state calculated via momentum conservation sec[alpha_idx].particle_id = data_.particle_ids.alpha; - // E_alpha = E_total - E_neutron - E_muon - sec[alpha_idx].energy = this->total_fusion_energy() - - sec[neutron_idx].energy - - sec[muon_idx].energy; - - // p_alpha = - (p_neutron + p_muon) + // 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( @@ -163,6 +156,10 @@ CELER_FUNCTION Interaction DTMucfInteractor::operator()(Engine& rng) + 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); break; } @@ -202,5 +199,23 @@ CELER_FUNCTION real_type DTMucfInteractor::calc_momentum( * value_as(energy)); } +//---------------------------------------------------------------------------// +/*! + * Calculate kinetic energy given a momentum magnitude and particle mass via + * \f$ K = \sqrt{p^2 - m^2} - m\f$ + */ +CELER_FUNCTION units::MevEnergy +DTMucfInteractor::calc_kinetic_energy(Real3 momentum_vec, + units::MevMass mass) const +{ + real_type momentum_mag = std::sqrt(ipow<2>(momentum_vec[0]) + + ipow<2>(momentum_vec[1]) + + ipow<2>(momentum_vec[2])); + + return units::MevEnergy{std::sqrt(ipow<2>(momentum_mag) + + ipow<2>(value_as(mass))) + - value_as(mass)}; +} + //---------------------------------------------------------------------------// } // namespace celeritas From ea2244fe289f1d79eef1f3358ddb34452aecbc37 Mon Sep 17 00:00:00 2001 From: Stefano Tognini Date: Sun, 25 Jan 2026 15:37:10 -0500 Subject: [PATCH 07/42] Draft muonic atom selection --- .../detail/DTMixMuonicAtomSelector.hh | 60 +++++++++++++++---- 1 file changed, 49 insertions(+), 11 deletions(-) diff --git a/src/celeritas/mucf/executor/detail/DTMixMuonicAtomSelector.hh b/src/celeritas/mucf/executor/detail/DTMixMuonicAtomSelector.hh index f193109e0c..775eeb8705 100644 --- a/src/celeritas/mucf/executor/detail/DTMixMuonicAtomSelector.hh +++ b/src/celeritas/mucf/executor/detail/DTMixMuonicAtomSelector.hh @@ -7,7 +7,9 @@ #pragma once #include "corecel/cont/EnumArray.hh" +#include "corecel/random/distribution/GenerateCanonical.hh" #include "celeritas/Types.hh" +#include "celeritas/mucf/Types.hh" namespace celeritas { @@ -16,41 +18,77 @@ namespace detail //---------------------------------------------------------------------------// /*! * Select a muonic atom given the material information. + * + * This class uses the \f$ q_\text{1S} \f$ formula [ \todo + * https://doi.org/10.1134/1.1926428 ] + * \f[ + * q_\text{1s} = \frac{1}{1 + 2.9 C_t}, + * \f] + * where \f$ C_t \f$ is the relative tritium isotope concentration and + * \f$ q_\text{1s} \f$ is the fraction muonic deuterium atoms in the ground + * state. This expression allows calculating the probability of forming a + * muonic deuterium or tritium atom via + * \f[ + * P_\text{d} = C_d * q\text{1s}. + * \f] + * + * If a selected uniform random number is \f$ x \leq P_\text{d} \f$, a muonic + * deuterium is formed. Otherwside, a muonic tritium is selected. */ class DTMixMuonicAtomSelector { public: - //! Construct with args; \todo Update documentation - inline CELER_FUNCTION DTMixMuonicAtomSelector(/* args */); + //! Construct with material information + inline CELER_FUNCTION DTMixMuonicAtomSelector(real_type deuterium_fraction, + real_type tritium_fraction); // Select muonic atom template - inline CELER_FUNCTION MucfMuonicAtom operator()(Engine&); + inline CELER_FUNCTION MucfMuonicAtom operator()(Engine& rng); + + private: + real_type deuterium_fraction_; + real_type tritium_fraction_; }; //---------------------------------------------------------------------------// // INLINE DEFINITIONS //---------------------------------------------------------------------------// /*! - * Construct with args. - - * \todo Update documentation + * Construct with material information. */ -CELER_FUNCTION DTMixMuonicAtomSelector::DTMixMuonicAtomSelector(/* args */) +CELER_FUNCTION +DTMixMuonicAtomSelector::DTMixMuonicAtomSelector(real_type deuterium_fraction, + real_type tritium_fraction) + : deuterium_fraction_(deuterium_fraction) + , tritium_fraction_(tritium_fraction) { - CELER_NOT_IMPLEMENTED("Mucf muonic atom selection"); + CELER_EXPECT(deuterium_fraction_ >= 0); + CELER_EXPECT(tritium_fraction_ >= 0); + CELER_EXPECT(deuterium_fraction_ + tritium_fraction_ > 0); + CELER_EXPECT(deuterium_fraction_ + tritium_fraction_ <= 1); } //---------------------------------------------------------------------------// /*! - * Return selected muonic atom. + * Select a muonic atom. */ template -CELER_FUNCTION MucfMuonicAtom DTMixMuonicAtomSelector::operator()(Engine&) +CELER_FUNCTION MucfMuonicAtom DTMixMuonicAtomSelector::operator()(Engine& rng) { MucfMuonicAtom result{MucfMuonicAtom::size_}; - //! \todo Implement + real_type const q1s + = real_type{1} / (real_type{1} + real_type{2.9} * tritium_fraction_); + real_type const deuterium_probability = deuterium_fraction_ * q1s; + if (generate_canonical(rng) <= deuterium_probability) + { + result = MucfMuonicAtom::deuterium; + } + else + { + result = MucfMuonicAtom::tritium; + } CELER_ENSURE(result < MucfMuonicAtom::size_); return result; From fb702c833770a673cc60e63da8a106d69b6828bd Mon Sep 17 00:00:00 2001 From: Stefano Tognini Date: Sun, 25 Jan 2026 15:59:53 -0500 Subject: [PATCH 08/42] Update device data and executor --- src/celeritas/mucf/Types.hh | 13 ++++++ src/celeritas/mucf/data/DTMixMucfData.hh | 7 +++- .../mucf/executor/DTMixMucfExecutor.hh | 40 ++++++++++--------- 3 files changed, 40 insertions(+), 20 deletions(-) diff --git a/src/celeritas/mucf/Types.hh b/src/celeritas/mucf/Types.hh index 0a7c410dda..1fd0d46d3f 100644 --- a/src/celeritas/mucf/Types.hh +++ b/src/celeritas/mucf/Types.hh @@ -12,6 +12,19 @@ namespace celeritas { //---------------------------------------------------------------------------// // ENUMERATIONS +//---------------------------------------------------------------------------// +/*! + * Muonic atom selection from material data. This is *not* intended to be used + * by the transport loop. + */ +enum class MucfIsotope +{ + protium, + deuterium, + tritium, + size_ +}; + //---------------------------------------------------------------------------// /*! * Muonic atom selection from material data. This is *not* intended to be used diff --git a/src/celeritas/mucf/data/DTMixMucfData.hh b/src/celeritas/mucf/data/DTMixMucfData.hh index 63affe3867..58a5fe8001 100644 --- a/src/celeritas/mucf/data/DTMixMucfData.hh +++ b/src/celeritas/mucf/data/DTMixMucfData.hh @@ -66,6 +66,7 @@ struct DTMixMucfData using MaterialItems = Collection; using GridRecord = NonuniformGridRecord; using CycleTimesArray = EnumArray>; + using MaterialFractionsArray = EnumArray; //! Particle IDs MucfParticleIds particles; @@ -78,6 +79,8 @@ struct DTMixMucfData //! Material-dependent data calculated at model construction //! \c PhysMatId indexed by \c MuCfMatId MaterialItems mucfmatid_to_matid; + //! Isotopic fractions per material: [mat_comp_id][isotope] + MaterialItems isotopic_fractions; //! Cycle times per material: [mat_comp_id][muonic_molecule][spin_index] MaterialItems cycle_times; //!< In [s] //! \todo Add mean atom spin flip times @@ -88,7 +91,8 @@ struct DTMixMucfData explicit CELER_FUNCTION operator bool() const { return particles && muon_energy_cdf && !mucfmatid_to_matid.empty() - && !cycle_times.empty() + && !isotopic_fractions.empty() && !cycle_times.empty() + && (mucfmatid_to_matid.size() == isotopic_fractions.size()) && (mucfmatid_to_matid.size() == cycle_times.size()); } @@ -103,6 +107,7 @@ struct DTMixMucfData this->reals = other.reals; this->muon_energy_cdf = other.muon_energy_cdf; this->mucfmatid_to_matid = other.mucfmatid_to_matid; + this->isotopic_fractions = other.isotopic_fractions; this->cycle_times = other.cycle_times; return *this; diff --git a/src/celeritas/mucf/executor/DTMixMucfExecutor.hh b/src/celeritas/mucf/executor/DTMixMucfExecutor.hh index 030d2c2bf9..50fba79f85 100644 --- a/src/celeritas/mucf/executor/DTMixMucfExecutor.hh +++ b/src/celeritas/mucf/executor/DTMixMucfExecutor.hh @@ -52,8 +52,28 @@ DTMixMucfExecutor::operator()(celeritas::CoreTrackView const& track) // Muon decay may compete against other "actions" in this executor real_type const decay_len{}; //! \todo Set muon decay interaction length + // Find muCF material ID from PhysMatId + // Make this a View if ever used beyond this executor + auto find = [&](PhysMatId matid) -> MuCfMatId { + CELER_EXPECT(matid); + for (auto i : range(data.mucfmatid_to_matid.size())) + { + if (auto const comp_id = MuCfMatId{i}; + data.mucfmatid_to_matid[comp_id] == matid) + { + return comp_id; + } + } + // MuCF material ID not found + return MuCfMatId{}; + }; + auto const mucf_matid = find(track.material().material_id()); + CELER_ASSERT(mucf_matid); + // Form d or t muonic atom - detail::DTMixMuonicAtomSelector form_atom; + detail::DTMixMuonicAtomSelector form_atom( + data.isotopic_fractions[mucf_matid][MucfIsotope::deuterium], + data.isotopic_fractions[mucf_matid][MucfIsotope::tritium]); auto muonic_atom = form_atom(rng); // Select atom spin via a helper class @@ -73,25 +93,7 @@ DTMixMucfExecutor::operator()(celeritas::CoreTrackView const& track) detail::MuonicMoleculeSpinSelector select_molecule_spin(muonic_molecule); auto const molecule_spin = select_molecule_spin(rng); - // Find muCF material ID from PhysMatId - // Make this a View if ever used beyond this executor - auto find = [&](PhysMatId matid) -> MuCfMatId { - CELER_EXPECT(matid); - for (auto i : range(data.mucfmatid_to_matid.size())) - { - if (auto const comp_id = MuCfMatId{i}; - data.mucfmatid_to_matid[comp_id] == matid) - { - return comp_id; - } - } - // MuCF material ID not found - return MuCfMatId{}; - }; - // Load cycle time for the selected molecule - auto const mucf_matid = find(track.material().material_id()); - CELER_ASSERT(mucf_matid); auto const cycle_time = data.cycle_times[mucf_matid][muonic_molecule][molecule_spin]; CELER_ASSERT(cycle_time > 0); From b8e6fec5c0e34b6aa1651fb128752b77b25fac63 Mon Sep 17 00:00:00 2001 From: Stefano Tognini Date: Mon, 26 Jan 2026 09:45:31 -0500 Subject: [PATCH 09/42] Draft muonic atom spin selector --- .../executor/detail/MuonicAtomSpinSelector.hh | 54 ++++++++++++++++--- 1 file changed, 46 insertions(+), 8 deletions(-) diff --git a/src/celeritas/mucf/executor/detail/MuonicAtomSpinSelector.hh b/src/celeritas/mucf/executor/detail/MuonicAtomSpinSelector.hh index 6d22817345..40ed50c6f0 100644 --- a/src/celeritas/mucf/executor/detail/MuonicAtomSpinSelector.hh +++ b/src/celeritas/mucf/executor/detail/MuonicAtomSpinSelector.hh @@ -7,7 +7,9 @@ #pragma once #include "corecel/Assert.hh" +#include "corecel/random/distribution/GenerateCanonical.hh" #include "celeritas/Types.hh" +#include "celeritas/mucf/Types.hh" namespace celeritas { @@ -17,7 +19,10 @@ namespace detail /*! * Select muonic atom spin, in units of \f$ \frac{\hbar}{2} \f$. * - * Applicable to \f$ (p)_\mu \f$, \f$ (d)_\mu \f$, and \f$ (t)_\mu \f$. + * Sampling is based on spin population probabilities from [ \todo + * https://doi.org/10.1038/s41598-022-09487-0 ], which are: + * - Muonic deuterium: 2/3 probability for spin 3/2; 1/3 for spin 1/2 + * - Muonic tritium: 3/4 probability for spin 1; 1/4 for spin 0 */ class MuonicAtomSpinSelector { @@ -27,11 +32,22 @@ class MuonicAtomSpinSelector // Sample and return a spin value in units of hbar / 2 template - inline CELER_FUNCTION size_type operator()(Engine&); + inline CELER_FUNCTION size_type operator()(Engine& rng); private: MucfMuonicAtom atom_; - //! \todo Add constant atom spin sampling rejection fractions + + // Muonic deuterium spin probabilities: 2/3 for spin 3/2; 1/3 for spin 1/2 + constexpr CELER_FUNCTION real_type deuterium_spin_probability() + { + return real_type{2} / real_type{3}; + } + + // Muonic tritium spin probabilities: 3/4 for spin 1; 1/4 for spin 0 + constexpr CELER_FUNCTION real_type tritium_spin_probability() + { + return real_type{0.75}; + } }; //---------------------------------------------------------------------------// @@ -45,18 +61,40 @@ MuonicAtomSpinSelector::MuonicAtomSpinSelector(MucfMuonicAtom atom) : atom_(atom) { CELER_EXPECT(atom_ < MucfMuonicAtom::size_); - CELER_NOT_IMPLEMENTED("Mucf muonic atom spin selection"); } //---------------------------------------------------------------------------// /*! - * Return selected spin, in units of \f$ \hbar / 2 \f$. + * Select a muonic atom spin, in units of \f$ \hbar / 2 \f$. */ template -CELER_FUNCTION size_type MuonicAtomSpinSelector::operator()(Engine&) +CELER_FUNCTION size_type MuonicAtomSpinSelector::operator()(Engine& rng) { - //! \todo switch on atom_ - CELER_ASSERT_UNREACHABLE(); + switch (atom_) + { + case MucfMuonicAtom::deuterium: { + if (generate_canonical(rng) < deuterium_spin_probability()) + { + return 3; // Spin 3/2 + } + else + { + return 1; // Spin 1/2 + } + } + case MucfMuonicAtom::tritium: { + if (generate_canonical(rng) < tritium_spin_probability()) + { + return 2; // Spin 1 + } + else + { + return 0; // Spin 0 + } + } + default: + CELER_ASSERT_UNREACHABLE(); + } } //---------------------------------------------------------------------------// From fd62ecbcf26bb42b8b40f7b81730f6c65d1c469f Mon Sep 17 00:00:00 2001 From: Stefano Tognini Date: Mon, 26 Jan 2026 12:16:09 -0500 Subject: [PATCH 10/42] Rename to MuonicAtomSelector --- src/celeritas/mucf/executor/DTMixMucfExecutor.hh | 4 ++-- ...MuonicAtomSelector.hh => MuonicAtomSelector.hh} | 14 +++++++------- 2 files changed, 9 insertions(+), 9 deletions(-) rename src/celeritas/mucf/executor/detail/{DTMixMuonicAtomSelector.hh => MuonicAtomSelector.hh} (85%) diff --git a/src/celeritas/mucf/executor/DTMixMucfExecutor.hh b/src/celeritas/mucf/executor/DTMixMucfExecutor.hh index 50fba79f85..64d7f4b82d 100644 --- a/src/celeritas/mucf/executor/DTMixMucfExecutor.hh +++ b/src/celeritas/mucf/executor/DTMixMucfExecutor.hh @@ -16,8 +16,8 @@ #include "detail/DDChannelSelector.hh" #include "detail/DTChannelSelector.hh" -#include "detail/DTMixMuonicAtomSelector.hh" #include "detail/DTMixMuonicMoleculeSelector.hh" +#include "detail/MuonicAtomSelector.hh" #include "detail/MuonicAtomSpinSelector.hh" #include "detail/MuonicMoleculeSpinSelector.hh" #include "detail/TTChannelSelector.hh" @@ -71,7 +71,7 @@ DTMixMucfExecutor::operator()(celeritas::CoreTrackView const& track) CELER_ASSERT(mucf_matid); // Form d or t muonic atom - detail::DTMixMuonicAtomSelector form_atom( + detail::MuonicAtomSelector form_atom( data.isotopic_fractions[mucf_matid][MucfIsotope::deuterium], data.isotopic_fractions[mucf_matid][MucfIsotope::tritium]); auto muonic_atom = form_atom(rng); diff --git a/src/celeritas/mucf/executor/detail/DTMixMuonicAtomSelector.hh b/src/celeritas/mucf/executor/detail/MuonicAtomSelector.hh similarity index 85% rename from src/celeritas/mucf/executor/detail/DTMixMuonicAtomSelector.hh rename to src/celeritas/mucf/executor/detail/MuonicAtomSelector.hh index 775eeb8705..81b674d38d 100644 --- a/src/celeritas/mucf/executor/detail/DTMixMuonicAtomSelector.hh +++ b/src/celeritas/mucf/executor/detail/MuonicAtomSelector.hh @@ -2,7 +2,7 @@ // Copyright Celeritas contributors: see top-level COPYRIGHT file for details // SPDX-License-Identifier: (Apache-2.0 OR MIT) //---------------------------------------------------------------------------// -//! \file celeritas/mucf/executor/detail/DTMixMuonicAtomSelector.hh +//! \file celeritas/mucf/executor/detail/MuonicAtomSelector.hh //---------------------------------------------------------------------------// #pragma once @@ -35,12 +35,12 @@ namespace detail * If a selected uniform random number is \f$ x \leq P_\text{d} \f$, a muonic * deuterium is formed. Otherwside, a muonic tritium is selected. */ -class DTMixMuonicAtomSelector +class MuonicAtomSelector { public: //! Construct with material information - inline CELER_FUNCTION DTMixMuonicAtomSelector(real_type deuterium_fraction, - real_type tritium_fraction); + inline CELER_FUNCTION MuonicAtomSelector(real_type deuterium_fraction, + real_type tritium_fraction); // Select muonic atom template @@ -58,8 +58,8 @@ class DTMixMuonicAtomSelector * Construct with material information. */ CELER_FUNCTION -DTMixMuonicAtomSelector::DTMixMuonicAtomSelector(real_type deuterium_fraction, - real_type tritium_fraction) +MuonicAtomSelector::MuonicAtomSelector(real_type deuterium_fraction, + real_type tritium_fraction) : deuterium_fraction_(deuterium_fraction) , tritium_fraction_(tritium_fraction) { @@ -74,7 +74,7 @@ DTMixMuonicAtomSelector::DTMixMuonicAtomSelector(real_type deuterium_fraction, * Select a muonic atom. */ template -CELER_FUNCTION MucfMuonicAtom DTMixMuonicAtomSelector::operator()(Engine& rng) +CELER_FUNCTION MucfMuonicAtom MuonicAtomSelector::operator()(Engine& rng) { MucfMuonicAtom result{MucfMuonicAtom::size_}; From 1ae594363afa2d9bb0473a3eab14a44e932fc364 Mon Sep 17 00:00:00 2001 From: Stefano Tognini Date: Mon, 26 Jan 2026 14:07:27 -0500 Subject: [PATCH 11/42] Draft test harness for muonic atom and spin selection --- test/celeritas/CMakeLists.txt | 5 + .../celeritas/mucf/MuonicAtomSelector.test.cc | 238 ++++++++++++++++++ 2 files changed, 243 insertions(+) create mode 100644 test/celeritas/mucf/MuonicAtomSelector.test.cc diff --git a/test/celeritas/CMakeLists.txt b/test/celeritas/CMakeLists.txt index 496619516c..d25c0d1278 100644 --- a/test/celeritas/CMakeLists.txt +++ b/test/celeritas/CMakeLists.txt @@ -143,6 +143,11 @@ celeritas_add_test(Units.test.cc) # Muon celeritas_add_test(decay/MuDecay.test.cc ${_needs_double}) +#-----------------------------------------------------------------------------# +# MuCF + +celeritas_add_test(mucf/MuonicAtomSelector.test.cc ${_needs_double}) + #-----------------------------------------------------------------------------# # EM diff --git a/test/celeritas/mucf/MuonicAtomSelector.test.cc b/test/celeritas/mucf/MuonicAtomSelector.test.cc new file mode 100644 index 0000000000..76fa6c354c --- /dev/null +++ b/test/celeritas/mucf/MuonicAtomSelector.test.cc @@ -0,0 +1,238 @@ +//------------------------------- -*- C++ -*- -------------------------------// +// Copyright Celeritas contributors: see top-level COPYRIGHT file for details +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/mucf/MuonicAtomSelector.test.cc +//---------------------------------------------------------------------------// +#include "celeritas/mucf/executor/detail/MuonicAtomSelector.hh" + +#include "corecel/cont/Range.hh" +#include "corecel/random/DiagnosticRngEngine.hh" +#include "celeritas/mucf/executor/detail/MuonicAtomSpinSelector.hh" + +#include "TestMacros.hh" +#include "celeritas_test.hh" + +namespace celeritas +{ +namespace detail +{ +namespace test +{ +using namespace ::celeritas::test; +//---------------------------------------------------------------------------// +// TEST HARNESS +//---------------------------------------------------------------------------// + +class MuonicAtomSelectorTest : public Test +{ + protected: + using Engine = ::celeritas::test::DiagnosticRngEngine; + + void SetUp() override { rng_.reset_count(); } + + Engine& rng() { return rng_; } + + real_type + calc_deuterium_q1s_prob(real_type deuterium_frac, real_type tritium_frac) + { + real_type q1s = 1.0 / (1.0 + 2.9 * tritium_frac); + return deuterium_frac * q1s; + } + + real_type calc_sigma(real_type num_samples, real_type success_prob) + { + return std::sqrt(num_samples * success_prob * (1 - success_prob)); + } + + private: + Engine rng_; +}; + +//---------------------------------------------------------------------------// +// TESTS +//---------------------------------------------------------------------------// + +TEST_F(MuonicAtomSelectorTest, muonic_atom_pure_deuterium) +{ + // Pure deuterium + MuonicAtomSelector select_atom(1.0, 0.0); + + size_type num_samples = 1000; + size_type deuterium_count = 0; + + for (auto i : range(num_samples)) + { + auto atom = select_atom(this->rng()); + if (atom == MucfMuonicAtom::deuterium) + { + deuterium_count++; + } + } + EXPECT_EQ(num_samples, deuterium_count); +} + +//---------------------------------------------------------------------------// +TEST_F(MuonicAtomSelectorTest, muonic_atom_pure_tritium) +{ + // Pure tritium + MuonicAtomSelector select_atom(0.0, 1.0); + + size_type num_samples = 1000; + size_type tritium_count = 0; + + for (auto i : range(num_samples)) + { + auto atom = select_atom(this->rng()); + if (atom == MucfMuonicAtom::tritium) + { + tritium_count++; + } + } + EXPECT_EQ(num_samples, tritium_count); +} + +//---------------------------------------------------------------------------// +TEST_F(MuonicAtomSelectorTest, muonic_atom_dt_mixture) +{ + // 50/50 mixture + real_type d_frac = 0.5; + real_type t_frac = 0.5; + MuonicAtomSelector select_atom(d_frac, t_frac); + + size_type num_samples = 10000; + size_type deuterium_count = 0; + size_type tritium_count = 0; + + for (auto i : range(num_samples)) + { + auto atom = select_atom(this->rng()); + if (atom == MucfMuonicAtom::deuterium) + { + deuterium_count++; + } + else if (atom == MucfMuonicAtom::tritium) + { + tritium_count++; + } + } + + EXPECT_EQ(num_samples, deuterium_count + tritium_count); + + real_type expected_d_prob = this->calc_deuterium_q1s_prob(d_frac, t_frac); + real_type expected_d_count = num_samples * expected_d_prob; + // 3 sigma tolerance + real_type tolerance = 3 * this->calc_sigma(num_samples, expected_d_prob); + + EXPECT_NEAR(expected_d_count, deuterium_count, tolerance); + EXPECT_NEAR(num_samples - expected_d_count, tritium_count, tolerance); +} + +//---------------------------------------------------------------------------// +TEST_F(MuonicAtomSelectorTest, muonic_atom_asymmetric_mixture) +{ + // 70/30 mixture + real_type d_frac = 0.7; + real_type t_frac = 0.3; + MuonicAtomSelector select_atom(d_frac, t_frac); + + size_type num_samples = 10000; + size_type deuterium_count = 0; + + for (auto i : range(num_samples)) + { + auto atom = select_atom(this->rng()); + if (atom == MucfMuonicAtom::deuterium) + { + deuterium_count++; + } + } + + real_type expected_d_prob = this->calc_deuterium_q1s_prob(d_frac, t_frac); + real_type expected_d_count = num_samples * expected_d_prob; + // 3 sigma tolerance + real_type tolerance = 3 * this->calc_sigma(num_samples, expected_d_prob); + + EXPECT_NEAR(expected_d_count, deuterium_count, tolerance); +} + +//---------------------------------------------------------------------------// +TEST_F(MuonicAtomSelectorTest, spin_selector_deuterium) +{ + MuonicAtomSpinSelector select_spin(MucfMuonicAtom::deuterium); + + size_type num_samples = 10000; + size_type spin_3_2_count = 0; // Spin 3/2 + size_type spin_1_2_count = 0; // Spin 1/2 + + for (auto i : range(num_samples)) + { + size_type spin = select_spin(this->rng()); + if (spin == 3) + { + spin_3_2_count++; + } + else if (spin == 1) + { + spin_1_2_count++; + } + else + { + FAIL() << "Unexpected spin value: " << spin; + } + } + + EXPECT_EQ(num_samples, spin_3_2_count + spin_1_2_count); + + // Expected: 2/3 for spin 3/2, 1/3 for spin 1/2 + real_type expected_3_2_prob = 2.0 / 3.0; + real_type expected_3_2_count = num_samples * expected_3_2_prob; + // 3 sigma tolerance + real_type tolerance = 3 * this->calc_sigma(num_samples, expected_3_2_prob); + + EXPECT_NEAR(expected_3_2_count, spin_3_2_count, tolerance); + EXPECT_NEAR(num_samples - expected_3_2_count, spin_1_2_count, tolerance); +} + +//---------------------------------------------------------------------------// +TEST_F(MuonicAtomSelectorTest, spin_selector_tritium) +{ + MuonicAtomSpinSelector select_spin(MucfMuonicAtom::tritium); + + size_type num_samples = 10000; + size_type spin_1_count = 0; // Spin 1 + size_type spin_0_count = 0; // Spin 0 + + for (auto i : range(num_samples)) + { + size_type spin = select_spin(this->rng()); + if (spin == 2) + { + spin_1_count++; + } + else if (spin == 0) + { + spin_0_count++; + } + else + { + FAIL() << "Unexpected spin value: " << spin; + } + } + + EXPECT_EQ(num_samples, spin_1_count + spin_0_count); + + // Expected: 3/4 for spin 1, 1/4 for spin 0 + real_type expected_1_prob = 0.75; + real_type expected_1_count = num_samples * expected_1_prob; + // 3 sigma tolerance + real_type tolerance = 3 * this->calc_sigma(num_samples, expected_1_prob); + + EXPECT_NEAR(expected_1_count, spin_1_count, tolerance); + EXPECT_NEAR(num_samples - expected_1_count, spin_0_count, tolerance); +} + +//---------------------------------------------------------------------------// +} // namespace test +} // namespace detail +} // namespace celeritas From 7233b656215e7fa63afe9d0ce063e19a98c26cf8 Mon Sep 17 00:00:00 2001 From: Stefano Tognini Date: Mon, 26 Jan 2026 15:08:42 -0500 Subject: [PATCH 12/42] Update test harness accordingly --- test/celeritas/mucf/DTMucfInteractor.test.cc | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/test/celeritas/mucf/DTMucfInteractor.test.cc b/test/celeritas/mucf/DTMucfInteractor.test.cc index 3947878701..02a9e9787a 100644 --- a/test/celeritas/mucf/DTMucfInteractor.test.cc +++ b/test/celeritas/mucf/DTMucfInteractor.test.cc @@ -93,13 +93,16 @@ class DTMucfInteractorTest : public MucfInteractorHostTestBase EXPECT_EQ(host_data.particle_ids.mu_minus, sec[1].particle_id); EXPECT_EQ(host_data.particle_ids.alpha, sec[2].particle_id); - // Check energy conservation (17.6 MeV total) + // Check approximate energy conservation + // The total kineticy energy is only very roughly 17.6 MeV due to + // simplistic sampling. See DTMucfInteractor documentation for + // details. real_type total_kinetic_energy = 0; for (auto const& sec : interaction.secondaries) { total_kinetic_energy += sec.energy.value(); } - EXPECT_SOFT_EQ(17.6, total_kinetic_energy); + EXPECT_SOFT_NEAR(17.6, total_kinetic_energy, 0.5); // Check momentum conservation // Momentum and energy conservation is not accurate (see the From 490d9f0f9abacb07f10cbb84c74be1f523177417 Mon Sep 17 00:00:00 2001 From: Stefano Tognini Date: Tue, 27 Jan 2026 11:07:51 -0500 Subject: [PATCH 13/42] Move documentation --- .../mucf/interactor/DTMucfInteractor.hh | 23 ++++++++----------- test/celeritas/mucf/DTMucfInteractor.test.cc | 6 ++--- 2 files changed, 13 insertions(+), 16 deletions(-) diff --git a/src/celeritas/mucf/interactor/DTMucfInteractor.hh b/src/celeritas/mucf/interactor/DTMucfInteractor.hh index dc2bc71656..90ae732769 100644 --- a/src/celeritas/mucf/interactor/DTMucfInteractor.hh +++ b/src/celeritas/mucf/interactor/DTMucfInteractor.hh @@ -25,6 +25,14 @@ namespace celeritas * Fusion channels: * - \f$ \alpha + \mu + n \f$ * - \f$ (\alpha)_\mu + n \f$ + * + * \warning This implementation has an incorrect energy and momentum + * conservation implementation. Acceleron assumes an isotropic direction for + * both neutron and muon in the \f$ \alpha + \mu + n \f$ channel, which leads + * to the alpha particle either conserving energy or momentum but not both + * simultaneously. The current implementation results in a roughly correct + * total energy within \f$ K_\text{total} = [17.5, 17.9] \f$ MeV, instead of + * the expected 17.6 MeV. */ class DTMucfInteractor { @@ -94,17 +102,6 @@ DTMucfInteractor::DTMucfInteractor(NativeCRef const& data, //---------------------------------------------------------------------------// /*! * Sample a dt muonic molecule fusion. - * - * \note Since secondaries come from an at rest interaction, their final state - * is a simple combination of random direction + momentum conservation. - * - * \warning This implementation has an incorrect energy and momentum - * conservation implementation. Acceleron assumes an isotropic direction for - * both neutron and muon in the \f$ \alpha + \mu + n \f$ channel, which leads - * to the alpha particle either conserving energy or momentum but not both - * simultaneously. The current implementation results in a roughly correct - * total energy within \f$ K_\text{total} = [17.5, 17.9] \f$ MeV, instead of - * the expected 17.6 MeV. */ template CELER_FUNCTION Interaction DTMucfInteractor::operator()(Engine& rng) @@ -188,7 +185,7 @@ CELER_FUNCTION Interaction DTMucfInteractor::operator()(Engine& rng) //---------------------------------------------------------------------------// /*! * Calculate momentum magnitude from particle energy and mass via - * \f$ p = \sqrt{K^2 + 2mK} \f$ + * \f$ p = \sqrt{K^2 + 2mK} \f$ . */ CELER_FUNCTION real_type DTMucfInteractor::calc_momentum( units::MevEnergy energy, units::MevMass mass) const @@ -202,7 +199,7 @@ CELER_FUNCTION real_type DTMucfInteractor::calc_momentum( //---------------------------------------------------------------------------// /*! * Calculate kinetic energy given a momentum magnitude and particle mass via - * \f$ K = \sqrt{p^2 - m^2} - m\f$ + * \f$ K = \sqrt{p^2 - m^2} - m\f$ . */ CELER_FUNCTION units::MevEnergy DTMucfInteractor::calc_kinetic_energy(Real3 momentum_vec, diff --git a/test/celeritas/mucf/DTMucfInteractor.test.cc b/test/celeritas/mucf/DTMucfInteractor.test.cc index 02a9e9787a..c7cd9fe3a7 100644 --- a/test/celeritas/mucf/DTMucfInteractor.test.cc +++ b/test/celeritas/mucf/DTMucfInteractor.test.cc @@ -187,7 +187,7 @@ TEST_F(DTMucfInteractorTest, alpha_muon_neutron) auto& rng = this->rng(); for (auto i : range(num_samples)) { - Interaction result = interact(rng); + auto result = interact(rng); this->validate_interaction(result, channel); } } @@ -208,7 +208,7 @@ TEST_F(DTMucfInteractorTest, muonicalpha_neutron) auto& rng = this->rng(); for (auto i : range(num_samples)) { - Interaction result = interact(rng); + auto result = interact(rng); this->validate_interaction(result, channel); } } @@ -230,7 +230,7 @@ TEST_F(DTMucfInteractorTest, stress_test) auto& rng = this->rng(); for (auto i : range(num_samples)) { - Interaction result = interact(rng); + auto result = interact(rng); total_avg_secondaries += result.secondaries.size(); } } From deb706db6951838e7ebad11993e92855b42c62a1 Mon Sep 17 00:00:00 2001 From: Stefano Tognini Date: Tue, 27 Jan 2026 14:00:18 -0500 Subject: [PATCH 14/42] Add `[[maybe_unused]]` to iterator counters in the test. --- test/celeritas/mucf/DTMucfInteractor.test.cc | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/celeritas/mucf/DTMucfInteractor.test.cc b/test/celeritas/mucf/DTMucfInteractor.test.cc index c7cd9fe3a7..42bbcf6008 100644 --- a/test/celeritas/mucf/DTMucfInteractor.test.cc +++ b/test/celeritas/mucf/DTMucfInteractor.test.cc @@ -185,7 +185,7 @@ TEST_F(DTMucfInteractorTest, alpha_muon_neutron) data_.host_ref(), channel, this->secondary_allocator()); auto& rng = this->rng(); - for (auto i : range(num_samples)) + for ([[maybe_unused]] auto i : range(num_samples)) { auto result = interact(rng); this->validate_interaction(result, channel); @@ -206,7 +206,7 @@ TEST_F(DTMucfInteractorTest, muonicalpha_neutron) data_.host_ref(), channel, this->secondary_allocator()); auto& rng = this->rng(); - for (auto i : range(num_samples)) + for ([[maybe_unused]] auto i : range(num_samples)) { auto result = interact(rng); this->validate_interaction(result, channel); @@ -228,7 +228,7 @@ TEST_F(DTMucfInteractorTest, stress_test) data_.host_ref(), channel, this->secondary_allocator()); auto& rng = this->rng(); - for (auto i : range(num_samples)) + for ([[maybe_unused]] auto i : range(num_samples)) { auto result = interact(rng); total_avg_secondaries += result.secondaries.size(); From 14d912539f63f0062cf00647b6ada1a252e4eb4e Mon Sep 17 00:00:00 2001 From: Stefano Tognini Date: Wed, 28 Jan 2026 10:43:00 -0500 Subject: [PATCH 15/42] Address comments --- src/celeritas/inp/MucfPhysics.cc | 53 +++++++++++-------- src/celeritas/mucf/data/DTMixMucfData.hh | 3 +- .../mucf/interactor/DTMucfInteractor.hh | 23 ++++---- .../mucf/MucfInteractorHostTestBase.cc | 29 ++++++---- 4 files changed, 63 insertions(+), 45 deletions(-) diff --git a/src/celeritas/inp/MucfPhysics.cc b/src/celeritas/inp/MucfPhysics.cc index fcc9047182..6a51bf5f5e 100644 --- a/src/celeritas/inp/MucfPhysics.cc +++ b/src/celeritas/inp/MucfPhysics.cc @@ -15,6 +15,8 @@ namespace //---------------------------------------------------------------------------// /*! * Muon energy CDF data for muon-catalyzed fusion. + * + * Data is extracted from https://doi.org/10.1103/PhysRevC.107.034607 . */ Grid mucf_muon_energy_cdf() { @@ -46,28 +48,30 @@ Grid mucf_muon_energy_cdf() 1.0, }; - // Energy [MeV] - cdf.y = {0., - 0.00048850540675768084, - 0.0008390389347819425, - 0.0012521213482687141, - 0.0017153033196164724, - 0.002253638712180777, - 0.002854653691809707, - 0.003606073540073316, - 0.004470346052913727, - 0.005560291219507215, - 0.006700556502915258, - 0.007953772477101693, - 0.009194596305637525, - 0.010849180562221111, - 0.012353474314071864, - 0.014045888515617822, - 0.015650634617544647, - 0.01738079707555165, - 0.019111008546659452, - 0.019976130619913615, - 0.08}; + // Muon energy [MeV] + cdf.y = { + 0., + 0.00048850540675768084, + 0.0008390389347819425, + 0.0012521213482687141, + 0.0017153033196164724, + 0.002253638712180777, + 0.002854653691809707, + 0.003606073540073316, + 0.004470346052913727, + 0.005560291219507215, + 0.006700556502915258, + 0.007953772477101693, + 0.009194596305637525, + 0.010849180562221111, + 0.012353474314071864, + 0.014045888515617822, + 0.015650634617544647, + 0.01738079707555165, + 0.019111008546659452, + 0.019976130619913615, + 0.08, + }; CELER_ENSURE(cdf); return cdf; @@ -108,6 +112,9 @@ std::vector mucf_atom_transfer_rates() //---------------------------------------------------------------------------// /*! * Initialize with hardcoded values. + * + * \note Default data replicated from Acceleron. + * \todo May need to be updated. */ MucfScalars MucfScalars::from_default() { @@ -123,7 +130,7 @@ MucfScalars MucfScalars::from_default() /*! * Construct hardcoded muon-catalyzed fusion physics data. * - * \todo Official tables are implemented directly in C++ in NK Labs Geant4 + * \todo Official tables are implemented directly in C++ in Acceleron's Geant4 * implementation. These likely will be loaded from external data files so that * muCF experimental data can be loaded at runtime. */ diff --git a/src/celeritas/mucf/data/DTMixMucfData.hh b/src/celeritas/mucf/data/DTMixMucfData.hh index 9f0ffbc66a..7d296b1fd8 100644 --- a/src/celeritas/mucf/data/DTMixMucfData.hh +++ b/src/celeritas/mucf/data/DTMixMucfData.hh @@ -110,7 +110,8 @@ struct DTMixMucfData MucfParticleMasses particle_masses; //! Muon CDF energy grid for sampling outgoing muCF muons - GridRecord muon_energy_cdf; //! \todo Verify energy unit + //! X-axis range is [0, 1) and y-axis is the outgoing muon energy in MeV + GridRecord muon_energy_cdf; Items reals; //!@{ diff --git a/src/celeritas/mucf/interactor/DTMucfInteractor.hh b/src/celeritas/mucf/interactor/DTMucfInteractor.hh index 90ae732769..fb7b131a91 100644 --- a/src/celeritas/mucf/interactor/DTMucfInteractor.hh +++ b/src/celeritas/mucf/interactor/DTMucfInteractor.hh @@ -75,12 +75,13 @@ class DTMucfInteractor } // Calculate momentum magnitude from energy and mass - inline CELER_FUNCTION real_type calc_momentum(units::MevEnergy energy, - units::MevMass mass) const; + inline CELER_FUNCTION real_type calc_momentum( + units::MevEnergy const energy, units::MevMass const mass) const; - // Calculate kinetic energy from momentum magnitude and mass + // Calculate kinetic energy from momentum and mass inline CELER_FUNCTION units::MevEnergy - calc_kinetic_energy(Real3 momentum_vec, units::MevMass mass) const; + calc_kinetic_energy(Real3 const& momentum_vec, + units::MevMass const mass) const; }; //---------------------------------------------------------------------------// @@ -119,6 +120,7 @@ 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); @@ -188,7 +190,7 @@ CELER_FUNCTION Interaction DTMucfInteractor::operator()(Engine& rng) * \f$ p = \sqrt{K^2 + 2mK} \f$ . */ CELER_FUNCTION real_type DTMucfInteractor::calc_momentum( - units::MevEnergy energy, units::MevMass mass) const + units::MevEnergy const energy, units::MevMass const mass) const { return std::sqrt(ipow<2>(value_as(energy)) @@ -198,17 +200,14 @@ CELER_FUNCTION real_type DTMucfInteractor::calc_momentum( //---------------------------------------------------------------------------// /*! - * Calculate kinetic energy given a momentum magnitude and particle mass via + * 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 momentum_vec, - units::MevMass mass) const +DTMucfInteractor::calc_kinetic_energy(Real3 const& momentum_vec, + units::MevMass const mass) const { - real_type momentum_mag = std::sqrt(ipow<2>(momentum_vec[0]) - + ipow<2>(momentum_vec[1]) - + ipow<2>(momentum_vec[2])); - + real_type momentum_mag = norm(momentum_vec); return units::MevEnergy{std::sqrt(ipow<2>(momentum_mag) + ipow<2>(value_as(mass))) - value_as(mass)}; diff --git a/test/celeritas/mucf/MucfInteractorHostTestBase.cc b/test/celeritas/mucf/MucfInteractorHostTestBase.cc index d0139b23bb..e39de6f2bb 100644 --- a/test/celeritas/mucf/MucfInteractorHostTestBase.cc +++ b/test/celeritas/mucf/MucfInteractorHostTestBase.cc @@ -7,6 +7,7 @@ #include "MucfInteractorHostTestBase.hh" #include "celeritas/Units.hh" +#include "celeritas/inp/MucfPhysics.hh" namespace celeritas { @@ -26,17 +27,27 @@ MucfInteractorHostBase::MucfInteractorHostBase() using units::Second; using InvSecond = RealQuantity>; - // Particle masses (from PDG) - constexpr MevMass muon_mass{105.6583745}; - constexpr MevMass protium_mass{938.272081}; - constexpr MevMass neutron_mass{939.565413}; - constexpr MevMass deuterium_mass{1875.612928}; - constexpr MevMass tritium_mass{2808.921132}; - constexpr MevMass alpha_mass{3727.379378}; - constexpr MevMass he3_mass{2808.391607}; + constexpr units::MevMass amu_mev{931.5}; // Convert from AMU to MeV + auto const scalars = inp::MucfScalars::from_default(); + + // Particle masses + // PDG, PRD 110, 030001, 2024 (https://doi.org/10.1103/PhysRevD.110.030001) + constexpr MevMass muon_mass{105.6583755}; + constexpr MevMass protium_mass{938.272088}; + constexpr MevMass neutron_mass{939.565420}; + // Acceleron default values + MevMass deuterium_mass{scalars.deuterium.value() * amu_mev}; + MevMass tritium_mass{scalars.tritium.value() * amu_mev}; + // CODATA 2022 (https://arxiv.org/pdf/2409.03787) + MevMass alpha_mass{3727.379}; + MevMass he3_mass{2808.391}; + // Decay constants + // Muon: PDG, 110, 030001, 2024 + // Tritium: NUBASE 2020, Chinese Physics C 45 030001 + // (https://iopscience.iop.org/article/10.1088/1674-1137/abddae) constexpr InvSecond muon_decay_constant{1 / 2.1969811e-6}; - constexpr InvSecond tritium_decay_constant{1 / 5.4927e8}; // ~12.3 years + constexpr InvSecond tritium_decay_constant{1 / 3.8879e+8}; // ParticlesParams used by the muCF process ParticleParams::Input par_inp = { From 3a8e3e3320a639fa2f8f96968e99631b830e3ca2 Mon Sep 17 00:00:00 2001 From: Stefano Tognini Date: Wed, 28 Jan 2026 13:41:47 -0500 Subject: [PATCH 16/42] Add [[maybe_unused]] to counters in the test --- .../celeritas/mucf/MuonicAtomSelector.test.cc | 59 ++++++++++--------- 1 file changed, 32 insertions(+), 27 deletions(-) diff --git a/test/celeritas/mucf/MuonicAtomSelector.test.cc b/test/celeritas/mucf/MuonicAtomSelector.test.cc index 76fa6c354c..7d1d39911e 100644 --- a/test/celeritas/mucf/MuonicAtomSelector.test.cc +++ b/test/celeritas/mucf/MuonicAtomSelector.test.cc @@ -56,12 +56,12 @@ class MuonicAtomSelectorTest : public Test TEST_F(MuonicAtomSelectorTest, muonic_atom_pure_deuterium) { // Pure deuterium - MuonicAtomSelector select_atom(1.0, 0.0); + MuonicAtomSelector select_atom(1, 0); - size_type num_samples = 1000; + size_type const num_samples = 100; size_type deuterium_count = 0; - for (auto i : range(num_samples)) + for ([[maybe_unused]] auto i : range(num_samples)) { auto atom = select_atom(this->rng()); if (atom == MucfMuonicAtom::deuterium) @@ -76,12 +76,12 @@ TEST_F(MuonicAtomSelectorTest, muonic_atom_pure_deuterium) TEST_F(MuonicAtomSelectorTest, muonic_atom_pure_tritium) { // Pure tritium - MuonicAtomSelector select_atom(0.0, 1.0); + MuonicAtomSelector select_atom(0, 1); - size_type num_samples = 1000; + size_type const num_samples = 100; size_type tritium_count = 0; - for (auto i : range(num_samples)) + for ([[maybe_unused]] auto i : range(num_samples)) { auto atom = select_atom(this->rng()); if (atom == MucfMuonicAtom::tritium) @@ -96,15 +96,15 @@ TEST_F(MuonicAtomSelectorTest, muonic_atom_pure_tritium) TEST_F(MuonicAtomSelectorTest, muonic_atom_dt_mixture) { // 50/50 mixture - real_type d_frac = 0.5; - real_type t_frac = 0.5; + real_type const d_frac = 0.5; + real_type const t_frac = 0.5; MuonicAtomSelector select_atom(d_frac, t_frac); - size_type num_samples = 10000; + size_type const num_samples = 10000; size_type deuterium_count = 0; size_type tritium_count = 0; - for (auto i : range(num_samples)) + for ([[maybe_unused]] auto i : range(num_samples)) { auto atom = select_atom(this->rng()); if (atom == MucfMuonicAtom::deuterium) @@ -119,10 +119,12 @@ TEST_F(MuonicAtomSelectorTest, muonic_atom_dt_mixture) EXPECT_EQ(num_samples, deuterium_count + tritium_count); - real_type expected_d_prob = this->calc_deuterium_q1s_prob(d_frac, t_frac); - real_type expected_d_count = num_samples * expected_d_prob; + real_type const expected_d_prob + = this->calc_deuterium_q1s_prob(d_frac, t_frac); + real_type const expected_d_count = num_samples * expected_d_prob; // 3 sigma tolerance - real_type tolerance = 3 * this->calc_sigma(num_samples, expected_d_prob); + real_type const tolerance + = 3 * this->calc_sigma(num_samples, expected_d_prob); EXPECT_NEAR(expected_d_count, deuterium_count, tolerance); EXPECT_NEAR(num_samples - expected_d_count, tritium_count, tolerance); @@ -132,14 +134,14 @@ TEST_F(MuonicAtomSelectorTest, muonic_atom_dt_mixture) TEST_F(MuonicAtomSelectorTest, muonic_atom_asymmetric_mixture) { // 70/30 mixture - real_type d_frac = 0.7; - real_type t_frac = 0.3; + real_type const d_frac = 0.7; + real_type const t_frac = 0.3; MuonicAtomSelector select_atom(d_frac, t_frac); - size_type num_samples = 10000; + size_type const num_samples = 10000; size_type deuterium_count = 0; - for (auto i : range(num_samples)) + for ([[maybe_unused]] auto i : range(num_samples)) { auto atom = select_atom(this->rng()); if (atom == MucfMuonicAtom::deuterium) @@ -148,10 +150,12 @@ TEST_F(MuonicAtomSelectorTest, muonic_atom_asymmetric_mixture) } } - real_type expected_d_prob = this->calc_deuterium_q1s_prob(d_frac, t_frac); - real_type expected_d_count = num_samples * expected_d_prob; + real_type const expected_d_prob + = this->calc_deuterium_q1s_prob(d_frac, t_frac); + real_type const expected_d_count = num_samples * expected_d_prob; // 3 sigma tolerance - real_type tolerance = 3 * this->calc_sigma(num_samples, expected_d_prob); + real_type const tolerance + = 3 * this->calc_sigma(num_samples, expected_d_prob); EXPECT_NEAR(expected_d_count, deuterium_count, tolerance); } @@ -161,11 +165,11 @@ TEST_F(MuonicAtomSelectorTest, spin_selector_deuterium) { MuonicAtomSpinSelector select_spin(MucfMuonicAtom::deuterium); - size_type num_samples = 10000; + size_type const num_samples = 10000; size_type spin_3_2_count = 0; // Spin 3/2 size_type spin_1_2_count = 0; // Spin 1/2 - for (auto i : range(num_samples)) + for ([[maybe_unused]] auto i : range(num_samples)) { size_type spin = select_spin(this->rng()); if (spin == 3) @@ -199,11 +203,11 @@ TEST_F(MuonicAtomSelectorTest, spin_selector_tritium) { MuonicAtomSpinSelector select_spin(MucfMuonicAtom::tritium); - size_type num_samples = 10000; + size_type const num_samples = 10000; size_type spin_1_count = 0; // Spin 1 size_type spin_0_count = 0; // Spin 0 - for (auto i : range(num_samples)) + for ([[maybe_unused]] auto i : range(num_samples)) { size_type spin = select_spin(this->rng()); if (spin == 2) @@ -223,10 +227,11 @@ TEST_F(MuonicAtomSelectorTest, spin_selector_tritium) EXPECT_EQ(num_samples, spin_1_count + spin_0_count); // Expected: 3/4 for spin 1, 1/4 for spin 0 - real_type expected_1_prob = 0.75; - real_type expected_1_count = num_samples * expected_1_prob; + real_type const expected_1_prob = 0.75; + real_type const expected_1_count = num_samples * expected_1_prob; // 3 sigma tolerance - real_type tolerance = 3 * this->calc_sigma(num_samples, expected_1_prob); + real_type const tolerance + = 3 * this->calc_sigma(num_samples, expected_1_prob); EXPECT_NEAR(expected_1_count, spin_1_count, tolerance); EXPECT_NEAR(num_samples - expected_1_count, spin_0_count, tolerance); From 8775c9f29888b41d4c6d6e6f543b662782d1ad23 Mon Sep 17 00:00:00 2001 From: Stefano Tognini Date: Wed, 28 Jan 2026 14:16:17 -0500 Subject: [PATCH 17/42] Improve documentation --- doc/CMakeLists.txt | 1 + doc/implementation/mucf-physics.rst | 7 +++++-- src/celeritas/mucf/executor/detail/MuonicAtomSelector.hh | 4 ++-- src/celeritas/mucf/executor/detail/TTChannelSelector.hh | 2 +- 4 files changed, 9 insertions(+), 5 deletions(-) diff --git a/doc/CMakeLists.txt b/doc/CMakeLists.txt index d797729cce..e2ea4d5458 100644 --- a/doc/CMakeLists.txt +++ b/doc/CMakeLists.txt @@ -270,6 +270,7 @@ file(GLOB _DOXYGEN_SOURCE "${PROJECT_SOURCE_DIR}/src/celeritas/em/xs/*.hh" "${PROJECT_SOURCE_DIR}/src/celeritas/mucf/data/*.hh" "${PROJECT_SOURCE_DIR}/src/celeritas/mucf/executor/*.hh" + "${PROJECT_SOURCE_DIR}/src/celeritas/mucf/executor/detail/*.hh" "${PROJECT_SOURCE_DIR}/src/celeritas/mucf/interactor/*.hh" "${PROJECT_SOURCE_DIR}/src/celeritas/mucf/model/*.hh" "${PROJECT_SOURCE_DIR}/src/celeritas/mucf/model/detail/*.hh" diff --git a/doc/implementation/mucf-physics.rst b/doc/implementation/mucf-physics.rst index 75e13de3bf..ff014d715d 100644 --- a/doc/implementation/mucf-physics.rst +++ b/doc/implementation/mucf-physics.rst @@ -133,11 +133,14 @@ host/device data via :cpp:class:`celeritas::detail::MucfMaterialInserter`. .. doxygenclass:: celeritas::detail::MucfMaterialInserter The main cycle is managed by the model's -:cpp:class:`celeritas::DTMixMucfExecutor`, with the Interactors reserved for -sampling final states of the outgoing secondaries. +:cpp:class:`celeritas::DTMixMucfExecutor`. The muonic atom selection is handled +by detail classes, while the Interactors are reserved for sampling final states +of the outgoing secondaries. .. note:: Only reactive channels are implemented. +.. doxygenclass:: celeritas::detail::MuonicAtomSelector +.. doxygenclass:: celeritas::detail::MuonicAtomSpinSelector .. doxygenclass:: celeritas::DDMucfInteractor .. doxygenclass:: celeritas::DTMucfInteractor .. doxygenclass:: celeritas::TTMucfInteractor diff --git a/src/celeritas/mucf/executor/detail/MuonicAtomSelector.hh b/src/celeritas/mucf/executor/detail/MuonicAtomSelector.hh index 81b674d38d..8eece1f035 100644 --- a/src/celeritas/mucf/executor/detail/MuonicAtomSelector.hh +++ b/src/celeritas/mucf/executor/detail/MuonicAtomSelector.hh @@ -27,9 +27,9 @@ namespace detail * where \f$ C_t \f$ is the relative tritium isotope concentration and * \f$ q_\text{1s} \f$ is the fraction muonic deuterium atoms in the ground * state. This expression allows calculating the probability of forming a - * muonic deuterium or tritium atom via + * muonic deuterium atom via * \f[ - * P_\text{d} = C_d * q\text{1s}. + * P_\text{d} = C_d \times q\text{1s}. * \f] * * If a selected uniform random number is \f$ x \leq P_\text{d} \f$, a muonic diff --git a/src/celeritas/mucf/executor/detail/TTChannelSelector.hh b/src/celeritas/mucf/executor/detail/TTChannelSelector.hh index 1f7c0b40ea..dd8372479a 100644 --- a/src/celeritas/mucf/executor/detail/TTChannelSelector.hh +++ b/src/celeritas/mucf/executor/detail/TTChannelSelector.hh @@ -2,7 +2,7 @@ // Copyright Celeritas contributors: see top-level COPYRIGHT file for details // SPDX-License-Identifier: (Apache-2.0 OR MIT) //---------------------------------------------------------------------------// -//! \file celeritas/mucf/executor/detail/TTChannelSelectionHelper.hh +//! \file celeritas/mucf/executor/detail/TTChannelSelector.hh //---------------------------------------------------------------------------// #pragma once From cc7703a2a5d35fc2f6ab9bca02e53fc56b8acb96 Mon Sep 17 00:00:00 2001 From: Stefano Tognini Date: Wed, 28 Jan 2026 14:45:57 -0500 Subject: [PATCH 18/42] Draft dt channel selector --- .../mucf/executor/detail/DTChannelSelector.hh | 40 +++++++++---------- 1 file changed, 20 insertions(+), 20 deletions(-) diff --git a/src/celeritas/mucf/executor/detail/DTChannelSelector.hh b/src/celeritas/mucf/executor/detail/DTChannelSelector.hh index 47e097cdfd..7eb75b9191 100644 --- a/src/celeritas/mucf/executor/detail/DTChannelSelector.hh +++ b/src/celeritas/mucf/executor/detail/DTChannelSelector.hh @@ -6,6 +6,7 @@ //---------------------------------------------------------------------------// #pragma once +#include "corecel/random/distribution/GenerateCanonical.hh" #include "celeritas/mucf/interactor/DTMucfInteractor.hh" namespace celeritas @@ -16,8 +17,9 @@ namespace detail /*! * Select final channel for muonic dt molecules. * - * This selection already accounts for sticking, as that is one of the possible - * outcomes. + * The selection is a simple selection based on a constant sticking fraction + * from [ \todo https://arxiv.org/abs/2112.08399 ], in which ~0.8% of the time + * the muonic alpha channel is selected. */ class DTChannelSelector { @@ -28,40 +30,38 @@ class DTChannelSelector //!@} public: - //! Construct with args; \todo Update documentation - inline CELER_FUNCTION DTChannelSelector(/* args */); + //! Default constructor + inline CELER_FUNCTION DTChannelSelector() = default; // Select fusion channel to be used by the interactor template - inline CELER_FUNCTION Channel operator()(Engine&); + inline CELER_FUNCTION Channel operator()(Engine& rng); + + private: + // Constant sticking fraction of dt fusion + inline CELER_FUNCTION real_type sticking_fraction() const + { + return 0.00857; + } }; //---------------------------------------------------------------------------// // INLINE DEFINITIONS //---------------------------------------------------------------------------// /*! - * Construct with args. - * - * \todo Update documentation - */ -CELER_FUNCTION DTChannelSelector::DTChannelSelector(/* args */) -{ - CELER_NOT_IMPLEMENTED("Mucf dt fusion channel selection"); -} - -//---------------------------------------------------------------------------// -/*! - * Return a sampled channel to be used as input in the dt muCF interactor. + * Return a selected fusion channel for the \f$ (dt)_\mu \f$ muonic molecule. * * \sa celeritas::DTMucfInteractor */ template -CELER_FUNCTION DTChannelSelector::Channel DTChannelSelector::operator()(Engine&) +CELER_FUNCTION DTChannelSelector::Channel +DTChannelSelector::operator()(Engine& rng) { Channel result{Channel::size_}; - //! \todo Implement - // Final channel selection already takes into account sticking. + generate_canonical(rng) > this->sticking_fraction() + ? result = Channel::alpha_muon_neutron + : result = Channel::muonicalpha_neutron; CELER_ENSURE(result < Channel::size_); return result; From d5c9295413b9dcfa240f35d46b1dd5142d3893fd Mon Sep 17 00:00:00 2001 From: Stefano Tognini Date: Wed, 28 Jan 2026 15:07:54 -0500 Subject: [PATCH 19/42] Draft dd channel selector --- .../mucf/executor/DTMixMucfExecutor.hh | 7 +- .../mucf/executor/detail/DDChannelSelector.hh | 71 +++++++++++++++---- .../mucf/executor/detail/DTChannelSelector.hh | 1 - .../mucf/interactor/DDMucfInteractor.hh | 4 +- 4 files changed, 64 insertions(+), 19 deletions(-) diff --git a/src/celeritas/mucf/executor/DTMixMucfExecutor.hh b/src/celeritas/mucf/executor/DTMixMucfExecutor.hh index 64d7f4b82d..3de120a9b1 100644 --- a/src/celeritas/mucf/executor/DTMixMucfExecutor.hh +++ b/src/celeritas/mucf/executor/DTMixMucfExecutor.hh @@ -44,7 +44,8 @@ DTMixMucfExecutor::operator()(celeritas::CoreTrackView const& track) auto elcomp_id = phys_step_view.element(); CELER_ASSERT(elcomp_id); - auto element = track.material().material_record().element_record(elcomp_id); + auto const& mat_record = track.material().material_record(); + auto element = mat_record.element_record(elcomp_id); CELER_ASSERT(element.atomic_number() == AtomicNumber{1}); // Must be H auto rng = track.rng(); @@ -118,7 +119,9 @@ DTMixMucfExecutor::operator()(celeritas::CoreTrackView const& track) case MucfMuonicMolecule::deuterium_deuterium: { // Return DD interaction DDMucfInteractor interact( - data, detail::DDChannelSelector()(rng), allocate_secondaries); + data, + detail::DDChannelSelector(mat_record.temperature())(rng), + allocate_secondaries); result = interact(rng); break; } diff --git a/src/celeritas/mucf/executor/detail/DDChannelSelector.hh b/src/celeritas/mucf/executor/detail/DDChannelSelector.hh index 05394c8019..a6c645b988 100644 --- a/src/celeritas/mucf/executor/detail/DDChannelSelector.hh +++ b/src/celeritas/mucf/executor/detail/DDChannelSelector.hh @@ -6,6 +6,7 @@ //---------------------------------------------------------------------------// #pragma once +#include "corecel/random/distribution/GenerateCanonical.hh" #include "celeritas/mucf/interactor/DDMucfInteractor.hh" namespace celeritas @@ -16,8 +17,9 @@ namespace detail /*! * Select final channel for muonic dd molecules. * - * This selection already accounts for sticking, as that is one of the possible - * outcomes. + * The branching ratio is temperature dependent and determines the probability + * of the outcome of the fusion ending in the \f$ ^3\text{He} \f$ channels + * versus the tritium channel. */ class DDChannelSelector { @@ -27,41 +29,82 @@ class DDChannelSelector using Channel = DDMucfInteractor::Channel; //!@} - public: - //! Construct with args; \todo Update documentation - inline CELER_FUNCTION DDChannelSelector(/* args */); + //! Construct with material temperature + inline CELER_FUNCTION DDChannelSelector(real_type material_temperature); // Select fusion channel to be used by the interactor template - inline CELER_FUNCTION Channel operator()(Engine&); + inline CELER_FUNCTION Channel operator()(Engine& rng); + + private: + real_type he3_probability_{}; + + // Constant sticking fraction between the two 3He channels + inline CELER_FUNCTION real_type sticking_fraction() const { return 0.122; } }; //---------------------------------------------------------------------------// // INLINE DEFINITIONS //---------------------------------------------------------------------------// /*! - * Construct with args. + * Construct with material temperature. * - * \todo Update documentation + * The temperature is used to calculate the branching ratio between either of + * the \f$ ^3\text{He} \f$ channels and the tritium channel. */ -CELER_FUNCTION DDChannelSelector::DDChannelSelector(/* args */) +CELER_FUNCTION +DDChannelSelector::DDChannelSelector(real_type material_temperature) { - CELER_NOT_IMPLEMENTED("Mucf dd fusion channel selection"); + CELER_EXPECT(material_temperature > 0); + + real_type branching_ratio{0}; + if (material_temperature < 50) + { + branching_ratio = 1; + } + else if (material_temperature < 100) + { + branching_ratio = 1.0088 * (material_temperature - 50); + } + else + { + branching_ratio = 1.44; + } + CELER_ASSERT(branching_ratio > 0); + + he3_probability_ = branching_ratio / (branching_ratio + 1); + CELER_ENSURE(he3_probability_ > 0 && he3_probability_ <= 1); } //---------------------------------------------------------------------------// /*! - * Return a sampled channel to be used as input in the dd muCF interactor. + * Return a selected fusion channel for the \f$ (dd)_\mu \f$ muonic molecule. * * \sa celeritas::DDMucfInteractor */ template -CELER_FUNCTION DDChannelSelector::Channel DDChannelSelector::operator()(Engine&) +CELER_FUNCTION DDChannelSelector::Channel +DDChannelSelector::operator()(Engine& rng) { Channel result{Channel::size_}; - //! \todo Implement - // Final channel selection already takes into account sticking. + if (generate_canonical(rng) < he3_probability_) + { + // Select between the two 3He channels + if (generate_canonical(rng) > this->sticking_fraction()) + { + result = Channel::helium3_muon_neutron; + } + else + { + result = Channel::muonichelium3_neutron; + } + } + else + { + // Select tritium channel + result = Channel::tritium_muon_proton; + } CELER_ENSURE(result < Channel::size_); return result; diff --git a/src/celeritas/mucf/executor/detail/DTChannelSelector.hh b/src/celeritas/mucf/executor/detail/DTChannelSelector.hh index 7eb75b9191..4909c054a9 100644 --- a/src/celeritas/mucf/executor/detail/DTChannelSelector.hh +++ b/src/celeritas/mucf/executor/detail/DTChannelSelector.hh @@ -29,7 +29,6 @@ class DTChannelSelector using Channel = DTMucfInteractor::Channel; //!@} - public: //! Default constructor inline CELER_FUNCTION DTChannelSelector() = default; diff --git a/src/celeritas/mucf/interactor/DDMucfInteractor.hh b/src/celeritas/mucf/interactor/DDMucfInteractor.hh index b75fa71c02..cf29b349c7 100644 --- a/src/celeritas/mucf/interactor/DDMucfInteractor.hh +++ b/src/celeritas/mucf/interactor/DDMucfInteractor.hh @@ -29,12 +29,12 @@ namespace celeritas class DDMucfInteractor { public: - //! \todo Implement muonichydrogen3_proton (\f$ (^3\text{H})_\mu + p \f$) + //! \todo Implement muonictritium_proton (\f$ (^3\text{H})_\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$ ^3\text{H} + \mu + p \f$ size_ }; From 6e952d0e5d054ce541cd909c785061e5d3f61458 Mon Sep 17 00:00:00 2001 From: Stefano Tognini Date: Wed, 28 Jan 2026 15:15:37 -0500 Subject: [PATCH 20/42] Draft tt channel selector --- .../mucf/executor/detail/DTChannelSelector.hh | 6 +-- .../mucf/executor/detail/TTChannelSelector.hh | 44 +++++++++---------- 2 files changed, 25 insertions(+), 25 deletions(-) diff --git a/src/celeritas/mucf/executor/detail/DTChannelSelector.hh b/src/celeritas/mucf/executor/detail/DTChannelSelector.hh index 4909c054a9..76246281f4 100644 --- a/src/celeritas/mucf/executor/detail/DTChannelSelector.hh +++ b/src/celeritas/mucf/executor/detail/DTChannelSelector.hh @@ -17,9 +17,9 @@ namespace detail /*! * Select final channel for muonic dt molecules. * - * The selection is a simple selection based on a constant sticking fraction - * from [ \todo https://arxiv.org/abs/2112.08399 ], in which ~0.8% of the time - * the muonic alpha channel is selected. + * The selection is based on a constant sticking fraction from + * [ \todo https://arxiv.org/abs/2112.08399 ], in which ~0.8% of the time the + * muonic alpha channel is selected. */ class DTChannelSelector { diff --git a/src/celeritas/mucf/executor/detail/TTChannelSelector.hh b/src/celeritas/mucf/executor/detail/TTChannelSelector.hh index dd8372479a..3c797a0374 100644 --- a/src/celeritas/mucf/executor/detail/TTChannelSelector.hh +++ b/src/celeritas/mucf/executor/detail/TTChannelSelector.hh @@ -6,6 +6,7 @@ //---------------------------------------------------------------------------// #pragma once +#include "corecel/random/distribution/GenerateCanonical.hh" #include "celeritas/mucf/interactor/TTMucfInteractor.hh" namespace celeritas @@ -16,8 +17,9 @@ namespace detail /*! * Select final channel for muonic tt molecules. * - * This selection already accounts for sticking, as that is one of the possible - * outcomes. + * The selection is based on a constant sticking fraction + * from [ \todo https://link.springer.com/article/10.1134/S1063776109020034 ], + * in which ~14% of the time the muonic alpha channel is selected. */ class TTChannelSelector { @@ -27,42 +29,40 @@ class TTChannelSelector using Channel = TTMucfInteractor::Channel; //!@} - public: - //! Construct with args; \todo Update documentation - inline CELER_FUNCTION TTChannelSelector(/* args */); + //! Default constructor + inline CELER_FUNCTION TTChannelSelector() = default; // Select fusion channel to be used by the interactor template - inline CELER_FUNCTION Channel operator()(Engine&); + inline CELER_FUNCTION Channel operator()(Engine& rng); + + private: + // Constant sticking fraction of tt fusion + inline CELER_FUNCTION real_type sticking_fraction() const { return 0.14; } }; //---------------------------------------------------------------------------// // INLINE DEFINITIONS //---------------------------------------------------------------------------// /*! - * Construct with args. - * - * \todo Update documentation - */ -CELER_FUNCTION TTChannelSelector::TTChannelSelector(/* args */) -{ - //! \todo Implement - CELER_NOT_IMPLEMENTED("Mucf tt fusion channel selection"); -} - -//---------------------------------------------------------------------------// -/*! - * Return a sampled channel to be used as input in the tt muCF interactor. + * Return a selected fusion channel for the \f$ (tt)_\mu \f$ muonic molecule. * * \sa celeritas::TTMucfInteractor */ template -CELER_FUNCTION TTChannelSelector::Channel TTChannelSelector::operator()(Engine&) +CELER_FUNCTION TTChannelSelector::Channel +TTChannelSelector::operator()(Engine& rng) { Channel result{Channel::size_}; - //! \todo Implement - // Final channel selection already takes into account sticking. + if (generate_canonical(rng) > sticking_fraction()) + { + result = Channel::alpha_muon_neutron_neutron; + } + else + { + result = Channel::muonicalpha_neutron_neutron; + } CELER_ENSURE(result < Channel::size_); return result; From b616562ba9965c45000cb5eef2ea7803714d3ec8 Mon Sep 17 00:00:00 2001 From: Stefano Tognini Date: Wed, 28 Jan 2026 16:15:52 -0500 Subject: [PATCH 21/42] Draft test harness --- .../mucf/executor/detail/DDChannelSelector.hh | 11 +- test/celeritas/CMakeLists.txt | 1 + test/celeritas/mucf/ChannelSelector.test.cc | 314 ++++++++++++++++++ .../celeritas/mucf/MuonicAtomSelector.test.cc | 2 +- 4 files changed, 322 insertions(+), 6 deletions(-) create mode 100644 test/celeritas/mucf/ChannelSelector.test.cc diff --git a/src/celeritas/mucf/executor/detail/DDChannelSelector.hh b/src/celeritas/mucf/executor/detail/DDChannelSelector.hh index a6c645b988..b3cf8453af 100644 --- a/src/celeritas/mucf/executor/detail/DDChannelSelector.hh +++ b/src/celeritas/mucf/executor/detail/DDChannelSelector.hh @@ -17,9 +17,10 @@ namespace detail /*! * Select final channel for muonic dd molecules. * - * The branching ratio is temperature dependent and determines the probability + * The material temperature determines the probability * of the outcome of the fusion ending in the \f$ ^3\text{He} \f$ channels - * versus the tritium channel. + * versus the tritium channel. If the outcome is a \f$ ^3\text{He} \f$ channel, + * a constant sticking fraction of 12.2% is used to define if sticking occurs. */ class DDChannelSelector { @@ -37,7 +38,7 @@ class DDChannelSelector inline CELER_FUNCTION Channel operator()(Engine& rng); private: - real_type he3_probability_{}; + real_type he3_probability_; // Constant sticking fraction between the two 3He channels inline CELER_FUNCTION real_type sticking_fraction() const { return 0.122; } @@ -58,11 +59,11 @@ DDChannelSelector::DDChannelSelector(real_type material_temperature) CELER_EXPECT(material_temperature > 0); real_type branching_ratio{0}; - if (material_temperature < 50) + if (material_temperature <= 50) { branching_ratio = 1; } - else if (material_temperature < 100) + else if (material_temperature <= 100) { branching_ratio = 1.0088 * (material_temperature - 50); } diff --git a/test/celeritas/CMakeLists.txt b/test/celeritas/CMakeLists.txt index 44082c1e21..163c974e6e 100644 --- a/test/celeritas/CMakeLists.txt +++ b/test/celeritas/CMakeLists.txt @@ -150,6 +150,7 @@ celeritas_add_test(decay/MuDecay.test.cc ${_needs_double}) #-----------------------------------------------------------------------------# # MuCF +celeritas_add_test(mucf/ChannelSelector.test.cc ${_needs_double}) celeritas_add_test(mucf/MuonicAtomSelector.test.cc ${_needs_double}) celeritas_add_test(mucf/DTMucfInteractor.test.cc ${_needs_double}) diff --git a/test/celeritas/mucf/ChannelSelector.test.cc b/test/celeritas/mucf/ChannelSelector.test.cc new file mode 100644 index 0000000000..b8ac7caf63 --- /dev/null +++ b/test/celeritas/mucf/ChannelSelector.test.cc @@ -0,0 +1,314 @@ +//------------------------------- -*- C++ -*- -------------------------------// +// Copyright Celeritas contributors: see top-level COPYRIGHT file for details +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/mucf/ChannelSelector.test.cc +//---------------------------------------------------------------------------// +#include "corecel/cont/Range.hh" +#include "corecel/random/DiagnosticRngEngine.hh" +#include "celeritas/mucf/executor/detail/DDChannelSelector.hh" +#include "celeritas/mucf/executor/detail/DTChannelSelector.hh" +#include "celeritas/mucf/executor/detail/TTChannelSelector.hh" + +#include "celeritas_test.hh" + +namespace celeritas +{ +namespace detail +{ +namespace test +{ +using namespace ::celeritas::test; +//---------------------------------------------------------------------------// +// TEST HARNESS +//---------------------------------------------------------------------------// + +class ChannelSelectorTest : public Test +{ + protected: + using Engine = DiagnosticRngEngine; + + void SetUp() override { rng_.reset_count(); } + + Engine& rng() { return rng_; } + + // Sticking fraction between the two dd --> 3He channels + real_type const dd_sticking_fraction() { return 0.122; } + // Sticking fraction for dt + real_type const dt_sticking_fraction() { return 0.00857; } + // Sticking fraction for tt + real_type const tt_sticking_fraction() { return 0.14; } + + // Calculate dd --> 3He channel probability from the branching ratio + real_type he3_probability(real_type branching_ratio) + { + return branching_ratio / (branching_ratio + 1); + } + + // Calculate sigma for the statistical tests + real_type calc_sigma(real_type num_samples, real_type success_prob) + { + return std::sqrt(num_samples * success_prob * (1 - success_prob)); + } + + private: + Engine rng_; +}; + +//---------------------------------------------------------------------------// +// TESTS +//---------------------------------------------------------------------------// + +TEST_F(ChannelSelectorTest, dd_channel_low_temperature) +{ + // dd fusion at T < 50 K: branching_ratio = 1 + real_type const temperature = 30.0; + real_type const branching_ratio = 1.0; + real_type const he3_probability = this->he3_probability(branching_ratio); + real_type const sticking_fraction = this->dd_sticking_fraction(); + + DDChannelSelector select_channel(temperature); + + int num_samples = 100000; + int helium3_count = 0; + int muonichelium3_count = 0; + int tritium_count = 0; + + for ([[maybe_unused]] auto i : range(num_samples)) + { + auto channel = select_channel(this->rng()); + if (channel == DDMucfInteractor::Channel::helium3_muon_neutron) + { + helium3_count++; + } + else if (channel == DDMucfInteractor::Channel::muonichelium3_neutron) + { + muonichelium3_count++; + } + else if (channel == DDMucfInteractor::Channel::tritium_muon_proton) + { + tritium_count++; + } + else + { + FAIL() << "Unexpected channel selected"; + } + } + + EXPECT_EQ(num_samples, helium3_count + muonichelium3_count + tritium_count); + + real_type expected_tritium_count = num_samples * (1 - he3_probability); + real_type expected_helium3_count = num_samples * he3_probability + * (1 - sticking_fraction); + real_type expected_muonichelium3_count = num_samples * he3_probability + * sticking_fraction; + // 3 sigma tolerance + real_type tolerance = 3 * this->calc_sigma(num_samples, he3_probability); + + EXPECT_NEAR(expected_tritium_count, tritium_count, tolerance); + EXPECT_NEAR(expected_helium3_count, helium3_count, tolerance); + EXPECT_NEAR(expected_muonichelium3_count, muonichelium3_count, tolerance); +} + +//---------------------------------------------------------------------------// +TEST_F(ChannelSelectorTest, dd_channel_mid_temperature) +{ + // DD fusion at 50 < T < 100 K: branching_ratio = 1.0088 * (T - 50) = 25.22 + real_type const temperature = 75.0; + real_type const branching_ratio = 1.0088 * (temperature - 50); + real_type const he3_probability = this->he3_probability(branching_ratio); + + DDChannelSelector select_channel(temperature); + + size_type const num_samples = 10000; + size_type he3_total_count = 0; + size_type tritium_count = 0; + + for ([[maybe_unused]] auto i : range(num_samples)) + { + auto channel = select_channel(this->rng()); + if (channel == DDMucfInteractor::Channel::helium3_muon_neutron + || channel == DDMucfInteractor::Channel::muonichelium3_neutron) + { + he3_total_count++; + } + else if (channel == DDMucfInteractor::Channel::tritium_muon_proton) + { + tritium_count++; + } + } + + EXPECT_EQ(num_samples, he3_total_count + tritium_count); + + real_type const expected_he3_count = num_samples * he3_probability; + real_type const tolerance + = 3 * this->calc_sigma(num_samples, he3_probability); + + EXPECT_NEAR(expected_he3_count, he3_total_count, tolerance); +} + +//---------------------------------------------------------------------------// +TEST_F(ChannelSelectorTest, dd_channel_high_temperature) +{ + // DD fusion at T < 300 K: branching_ratio = 1.44 + real_type const temperature = 300; + real_type const branching_ratio = 1.44; + real_type const he3_probability = this->he3_probability(branching_ratio); + + DDChannelSelector select_channel(temperature); + + size_type const num_samples = 10000; + size_type he3_total_count = 0; + size_type tritium_count = 0; + + for ([[maybe_unused]] auto i : range(num_samples)) + { + auto channel = select_channel(this->rng()); + if (channel == DDMucfInteractor::Channel::helium3_muon_neutron + || channel == DDMucfInteractor::Channel::muonichelium3_neutron) + { + he3_total_count++; + } + else if (channel == DDMucfInteractor::Channel::tritium_muon_proton) + { + tritium_count++; + } + } + + EXPECT_EQ(num_samples, he3_total_count + tritium_count); + + real_type const expected_he3_count = num_samples * he3_probability; + // 3 sigma tolerance + real_type const tolerance + = 3 * this->calc_sigma(num_samples, he3_probability); + + EXPECT_NEAR(expected_he3_count, he3_total_count, tolerance); +} + +//---------------------------------------------------------------------------// +TEST_F(ChannelSelectorTest, dd_sticking_fraction_within_he3) +{ + // Test that when He3 channel is selected, sticking fraction is 12.2% + real_type const temperature = 300; + real_type const sticking_fraction = this->dd_sticking_fraction(); + + DDChannelSelector select_channel(temperature); + + size_type const num_samples = 10000; + size_type helium3_count = 0; + size_type muonichelium3_count = 0; + + // Only count He3 channels + for ([[maybe_unused]] auto i : range(num_samples)) + { + auto channel = select_channel(this->rng()); + if (channel == DDMucfInteractor::Channel::helium3_muon_neutron) + { + helium3_count++; + } + else if (channel == DDMucfInteractor::Channel::muonichelium3_neutron) + { + muonichelium3_count++; + } + } + + size_type total_he3 = helium3_count + muonichelium3_count; + EXPECT_GT(total_he3, 0); + + real_type expected_muonichelium3 = total_he3 * sticking_fraction; + // 3 sigma tolerance + real_type const tolerance + = 3 * this->calc_sigma(total_he3, sticking_fraction); + + EXPECT_NEAR(expected_muonichelium3, muonichelium3_count, tolerance); +} + +//---------------------------------------------------------------------------// +TEST_F(ChannelSelectorTest, dt_channel) +{ + // DT fusion: ~0.8% sticking + DTChannelSelector select_channel; + + size_type const num_samples = 100000; + size_type alpha_count = 0; + size_type muonicalpha_count = 0; + + for ([[maybe_unused]] auto i : range(num_samples)) + { + auto channel = select_channel(this->rng()); + if (channel == DTMucfInteractor::Channel::alpha_muon_neutron) + { + alpha_count++; + } + else if (channel == DTMucfInteractor::Channel::muonicalpha_neutron) + { + muonicalpha_count++; + } + else + { + FAIL() << "Unexpected channel selected"; + } + } + + EXPECT_EQ(num_samples, alpha_count + muonicalpha_count); + + real_type const sticking_fraction = this->dt_sticking_fraction(); + real_type const expected_muonicalpha_count = num_samples + * sticking_fraction; + real_type const expected_alpha_count = num_samples + * (1 - sticking_fraction); + // 3 sigma tolerance + real_type const tolerance + = 3 * this->calc_sigma(num_samples, sticking_fraction); + + EXPECT_NEAR(expected_muonicalpha_count, muonicalpha_count, tolerance); + EXPECT_NEAR(expected_alpha_count, alpha_count, tolerance); +} + +//---------------------------------------------------------------------------// +TEST_F(ChannelSelectorTest, tt_channel) +{ + // TT fusion: ~14% sticking + TTChannelSelector select_channel; + + size_type const num_samples = 10000; + size_type alpha_count = 0; + size_type muonicalpha_count = 0; + + for ([[maybe_unused]] auto i : range(num_samples)) + { + auto channel = select_channel(this->rng()); + if (channel == TTMucfInteractor::Channel::alpha_muon_neutron_neutron) + { + alpha_count++; + } + else if (channel + == TTMucfInteractor::Channel::muonicalpha_neutron_neutron) + { + muonicalpha_count++; + } + else + { + FAIL() << "Unexpected channel selected"; + } + } + + EXPECT_EQ(num_samples, alpha_count + muonicalpha_count); + + real_type const sticking_fraction = this->tt_sticking_fraction(); + real_type const expected_muonicalpha_count = num_samples + * sticking_fraction; + real_type const expected_alpha_count = num_samples + * (1 - sticking_fraction); + + real_type const tolerance + = 3 * this->calc_sigma(num_samples, sticking_fraction); + + EXPECT_NEAR(expected_muonicalpha_count, muonicalpha_count, tolerance); + EXPECT_NEAR(expected_alpha_count, alpha_count, tolerance); +} + +//---------------------------------------------------------------------------// +} // namespace test +} // namespace detail +} // namespace celeritas diff --git a/test/celeritas/mucf/MuonicAtomSelector.test.cc b/test/celeritas/mucf/MuonicAtomSelector.test.cc index 7d1d39911e..14a6443cbc 100644 --- a/test/celeritas/mucf/MuonicAtomSelector.test.cc +++ b/test/celeritas/mucf/MuonicAtomSelector.test.cc @@ -27,7 +27,7 @@ using namespace ::celeritas::test; class MuonicAtomSelectorTest : public Test { protected: - using Engine = ::celeritas::test::DiagnosticRngEngine; + using Engine = DiagnosticRngEngine; void SetUp() override { rng_.reset_count(); } From 2048637ec642962c42cebf31f128592363934d7a Mon Sep 17 00:00:00 2001 From: Stefano Tognini Date: Wed, 28 Jan 2026 16:16:05 -0500 Subject: [PATCH 22/42] Expand documentation --- doc/implementation/mucf-physics.rst | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/doc/implementation/mucf-physics.rst b/doc/implementation/mucf-physics.rst index ff014d715d..97b9ca751c 100644 --- a/doc/implementation/mucf-physics.rst +++ b/doc/implementation/mucf-physics.rst @@ -133,14 +133,17 @@ host/device data via :cpp:class:`celeritas::detail::MucfMaterialInserter`. .. doxygenclass:: celeritas::detail::MucfMaterialInserter The main cycle is managed by the model's -:cpp:class:`celeritas::DTMixMucfExecutor`. The muonic atom selection is handled -by detail classes, while the Interactors are reserved for sampling final states -of the outgoing secondaries. +:cpp:class:`celeritas::DTMixMucfExecutor`. The muonic atom, spin, fusion channel +selections are handled by detail classes, while the Interactors are reserved for +sampling final states of the outgoing secondaries. .. note:: Only reactive channels are implemented. .. doxygenclass:: celeritas::detail::MuonicAtomSelector .. doxygenclass:: celeritas::detail::MuonicAtomSpinSelector +.. doxygenclass:: celeritas::detail::DDChannelSelector +.. doxygenclass:: celeritas::detail::DTChannelSelector +.. doxygenclass:: celeritas::detail::TTChannelSelector .. doxygenclass:: celeritas::DDMucfInteractor .. doxygenclass:: celeritas::DTMucfInteractor .. doxygenclass:: celeritas::TTMucfInteractor From 0e593ad432ae4a266bae221c8e3a09d3dbfd74fe Mon Sep 17 00:00:00 2001 From: Stefano Tognini Date: Thu, 29 Jan 2026 09:06:29 -0500 Subject: [PATCH 23/42] Add missing reference --- src/celeritas/mucf/executor/detail/DDChannelSelector.hh | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/celeritas/mucf/executor/detail/DDChannelSelector.hh b/src/celeritas/mucf/executor/detail/DDChannelSelector.hh index b3cf8453af..09d20698fe 100644 --- a/src/celeritas/mucf/executor/detail/DDChannelSelector.hh +++ b/src/celeritas/mucf/executor/detail/DDChannelSelector.hh @@ -20,7 +20,8 @@ namespace detail * The material temperature determines the probability * of the outcome of the fusion ending in the \f$ ^3\text{He} \f$ channels * versus the tritium channel. If the outcome is a \f$ ^3\text{He} \f$ channel, - * a constant sticking fraction of 12.2% is used to define if sticking occurs. + * a constant sticking fraction of 12.2% is used to define if sticking occurs + * [ \todo https://doi.org/10.1134/S106377961102002X ]. */ class DDChannelSelector { From 0758cb6770431123265347ed0bea6f06a5872de3 Mon Sep 17 00:00:00 2001 From: Stefano Tognini Date: Thu, 19 Feb 2026 13:59:00 -0500 Subject: [PATCH 24/42] Form molecule, define final cycle time, and set track time accordingly --- src/celeritas/mucf/data/DTMixMucfData.hh | 10 +- .../mucf/executor/DTMixMucfExecutor.hh | 43 ++----- .../detail/DTMixMuonicMoleculeSelector.hh | 108 +++++++++++++++--- .../mucf/model/detail/MucfMaterialInserter.cc | 42 +++---- .../mucf/model/detail/MucfMaterialInserter.hh | 10 +- test/celeritas/mucf/DTMixMucfModel.test.cc | 15 +-- 6 files changed, 143 insertions(+), 85 deletions(-) diff --git a/src/celeritas/mucf/data/DTMixMucfData.hh b/src/celeritas/mucf/data/DTMixMucfData.hh index ff6aeb7c66..2e9d277b6f 100644 --- a/src/celeritas/mucf/data/DTMixMucfData.hh +++ b/src/celeritas/mucf/data/DTMixMucfData.hh @@ -108,7 +108,7 @@ struct DTMixMucfData template using MaterialItems = Collection; using GridRecord = NonuniformGridRecord; - using CycleTimesArray = EnumArray>; + using CycleRatesArray = EnumArray>; using MaterialFractionsArray = EnumArray; //! Particles @@ -126,8 +126,8 @@ struct DTMixMucfData MaterialItems mucfmatid_to_matid; //! Isotopic fractions per material: [mat_comp_id][isotope] MaterialItems isotopic_fractions; - //! Cycle times per material: [mat_comp_id][muonic_molecule][spin_index] - MaterialItems cycle_times; //!< In [s] + //! Cycle rates per material: [mat_comp_id][muonic_molecule][spin_index] + MaterialItems cycle_rates; //!< In [1/s] //! \todo Add mean atom spin flip times //! \todo Add mean atom transfer times //!@} @@ -137,7 +137,7 @@ struct DTMixMucfData { return particle_ids && particle_masses && muon_energy_cdf && !reals.empty() && !mucfmatid_to_matid.empty() - && !isotopic_fractions.empty() && !cycle_times.empty(); + && !isotopic_fractions.empty() && !cycle_rates.empty(); } //! Assign from another set of data @@ -153,7 +153,7 @@ struct DTMixMucfData this->reals = other.reals; this->mucfmatid_to_matid = other.mucfmatid_to_matid; this->isotopic_fractions = other.isotopic_fractions; - this->cycle_times = other.cycle_times; + this->cycle_rates = other.cycle_rates; return *this; } diff --git a/src/celeritas/mucf/executor/DTMixMucfExecutor.hh b/src/celeritas/mucf/executor/DTMixMucfExecutor.hh index 5431f59b55..f604dd5b19 100644 --- a/src/celeritas/mucf/executor/DTMixMucfExecutor.hh +++ b/src/celeritas/mucf/executor/DTMixMucfExecutor.hh @@ -48,11 +48,6 @@ DTMixMucfExecutor::operator()(celeritas::CoreTrackView const& track) auto element = mat_record.element_record(elcomp_id); CELER_ASSERT(element.atomic_number() == AtomicNumber{1}); // Must be H - auto rng = track.rng(); - - // Muon decay may compete against other "actions" in this executor - real_type const decay_len{}; //! \todo Set muon decay interaction length - // Find muCF material ID from PhysMatId // Make this a View if ever used beyond this executor auto find = [&](PhysMatId matid) -> MuCfMatId { @@ -71,13 +66,12 @@ DTMixMucfExecutor::operator()(celeritas::CoreTrackView const& track) auto const mucf_matid = find(track.material().material_id()); CELER_ASSERT(mucf_matid); - // Form d or t muonic atom - detail::MuonicAtomSelector form_atom( - data.isotopic_fractions[mucf_matid][MucfIsotope::deuterium]); - auto muonic_atom = form_atom(rng); + auto rng = track.rng(); - // Select atom spin via a helper class - detail::MuonicAtomSpinSelector select_atom_spin(muonic_atom); + // Form d or t muonic atom + auto muonic_atom = detail::MuonicAtomSelector( + data.isotopic_fractions[mucf_matid][MucfIsotope::deuterium])(rng); + auto atom_spin = detail::MuonicAtomSpinSelector(muonic_atom)(rng); // { // Competing at-rest processes which add to the total track time @@ -86,31 +80,16 @@ DTMixMucfExecutor::operator()(celeritas::CoreTrackView const& track) // } // Form dd, dt, or tt muonic molecule - detail::DTMixMuonicMoleculeSelector form_muonic_molecule; - auto muonic_molecule = form_muonic_molecule(rng); - - // Select molecule spin - detail::MuonicMoleculeSpinSelector select_molecule_spin(muonic_molecule); - auto const molecule_spin = select_molecule_spin(rng); - - // Load cycle time for the selected molecule - auto const cycle_time - = data.cycle_times[mucf_matid][muonic_molecule][molecule_spin]; - CELER_ASSERT(cycle_time > 0); - - // Check if muon decays before fusion happens - real_type const mucf_len = cycle_time * track.sim().step_length(); - if (decay_len < mucf_len) - { - // Muon decays and halts the interaction - //! \todo Update track time and return muon decay interactor - } + auto [muonic_molecule, cycle_time] = detail::DTMixMuonicMoleculeSelector( + muonic_atom, + atom_spin, + data.isotopic_fractions[mucf_matid], + data.cycle_rates[mucf_matid])(rng); - //! \todo Correct track time update? Or should be done in Interactors? + // Update track time according to the sampled cycle time track.sim().add_time(cycle_time); // Fuse molecule and generate secondaries - //! \todo Maybe move the channel selectors into the interactors auto allocate_secondaries = phys_step_view.make_secondary_allocator(); Interaction result; switch (muonic_molecule) diff --git a/src/celeritas/mucf/executor/detail/DTMixMuonicMoleculeSelector.hh b/src/celeritas/mucf/executor/detail/DTMixMuonicMoleculeSelector.hh index cc4c7bd30b..e42a96720d 100644 --- a/src/celeritas/mucf/executor/detail/DTMixMuonicMoleculeSelector.hh +++ b/src/celeritas/mucf/executor/detail/DTMixMuonicMoleculeSelector.hh @@ -6,7 +6,9 @@ //---------------------------------------------------------------------------// #pragma once +#include "corecel/random/distribution/ExponentialDistribution.hh" #include "celeritas/Types.hh" +#include "celeritas/mucf/Types.hh" namespace celeritas { @@ -23,12 +25,38 @@ namespace detail class DTMixMuonicMoleculeSelector { public: - //! Construct with args; \todo Update documentation - inline CELER_FUNCTION DTMixMuonicMoleculeSelector(/* args */); + //!@{ + //! \name Type aliases + using CycleRatesArray = EnumArray>; + using MaterialFractionsArray = EnumArray; + using HalfSpinInt = units::HalfSpinInt; + //!@} + + // Result of molecule and cycle time selection + struct Result + { + MucfMuonicMolecule molecule; + real_type cycle_time; // [s] + }; + + //! Construct with muonic atom and material information + inline CELER_FUNCTION + DTMixMuonicMoleculeSelector(MucfMuonicAtom atom, + HalfSpinInt spin, + MaterialFractionsArray material_fractions, + CycleRatesArray cycle_rates); // Select muonic molecule template - inline CELER_FUNCTION MucfMuonicMolecule operator()(Engine&); + inline CELER_FUNCTION Result operator()(Engine&); + + private: + MucfMuonicAtom atom_; + HalfSpinInt spin_; + MaterialFractionsArray material_fractions_; + CycleRatesArray cycle_rates_; + + inline CELER_FUNCTION size_type spin_state_index(HalfSpinInt spin) const; }; //---------------------------------------------------------------------------// @@ -40,26 +68,80 @@ class DTMixMuonicMoleculeSelector * \todo Update documentation */ CELER_FUNCTION -DTMixMuonicMoleculeSelector::DTMixMuonicMoleculeSelector(/* args */) +DTMixMuonicMoleculeSelector::DTMixMuonicMoleculeSelector( + MucfMuonicAtom atom, + HalfSpinInt spin, + MaterialFractionsArray material_fractions, + CycleRatesArray cycle_rates) + : atom_(atom) + , spin_(spin) + , material_fractions_(material_fractions) + , cycle_rates_(cycle_rates) { - //! \todo Implement - CELER_NOT_IMPLEMENTED("Mucf muonic molecule selection"); + CELER_EXPECT(atom < MucfMuonicAtom::size_); + CELER_EXPECT(material_fractions_[MucfIsotope::deuterium] + + material_fractions_[MucfIsotope::tritium] + > 0); } //---------------------------------------------------------------------------// /*! - * Return selected muonic molecule. + * Return a muonic molecule by selecting the shortest sampled molecule + * interaction time (\f$ T = -\ln(r) \times \tau_\text{cycle} \f$). */ template -CELER_FUNCTION MucfMuonicMolecule -DTMixMuonicMoleculeSelector::operator()(Engine&) +CELER_FUNCTION DTMixMuonicMoleculeSelector::Result +DTMixMuonicMoleculeSelector::operator()(Engine& rng) { - MucfMuonicMolecule result{MucfMuonicMolecule::size_}; + using MMM = MucfMuonicMolecule; - //! \todo Implement + real_type shortest_time = std::numeric_limits::max(); + MucfMuonicMolecule molecule{MMM::size_}; + for (auto mol : range(MMM::size_)) + { + auto const& cycle_rate + = cycle_rates_[mol][this->spin_state_index(spin_)]; + + auto sampled_time = ExponentialDistribution(cycle_rate)(rng); + if (sampled_time < shortest_time) + { + shortest_time = sampled_time; + molecule = mol; + } + } + + CELER_ENSURE(molecule < MMM::size_); + return {molecule, shortest_time}; +} + +//---------------------------------------------------------------------------// +/*! + * Return correct \c CycleRatesArray spin state index for each molecule. + * + * Given the limited number of states, this is manually hardcoded during data + * construction by \c MucfMaterialInserter::calc_[dd|dt|tt]_cycle . + */ +CELER_FUNCTION size_type +DTMixMuonicMoleculeSelector::spin_state_index(HalfSpinInt spin) const +{ + switch (atom_) + { + case MucfMuonicAtom::deuterium: { + // Muonic deuterium has spin states 1/2 and 3/2, which correspond + // to indices 0 and 1, respectively + CELER_EXPECT(spin == HalfSpinInt{1} || spin == HalfSpinInt{3}); + return (spin == HalfSpinInt{1}) ? 0 : 1; + } - CELER_ENSURE(result < MucfMuonicMolecule::size_); - return result; + case MucfMuonicAtom::tritium: { + // Muonic tritium has spin states 0 and 1, which correspond to + // indices 0 and 1, respectively + CELER_EXPECT(spin == HalfSpinInt{0} || spin == HalfSpinInt{1}); + return (spin == HalfSpinInt{0}) ? 0 : 1; + } + default: + CELER_ASSERT_UNREACHABLE(); + } } //---------------------------------------------------------------------------// diff --git a/src/celeritas/mucf/model/detail/MucfMaterialInserter.cc b/src/celeritas/mucf/model/detail/MucfMaterialInserter.cc index f4d9930249..e6a7eaac1b 100644 --- a/src/celeritas/mucf/model/detail/MucfMaterialInserter.cc +++ b/src/celeritas/mucf/model/detail/MucfMaterialInserter.cc @@ -20,12 +20,12 @@ MucfMaterialInserter::MucfMaterialInserter(HostVal* host_data, inp::MucfPhysics const& data) : mucfmatid_to_matid_(&host_data->mucfmatid_to_matid) , isotopic_fractions_(&host_data->isotopic_fractions) - , cycle_times_(&host_data->cycle_times) + , cycle_rates_(&host_data->cycle_rates) , data_(data) { CELER_EXPECT(data_); - // Initialize interpolators for cycle time tables + // Initialize interpolators for cycle rate tables for (auto const& cycle_data : data_.cycle_rates) { // Use emplace to avoid copy/move of InterpolatorHelper objects @@ -43,7 +43,7 @@ bool MucfMaterialInserter::operator()(MaterialView const& material) using LhdArray = EquilibrateDensitiesSolver::LhdArray; MaterialFractionsArray isotopic_fractions; - CycleTimesArray cycle_times; + CycleRatesArray cycle_rates; LhdArray lhd_densities{}; auto from_mass_number = [&](AtomicMassNumber mass) -> MucfIsotope { @@ -94,19 +94,19 @@ bool MucfMaterialInserter::operator()(MaterialView const& material) if (lhd_densities[MucfIsotope::deuterium]) { - cycle_times[MucfMuonicMolecule::deuterium_deuterium] + cycle_rates[MucfMuonicMolecule::deuterium_deuterium] = this->calc_dd_cycle(equilibrium_densities, material.temperature()); } if (lhd_densities[MucfIsotope::tritium]) { - cycle_times[MucfMuonicMolecule::tritium_tritium] = this->calc_tt_cycle( + cycle_rates[MucfMuonicMolecule::tritium_tritium] = this->calc_tt_cycle( equilibrium_densities, material.temperature()); } if (lhd_densities[MucfIsotope::deuterium] && lhd_densities[MucfIsotope::tritium]) { - cycle_times[MucfMuonicMolecule::deuterium_tritium] + cycle_rates[MucfMuonicMolecule::deuterium_tritium] = this->calc_dt_cycle(equilibrium_densities, material.temperature()); } @@ -114,7 +114,7 @@ bool MucfMaterialInserter::operator()(MaterialView const& material) // Add muCF material to the model's host/device data mucfmatid_to_matid_.push_back(material.material_id()); isotopic_fractions_.push_back(std::move(isotopic_fractions)); - cycle_times_.push_back(std::move(cycle_times)); + cycle_rates_.push_back(std::move(cycle_rates)); //! \todo Store mean atom spin flip and transfer times @@ -123,7 +123,7 @@ bool MucfMaterialInserter::operator()(MaterialView const& material) //---------------------------------------------------------------------------// /*! - * Calculate dd muonic molecules cycle times. + * Calculate dd muonic molecules cycle rates. * * F = 1/2 and F = 3/2 are the reactive spin states for dd fusion. */ @@ -143,10 +143,8 @@ MucfMaterialInserter::calc_dd_cycle(EquilibriumArray const& eq_dens, = this->interpolator(CTT::deuterium_deuterium, HalfSpinInt{3}); MoleculeCycles result; - result[0] = real_type{1} - / (dd_dens * dd_1_over_2_interpolate(temperature)); // F = 1/2 - result[1] = real_type{1} - / (dd_dens * dd_3_over_2_interpolate(temperature)); // F = 3/2 + result[0] = (dd_dens * dd_1_over_2_interpolate(temperature)); // F = 1/2 + result[1] = (dd_dens * dd_3_over_2_interpolate(temperature)); // F = 3/2 CELER_ENSURE(result[0] >= 0 && result[1] >= 0); return result; @@ -154,7 +152,7 @@ MucfMaterialInserter::calc_dd_cycle(EquilibriumArray const& eq_dens, //---------------------------------------------------------------------------// /*! - * Calculate dt muonic molecules cycle times. + * Calculate dt muonic molecules cycle rates. * * F = 0 and F = 1 are the reactive spin states for dt fusion. */ @@ -189,14 +187,12 @@ MucfMaterialInserter::calc_dt_cycle(EquilibriumArray const& eq_dens, // Interpolate over rates, store final cycle time (1/rate) MoleculeCycles result; - result[0] = real_type{1} - / (hd_dens * hd0_interpolate(temperature) - + dd_dens * dd0_interpolate(temperature) - + dt_dens * dt0_interpolate(temperature)); // F = 0 - result[1] = real_type{1} - / (hd_dens * hd1_interpolate(temperature) - + dd_dens * dd1_interpolate(temperature) - + dt_dens * dt1_interpolate(temperature)); // F = 1 + result[0] = hd_dens * hd0_interpolate(temperature) + + dd_dens * dd0_interpolate(temperature) + + dt_dens * dt0_interpolate(temperature); // F = 0 + result[1] = hd_dens * hd1_interpolate(temperature) + + dd_dens * dd1_interpolate(temperature) + + dt_dens * dt1_interpolate(temperature); // F = 1 CELER_ENSURE(result[0] >= 0 && result[1] >= 0); return result; @@ -204,7 +200,7 @@ MucfMaterialInserter::calc_dt_cycle(EquilibriumArray const& eq_dens, //---------------------------------------------------------------------------// /*! - * Calculate tt muonic molecules cycle times. + * Calculate tt muonic molecules cycle rates. * * F = 1/2 is the only reactive spin state for tt fusion. */ @@ -221,7 +217,7 @@ MucfMaterialInserter::calc_tt_cycle(EquilibriumArray const& eq_dens, = this->interpolator(CTT::tritium_tritium, HalfSpinInt{1}); MoleculeCycles result; - result[0] = real_type{1} / (tt_dens * tt_interpolate(temperature)); + result[0] = tt_dens * tt_interpolate(temperature); // F = 1/2 CELER_ENSURE(result[0] >= 0 && result[1] == 0); return result; diff --git a/src/celeritas/mucf/model/detail/MucfMaterialInserter.hh b/src/celeritas/mucf/model/detail/MucfMaterialInserter.hh index 113605bb57..076c481b7e 100644 --- a/src/celeritas/mucf/model/detail/MucfMaterialInserter.hh +++ b/src/celeritas/mucf/model/detail/MucfMaterialInserter.hh @@ -47,7 +47,7 @@ class MucfMaterialInserter private: using MoleculeCycles = Array; - using CycleTimesArray = EnumArray; + using CycleRatesArray = EnumArray; using EquilibriumArray = EquilibrateDensitiesSolver::EquilibriumArray; using MaterialFractionsArray = EnumArray; using AtomicMassNumber = AtomicNumber; @@ -61,7 +61,7 @@ class MucfMaterialInserter CollectionBuilder mucfmatid_to_matid_; CollectionBuilder isotopic_fractions_; - CollectionBuilder cycle_times_; + CollectionBuilder cycle_rates_; // Const data std::map const mass_isotope_map_{ {AtomicMassNumber{1}, MucfIsotope::protium}, @@ -73,15 +73,15 @@ class MucfMaterialInserter //// HELPER FUNCTIONS //// - // Calculate mean fusion cycle times for dd muonic molecules + // Calculate mean fusion cycle rates for dd muonic molecules Array calc_dd_cycle(EquilibriumArray const& eq_dens, real_type const temperature); - // Calculate mean fusion cycle times for dt muonic molecules + // Calculate mean fusion cycle rates for dt muonic molecules Array calc_dt_cycle(EquilibriumArray const& eq_dens, real_type const temperature); - // Calculate mean fusion cycle times for tt muonic molecules + // Calculate mean fusion cycle rates for tt muonic molecules Array calc_tt_cycle(EquilibriumArray const& eq_dens, real_type const temperature); diff --git a/test/celeritas/mucf/DTMixMucfModel.test.cc b/test/celeritas/mucf/DTMixMucfModel.test.cc index f51f74f28c..19767f3d39 100644 --- a/test/celeritas/mucf/DTMixMucfModel.test.cc +++ b/test/celeritas/mucf/DTMixMucfModel.test.cc @@ -83,22 +83,23 @@ TEST_F(DTMixMucfModelTest, data) // TT static double const expected_tt_1_over_2_cycle_time{1.4056833511329384e-06}; - auto const& cycles = data.cycle_times; + auto const& rates = data.cycle_rates; // DD cycle times EXPECT_SOFT_EQ(expected_dd_1_over_2_cycle_time, - cycles[MuCfMatId{0}][Molecule::deuterium_deuterium][0]); + 1 / rates[MuCfMatId{0}][Molecule::deuterium_deuterium][0]); EXPECT_SOFT_EQ(expected_dd_3_over_2_cycle_time, - cycles[MuCfMatId{0}][Molecule::deuterium_deuterium][1]); + 1 / rates[MuCfMatId{0}][Molecule::deuterium_deuterium][1]); // DT cycle times EXPECT_SOFT_EQ(expected_dt_0_cycle_time, - cycles[MuCfMatId{0}][Molecule::deuterium_tritium][0]); + 1 / rates[MuCfMatId{0}][Molecule::deuterium_tritium][0]); EXPECT_SOFT_EQ(expected_dt_1_cycle_time, - cycles[MuCfMatId{0}][Molecule::deuterium_tritium][1]); + 1 / rates[MuCfMatId{0}][Molecule::deuterium_tritium][1]); // TT cycle times EXPECT_SOFT_EQ(expected_tt_1_over_2_cycle_time, - cycles[MuCfMatId{0}][Molecule::tritium_tritium][0]); - EXPECT_SOFT_EQ(0, cycles[MuCfMatId{0}][Molecule::tritium_tritium][1]); + 1 / rates[MuCfMatId{0}][Molecule::tritium_tritium][0]); + EXPECT_SOFT_EQ(numeric_limits::infinity(), + 1 / rates[MuCfMatId{0}][Molecule::tritium_tritium][1]); } //---------------------------------------------------------------------------// From d91a3047267de621e3b1005fbbf4f3c6671cb5eb Mon Sep 17 00:00:00 2001 From: Stefano Tognini Date: Thu, 19 Feb 2026 19:03:07 -0500 Subject: [PATCH 25/42] Add test; Refactor molecule selector; wip: fix array index selection --- .../mucf/executor/DTMixMucfExecutor.hh | 12 +- .../detail/DTMixMuonicMoleculeSelector.hh | 149 ---------------- .../executor/detail/MuonicMoleculeSelector.hh | 167 ++++++++++++++++++ .../detail/MuonicMoleculeSpinSelector.hh | 64 ------- test/celeritas/CMakeLists.txt | 1 + .../mucf/MuonicMoleculeSelector.test.cc | 132 ++++++++++++++ 6 files changed, 305 insertions(+), 220 deletions(-) delete mode 100644 src/celeritas/mucf/executor/detail/DTMixMuonicMoleculeSelector.hh create mode 100644 src/celeritas/mucf/executor/detail/MuonicMoleculeSelector.hh delete mode 100644 src/celeritas/mucf/executor/detail/MuonicMoleculeSpinSelector.hh create mode 100644 test/celeritas/mucf/MuonicMoleculeSelector.test.cc diff --git a/src/celeritas/mucf/executor/DTMixMucfExecutor.hh b/src/celeritas/mucf/executor/DTMixMucfExecutor.hh index f604dd5b19..d92736b984 100644 --- a/src/celeritas/mucf/executor/DTMixMucfExecutor.hh +++ b/src/celeritas/mucf/executor/DTMixMucfExecutor.hh @@ -16,10 +16,9 @@ #include "detail/DDChannelSelector.hh" #include "detail/DTChannelSelector.hh" -#include "detail/DTMixMuonicMoleculeSelector.hh" #include "detail/MuonicAtomSelector.hh" #include "detail/MuonicAtomSpinSelector.hh" -#include "detail/MuonicMoleculeSpinSelector.hh" +#include "detail/MuonicMoleculeSelector.hh" #include "detail/TTChannelSelector.hh" namespace celeritas @@ -48,6 +47,8 @@ DTMixMucfExecutor::operator()(celeritas::CoreTrackView const& track) auto element = mat_record.element_record(elcomp_id); CELER_ASSERT(element.atomic_number() == AtomicNumber{1}); // Must be H + //! \todo Make sure that at this point we selected d or t already + // Find muCF material ID from PhysMatId // Make this a View if ever used beyond this executor auto find = [&](PhysMatId matid) -> MuCfMatId { @@ -80,11 +81,8 @@ DTMixMucfExecutor::operator()(celeritas::CoreTrackView const& track) // } // Form dd, dt, or tt muonic molecule - auto [muonic_molecule, cycle_time] = detail::DTMixMuonicMoleculeSelector( - muonic_atom, - atom_spin, - data.isotopic_fractions[mucf_matid], - data.cycle_rates[mucf_matid])(rng); + auto [muonic_molecule, cycle_time] = detail::MuonicMoleculeSelector( + muonic_atom, atom_spin, data.cycle_rates[mucf_matid])(rng); // Update track time according to the sampled cycle time track.sim().add_time(cycle_time); diff --git a/src/celeritas/mucf/executor/detail/DTMixMuonicMoleculeSelector.hh b/src/celeritas/mucf/executor/detail/DTMixMuonicMoleculeSelector.hh deleted file mode 100644 index e42a96720d..0000000000 --- a/src/celeritas/mucf/executor/detail/DTMixMuonicMoleculeSelector.hh +++ /dev/null @@ -1,149 +0,0 @@ -//------------------------------- -*- C++ -*- -------------------------------// -// Copyright Celeritas contributors: see top-level COPYRIGHT file for details -// SPDX-License-Identifier: (Apache-2.0 OR MIT) -//---------------------------------------------------------------------------// -//! \file celeritas/mucf/executor/detail/DTMixMuonicMoleculeSelector.hh -//---------------------------------------------------------------------------// -#pragma once - -#include "corecel/random/distribution/ExponentialDistribution.hh" -#include "celeritas/Types.hh" -#include "celeritas/mucf/Types.hh" - -namespace celeritas -{ -namespace detail -{ -//---------------------------------------------------------------------------// -/*! - * Select a muonic molecule by calculating the interaction lengths of the - * possible molecule formations. - * - * This is the equivalent of Geant4's - * \c G4VRestProcess::AtRestGetPhysicalInteractionLength - */ -class DTMixMuonicMoleculeSelector -{ - public: - //!@{ - //! \name Type aliases - using CycleRatesArray = EnumArray>; - using MaterialFractionsArray = EnumArray; - using HalfSpinInt = units::HalfSpinInt; - //!@} - - // Result of molecule and cycle time selection - struct Result - { - MucfMuonicMolecule molecule; - real_type cycle_time; // [s] - }; - - //! Construct with muonic atom and material information - inline CELER_FUNCTION - DTMixMuonicMoleculeSelector(MucfMuonicAtom atom, - HalfSpinInt spin, - MaterialFractionsArray material_fractions, - CycleRatesArray cycle_rates); - - // Select muonic molecule - template - inline CELER_FUNCTION Result operator()(Engine&); - - private: - MucfMuonicAtom atom_; - HalfSpinInt spin_; - MaterialFractionsArray material_fractions_; - CycleRatesArray cycle_rates_; - - inline CELER_FUNCTION size_type spin_state_index(HalfSpinInt spin) const; -}; - -//---------------------------------------------------------------------------// -// INLINE DEFINITIONS -//---------------------------------------------------------------------------// -/*! - * Construct with args. - * - * \todo Update documentation - */ -CELER_FUNCTION -DTMixMuonicMoleculeSelector::DTMixMuonicMoleculeSelector( - MucfMuonicAtom atom, - HalfSpinInt spin, - MaterialFractionsArray material_fractions, - CycleRatesArray cycle_rates) - : atom_(atom) - , spin_(spin) - , material_fractions_(material_fractions) - , cycle_rates_(cycle_rates) -{ - CELER_EXPECT(atom < MucfMuonicAtom::size_); - CELER_EXPECT(material_fractions_[MucfIsotope::deuterium] - + material_fractions_[MucfIsotope::tritium] - > 0); -} - -//---------------------------------------------------------------------------// -/*! - * Return a muonic molecule by selecting the shortest sampled molecule - * interaction time (\f$ T = -\ln(r) \times \tau_\text{cycle} \f$). - */ -template -CELER_FUNCTION DTMixMuonicMoleculeSelector::Result -DTMixMuonicMoleculeSelector::operator()(Engine& rng) -{ - using MMM = MucfMuonicMolecule; - - real_type shortest_time = std::numeric_limits::max(); - MucfMuonicMolecule molecule{MMM::size_}; - for (auto mol : range(MMM::size_)) - { - auto const& cycle_rate - = cycle_rates_[mol][this->spin_state_index(spin_)]; - - auto sampled_time = ExponentialDistribution(cycle_rate)(rng); - if (sampled_time < shortest_time) - { - shortest_time = sampled_time; - molecule = mol; - } - } - - CELER_ENSURE(molecule < MMM::size_); - return {molecule, shortest_time}; -} - -//---------------------------------------------------------------------------// -/*! - * Return correct \c CycleRatesArray spin state index for each molecule. - * - * Given the limited number of states, this is manually hardcoded during data - * construction by \c MucfMaterialInserter::calc_[dd|dt|tt]_cycle . - */ -CELER_FUNCTION size_type -DTMixMuonicMoleculeSelector::spin_state_index(HalfSpinInt spin) const -{ - switch (atom_) - { - case MucfMuonicAtom::deuterium: { - // Muonic deuterium has spin states 1/2 and 3/2, which correspond - // to indices 0 and 1, respectively - CELER_EXPECT(spin == HalfSpinInt{1} || spin == HalfSpinInt{3}); - return (spin == HalfSpinInt{1}) ? 0 : 1; - } - - case MucfMuonicAtom::tritium: { - // Muonic tritium has spin states 0 and 1, which correspond to - // indices 0 and 1, respectively - CELER_EXPECT(spin == HalfSpinInt{0} || spin == HalfSpinInt{1}); - return (spin == HalfSpinInt{0}) ? 0 : 1; - } - default: - CELER_ASSERT_UNREACHABLE(); - } -} - -//---------------------------------------------------------------------------// -} // namespace detail -} // namespace celeritas diff --git a/src/celeritas/mucf/executor/detail/MuonicMoleculeSelector.hh b/src/celeritas/mucf/executor/detail/MuonicMoleculeSelector.hh new file mode 100644 index 0000000000..ef855262ea --- /dev/null +++ b/src/celeritas/mucf/executor/detail/MuonicMoleculeSelector.hh @@ -0,0 +1,167 @@ +//------------------------------- -*- C++ -*- -------------------------------// +// Copyright Celeritas contributors: see top-level COPYRIGHT file for details +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/mucf/executor/detail/MuonicMoleculeSelector.hh +//---------------------------------------------------------------------------// +#pragma once + +#include "corecel/cont/EnumArray.hh" +#include "corecel/cont/Range.hh" +#include "corecel/random/distribution/ExponentialDistribution.hh" +#include "celeritas/Quantities.hh" +#include "celeritas/Types.hh" +#include "celeritas/mucf/Types.hh" + +namespace celeritas +{ +namespace detail +{ +//---------------------------------------------------------------------------// +/*! + * Select a muonic molecule and its final cycle time. + * + * The muonic molecule is selected by sampling the shortest fusion cycle time + * (\f$ T = -\ln(r) \times \tau_\text{cycle} \f$). + * + * \note The cycle rates cached in the model data are calculated based on the + * material definition, and thus are propagated to the final cycle time sampled + * here and used in the fusion process. + */ +class MuonicMoleculeSelector +{ + public: + //!@{ + //! \name Type aliases + using CycleRatesArray = EnumArray>; + using HalfSpinInt = units::HalfSpinInt; + //!@} + + // Result of molecule and cycle time selection + struct Result + { + MucfMuonicMolecule molecule; + real_type cycle_time; //!< in [s] + }; + + //! Construct with muonic atom and material information + inline CELER_FUNCTION MuonicMoleculeSelector(MucfMuonicAtom atom, + HalfSpinInt spin, + CycleRatesArray cycle_rates); + + // Select muonic molecule and sample its final cycle time + template + inline CELER_FUNCTION Result operator()(Engine& rng); + + private: + MucfMuonicAtom atom_; + HalfSpinInt spin_; + CycleRatesArray cycle_rates_; + + // Return the correct array index for each reactive spin state + inline CELER_FUNCTION size_type spin_state_index(MucfMuonicMolecule mol, + HalfSpinInt spin) const; + + // Sample the final cycle time for a given molecule + template + inline CELER_FUNCTION real_type sample_exp_time(MucfMuonicMolecule molecule, + Engine& rng); +}; + +//---------------------------------------------------------------------------// +// INLINE DEFINITIONS +//---------------------------------------------------------------------------// +/*! + * Construct with muonic atom and material information. + */ +CELER_FUNCTION +MuonicMoleculeSelector::MuonicMoleculeSelector(MucfMuonicAtom atom, + HalfSpinInt spin, + CycleRatesArray cycle_rates) + : atom_(atom), spin_(spin), cycle_rates_(cycle_rates) +{ + CELER_EXPECT(atom < MucfMuonicAtom::size_); +} + +//---------------------------------------------------------------------------// +/*! + * Return a muonic molecule and its sampled cycle time. + */ +template +CELER_FUNCTION MuonicMoleculeSelector::Result +MuonicMoleculeSelector::operator()(Engine& rng) +{ + using MMM = MucfMuonicMolecule; + using CycleTimeArray = EnumArray; + + // Select final molecule and its respective cycle time + CycleTimeArray sampled_times; + real_type shortest_time = std::numeric_limits::infinity(); + MMM molecule{MMM::size_}; + for (auto mol : range(MMM::size_)) + { + sampled_times[mol] = this->sample_exp_time(mol, rng); + if (sampled_times[mol] < shortest_time) + { + shortest_time = sampled_times[mol]; + molecule = mol; + } + } + + CELER_ENSURE(molecule < MMM::size_); + return {molecule, shortest_time}; +} + +//---------------------------------------------------------------------------// +/*! + * Return correct \c CycleRatesArray spin state index for each molecule. + * + * Given the limited number of states, this is manually hardcoded during data + * construction by \c MucfMaterialInserter::calc_[dd|dt|tt]_cycle . + */ +CELER_FUNCTION size_type MuonicMoleculeSelector::spin_state_index( + MucfMuonicMolecule mol, HalfSpinInt spin) const +{ + // The array index is wonky and I have to fix that... + switch (atom_) + { + case MucfMuonicAtom::deuterium: { + // Muonic deuterium has spin states 1/2 and 3/2, which correspond + // to indices 0 and 1, respectively + CELER_EXPECT(spin == HalfSpinInt{1} || spin == HalfSpinInt{3}); + return (spin == HalfSpinInt{1}) ? 0 : 1; + } + + case MucfMuonicAtom::tritium: { + // Muonic tritium has spin states 0 and 1, which correspond to + // indices 0 and 1, respectively + CELER_EXPECT(spin == HalfSpinInt{0} || spin == HalfSpinInt{1}); + return (spin == HalfSpinInt{0}) ? 0 : 1; + } + default: + CELER_ASSERT_UNREACHABLE(); + } +} + +//---------------------------------------------------------------------------// +/*! + * Sample the final exponential distribution for a given cycle rate using + * \f$ T = -\ln(r) \times \tau_\text{cycle} \f$ . If the cached rate is zero, + * the returned time is set to infinity, which effectively removes that + * molecule as a possible selection. + */ +template +CELER_FUNCTION real_type MuonicMoleculeSelector::sample_exp_time( + MucfMuonicMolecule molecule, Engine& rng) +{ + CELER_EXPECT(molecule < MucfMuonicMolecule::size_); + auto const rate = cycle_rates_[molecule][this->spin_state_index(spin_)]; + + // Return infinity when the cached rate is zero + return (rate > 0) ? ExponentialDistribution(rate)(rng) + : std::numeric_limits::infinity(); +} + +//---------------------------------------------------------------------------// +} // namespace detail +} // namespace celeritas diff --git a/src/celeritas/mucf/executor/detail/MuonicMoleculeSpinSelector.hh b/src/celeritas/mucf/executor/detail/MuonicMoleculeSpinSelector.hh deleted file mode 100644 index 67ae8d4e29..0000000000 --- a/src/celeritas/mucf/executor/detail/MuonicMoleculeSpinSelector.hh +++ /dev/null @@ -1,64 +0,0 @@ -//------------------------------- -*- C++ -*- -------------------------------// -// Copyright Celeritas contributors: see top-level COPYRIGHT file for details -// SPDX-License-Identifier: (Apache-2.0 OR MIT) -//---------------------------------------------------------------------------// -//! \file celeritas/mucf/executor/detail/MuonicMoleculeSpinSelector.hh -//---------------------------------------------------------------------------// -#pragma once - -#include "corecel/Assert.hh" -#include "celeritas/Types.hh" - -namespace celeritas -{ -namespace detail -{ -//---------------------------------------------------------------------------// -/*! - * Select muonic molecule spin, in units of \f$ \frac{\hbar}{2} \f$. - * - * Applicable to \f$ (dd)_\mu \f$, * \f$ (dt)_\mu \f$ , and \f$ (tt)_\mu \f$. - */ -class MuonicMoleculeSpinSelector -{ - public: - // Construct with muonic molecule type - inline CELER_FUNCTION - MuonicMoleculeSpinSelector(MucfMuonicMolecule molecule); - - // Sample and return a spin value in units of hbar / 2 - template - inline CELER_FUNCTION size_type operator()(Engine&); - - private: - MucfMuonicMolecule molecule_; -}; - -//---------------------------------------------------------------------------// -// INLINE DEFINITIONS -//---------------------------------------------------------------------------// -/*! - * Construct with muonic molecule. - */ -CELER_FUNCTION MuonicMoleculeSpinSelector::MuonicMoleculeSpinSelector( - MucfMuonicMolecule molecule) - : molecule_(molecule) -{ - CELER_EXPECT(molecule_ < MucfMuonicMolecule::size_); - CELER_NOT_IMPLEMENTED("Mucf muonic molecule spin selection"); -} - -//---------------------------------------------------------------------------// -/*! - * Return selected spin, in units of \f$ \hbar / 2 \f$. - */ -template -CELER_FUNCTION size_type MuonicMoleculeSpinSelector::operator()(Engine&) -{ - //! \todo switch on molecule_ - CELER_ASSERT_UNREACHABLE(); -} - -//---------------------------------------------------------------------------// -} // namespace detail -} // namespace celeritas diff --git a/test/celeritas/CMakeLists.txt b/test/celeritas/CMakeLists.txt index 1998471b81..659881cf92 100644 --- a/test/celeritas/CMakeLists.txt +++ b/test/celeritas/CMakeLists.txt @@ -156,6 +156,7 @@ celeritas_add_test(mucf/DDMucfInteractor.test.cc) celeritas_add_test(mucf/DTMucfInteractor.test.cc) celeritas_add_test(mucf/DTMixMucfModel.test.cc) celeritas_add_test(mucf/MuonicAtomSelector.test.cc) +celeritas_add_test(mucf/MuonicMoleculeSelector.test.cc) #-----------------------------------------------------------------------------# # EM diff --git a/test/celeritas/mucf/MuonicMoleculeSelector.test.cc b/test/celeritas/mucf/MuonicMoleculeSelector.test.cc new file mode 100644 index 0000000000..b3d2bb3b69 --- /dev/null +++ b/test/celeritas/mucf/MuonicMoleculeSelector.test.cc @@ -0,0 +1,132 @@ +//------------------------------- -*- C++ -*- -------------------------------// +// Copyright Celeritas contributors: see top-level COPYRIGHT file for details +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/mucf/MuonicMoleculeSelector.test.cc +//---------------------------------------------------------------------------// +#include "celeritas/mucf/executor/detail/MuonicMoleculeSelector.hh" + +#include "corecel/cont/Range.hh" +#include "corecel/io/Logger.hh" +#include "corecel/random/DiagnosticRngEngine.hh" + +#include "MucfInteractorHostTestBase.hh" +#include "celeritas_test.hh" + +namespace celeritas +{ +namespace detail +{ +namespace test +{ +using namespace ::celeritas::test; +//---------------------------------------------------------------------------// +// TEST HARNESS +//---------------------------------------------------------------------------// + +class MuonicMoleculeSelectorTest + : public ::celeritas::test::MucfInteractorHostTestBase +{ + protected: + using Base = ::celeritas::test::MucfInteractorHostTestBase; + using Engine = DiagnosticRngEngine; + using MMM = MucfMuonicMolecule; + using MMA = MucfMuonicAtom; + using HalfSpinInt = MuonicMoleculeSelector::HalfSpinInt; + using CycleRatesArray = MuonicMoleculeSelector::CycleRatesArray; + using MoleculeCountArray = EnumArray; + + void SetUp() override + { + data_ = this->host_data(); + rng_.reset_count(); + + auto const mucfmatid = MuCfMatId{0}; + cycle_rates_ = data_.cycle_rates[mucfmatid]; + } + + MoleculeCountArray sample_molecules(MucfMuonicAtom atom, + HalfSpinInt spin, + CycleRatesArray cycle_rates, + size_type num_samples) + { + MuonicMoleculeSelector select_molecule(atom, spin, cycle_rates); + + MoleculeCountArray result; + for ([[maybe_unused]] auto i : range(num_samples)) + { + auto sampled = select_molecule(rng_); + EXPECT_LT(sampled.molecule, MucfMuonicMolecule::size_); + result[sampled.molecule]++; + } + return result; + } + + protected: + HostCRef data_; + Engine rng_; + CycleRatesArray cycle_rates_; +}; + +//---------------------------------------------------------------------------// +// TESTS +//---------------------------------------------------------------------------// +TEST_F(MuonicMoleculeSelectorTest, deuterium_spin_selection) +{ + size_type const num_samples{1}; + CycleRatesArray rates; + + //// DD molecule ///// + { + // Make DT and TT rates zero + rates[MMM::deuterium_tritium] = {0, 0}; + rates[MMM::tritium_tritium] = {0, 0}; + + // Test DD, F = 1/2 + rates[MMM::deuterium_deuterium] = {1e6, 0}; + auto const dd_1_result = MuonicMoleculeSelector( + MMA::deuterium, HalfSpinInt{1}, rates)(rng_); + + EXPECT_EQ(dd_1_result.molecule, MMM::deuterium_deuterium); + EXPECT_SOFT_EQ(dd_1_result.cycle_time, 1.9989533630516e-06); + + // Test DD, F = 3/2 + rates[MMM::deuterium_deuterium] = {0, 1e6}; + auto const dd_3_result = MuonicMoleculeSelector( + MMA::deuterium, HalfSpinInt{3}, rates)(rng_); + + EXPECT_EQ(dd_3_result.molecule, MMM::deuterium_deuterium); + EXPECT_SOFT_EQ(dd_3_result.cycle_time, 1.8031326676554e-07); + } + + //// DT molecule ///// + { + // Make DD and TT rates zero + rates[MMM::deuterium_deuterium] = {0, 0}; + rates[MMM::tritium_tritium] = {0, 0}; + + // Test DT, F = 0 + rates[MMM::deuterium_tritium] = {1e6, 0}; + auto const dt_0_result = MuonicMoleculeSelector( + MMA::deuterium, HalfSpinInt{0}, rates)(rng_); + + EXPECT_EQ(dt_0_result.molecule, MMM::deuterium_tritium); + EXPECT_SOFT_EQ(dt_0_result.cycle_time, 1.9989533630516e-06); + + // Test DT, F = 1 + rates[MMM::deuterium_tritium] = {0, 1e6}; + auto const dt_1_result = MuonicMoleculeSelector( + MMA::deuterium, HalfSpinInt{1}, rates)(rng_); + + EXPECT_EQ(dt_1_result.molecule, MMM::deuterium_tritium); + EXPECT_SOFT_EQ(dt_1_result.cycle_time, 1.8031326676554e-07); + } +} + +//---------------------------------------------------------------------------// +TEST_F(MuonicMoleculeSelectorTest, tritium_spin_selection) {} + +//---------------------------------------------------------------------------// +} // namespace test +} // namespace detail +} // namespace celeritas From 4fed2c6d7fd4ffe8f5cc638c8ad6f66d474d7293 Mon Sep 17 00:00:00 2001 From: Stefano Tognini Date: Fri, 20 Feb 2026 12:59:13 -0500 Subject: [PATCH 26/42] Finish test harness; squash bugs; refactor `MuCfMatId` -> `MucfMatId` --- src/celeritas/mucf/Types.hh | 2 +- src/celeritas/mucf/data/DTMixMucfData.hh | 4 +- .../mucf/executor/DTMixMucfExecutor.hh | 6 +- .../executor/detail/MuonicMoleculeSelector.hh | 131 +++++++++------- .../mucf/model/detail/MucfMaterialInserter.hh | 6 +- test/celeritas/mucf/DTMixMucfModel.test.cc | 16 +- .../mucf/MuonicMoleculeSelector.test.cc | 147 +++++++++++------- 7 files changed, 181 insertions(+), 131 deletions(-) diff --git a/src/celeritas/mucf/Types.hh b/src/celeritas/mucf/Types.hh index 23a02637c0..c986fcf900 100644 --- a/src/celeritas/mucf/Types.hh +++ b/src/celeritas/mucf/Types.hh @@ -55,7 +55,7 @@ enum class MucfMuonicMolecule //---------------------------------------------------------------------------// //! Opaque index of a muCF material component -using MuCfMatId = OpaqueId; +using MucfMatId = OpaqueId; //---------------------------------------------------------------------------// } // namespace celeritas diff --git a/src/celeritas/mucf/data/DTMixMucfData.hh b/src/celeritas/mucf/data/DTMixMucfData.hh index 2e9d277b6f..255effccf6 100644 --- a/src/celeritas/mucf/data/DTMixMucfData.hh +++ b/src/celeritas/mucf/data/DTMixMucfData.hh @@ -106,7 +106,7 @@ struct DTMixMucfData template using Items = Collection; template - using MaterialItems = Collection; + using MaterialItems = Collection; using GridRecord = NonuniformGridRecord; using CycleRatesArray = EnumArray>; using MaterialFractionsArray = EnumArray; @@ -122,7 +122,7 @@ struct DTMixMucfData //!@{ //! Material-dependent data calculated at model construction - //! \c PhysMatId indexed by \c MuCfMatId + //! \c PhysMatId indexed by \c MucfMatId MaterialItems mucfmatid_to_matid; //! Isotopic fractions per material: [mat_comp_id][isotope] MaterialItems isotopic_fractions; diff --git a/src/celeritas/mucf/executor/DTMixMucfExecutor.hh b/src/celeritas/mucf/executor/DTMixMucfExecutor.hh index d92736b984..000f7ed235 100644 --- a/src/celeritas/mucf/executor/DTMixMucfExecutor.hh +++ b/src/celeritas/mucf/executor/DTMixMucfExecutor.hh @@ -51,18 +51,18 @@ DTMixMucfExecutor::operator()(celeritas::CoreTrackView const& track) // Find muCF material ID from PhysMatId // Make this a View if ever used beyond this executor - auto find = [&](PhysMatId matid) -> MuCfMatId { + auto find = [&](PhysMatId matid) -> MucfMatId { CELER_EXPECT(matid); for (auto i : range(data.mucfmatid_to_matid.size())) { - if (auto const comp_id = MuCfMatId{i}; + if (auto const comp_id = MucfMatId{i}; data.mucfmatid_to_matid[comp_id] == matid) { return comp_id; } } // MuCF material ID not found - return MuCfMatId{}; + return MucfMatId{}; }; auto const mucf_matid = find(track.material().material_id()); CELER_ASSERT(mucf_matid); diff --git a/src/celeritas/mucf/executor/detail/MuonicMoleculeSelector.hh b/src/celeritas/mucf/executor/detail/MuonicMoleculeSelector.hh index ef855262ea..aaea0244fe 100644 --- a/src/celeritas/mucf/executor/detail/MuonicMoleculeSelector.hh +++ b/src/celeritas/mucf/executor/detail/MuonicMoleculeSelector.hh @@ -40,8 +40,15 @@ class MuonicMoleculeSelector // Result of molecule and cycle time selection struct Result { - MucfMuonicMolecule molecule; - real_type cycle_time; //!< in [s] + MucfMuonicMolecule molecule{MucfMuonicMolecule::size_}; + real_type cycle_time{std::numeric_limits::max()}; + + //! Check whether the data are assigned + explicit operator bool() const + { + return molecule < MucfMuonicMolecule::size_ && cycle_time > 0 + && cycle_time < std::numeric_limits::max(); + } }; //! Construct with muonic atom and material information @@ -55,13 +62,9 @@ class MuonicMoleculeSelector private: MucfMuonicAtom atom_; - HalfSpinInt spin_; + size_type cycle_rate_index_; CycleRatesArray cycle_rates_; - // Return the correct array index for each reactive spin state - inline CELER_FUNCTION size_type spin_state_index(MucfMuonicMolecule mol, - HalfSpinInt spin) const; - // Sample the final cycle time for a given molecule template inline CELER_FUNCTION real_type sample_exp_time(MucfMuonicMolecule molecule, @@ -72,15 +75,48 @@ class MuonicMoleculeSelector // INLINE DEFINITIONS //---------------------------------------------------------------------------// /*! - * Construct with muonic atom and material information. + * Construct with muonic atom and cycle rate information. + * + * The correct cycle rate array index is determined at construction by the + * combination of atom+spin. Given the limited number of states, the position + * in the index is set in \c MucfMaterialInserter::calc_[dd|dt|tt]_cycle + * manually. + * + * \note Current implementation has no safe method to access cycle rate data + * for each spin state. */ CELER_FUNCTION MuonicMoleculeSelector::MuonicMoleculeSelector(MucfMuonicAtom atom, HalfSpinInt spin, CycleRatesArray cycle_rates) - : atom_(atom), spin_(spin), cycle_rates_(cycle_rates) + : atom_(atom), cycle_rates_(cycle_rates) { CELER_EXPECT(atom < MucfMuonicAtom::size_); + + // Check that the spin value is valid for the given atom type and set the + // cycle rate array index accordingly + switch (atom_) + { + case MucfMuonicAtom::deuterium: { + // Related to DD states F = 1/2 and F = 3/2 + CELER_EXPECT(spin == HalfSpinInt{1} || spin == HalfSpinInt{3}); + // F = 1/2 and F = 3/2 correspond to indices 0 and 1, respectively + cycle_rate_index_ = (spin == HalfSpinInt{1}) ? 0 : 1; + break; + } + case MucfMuonicAtom::tritium: { + // Related to DT states F = 0 and F = 1, and TT state F = 1/2 + CELER_EXPECT(spin == HalfSpinInt{0} || spin == HalfSpinInt{2} + || spin == HalfSpinInt{1}); + // DT: F = 0 and F = 1 correspond to indices 0 and 1, respectively + // TT: F = 1/2 corresponds to index 0 + cycle_rate_index_ + = (spin == HalfSpinInt{0} || spin == HalfSpinInt{1}) ? 0 : 1; + break; + } + default: + CELER_ASSERT_UNREACHABLE(); + } } //---------------------------------------------------------------------------// @@ -94,72 +130,57 @@ MuonicMoleculeSelector::operator()(Engine& rng) using MMM = MucfMuonicMolecule; using CycleTimeArray = EnumArray; - // Select final molecule and its respective cycle time CycleTimeArray sampled_times; - real_type shortest_time = std::numeric_limits::infinity(); - MMM molecule{MMM::size_}; - for (auto mol : range(MMM::size_)) + auto const inf = std::numeric_limits::max(); + if (atom_ == MucfMuonicAtom::deuterium) { - sampled_times[mol] = this->sample_exp_time(mol, rng); - if (sampled_times[mol] < shortest_time) - { - shortest_time = sampled_times[mol]; - molecule = mol; - } + // DD fusion is only triggered by a muonic deuterium + sampled_times[MMM::deuterium_deuterium] + = this->sample_exp_time(MMM::deuterium_deuterium, rng); + sampled_times[MMM::deuterium_tritium] = inf; + sampled_times[MMM::tritium_tritium] = inf; } - - CELER_ENSURE(molecule < MMM::size_); - return {molecule, shortest_time}; -} - -//---------------------------------------------------------------------------// -/*! - * Return correct \c CycleRatesArray spin state index for each molecule. - * - * Given the limited number of states, this is manually hardcoded during data - * construction by \c MucfMaterialInserter::calc_[dd|dt|tt]_cycle . - */ -CELER_FUNCTION size_type MuonicMoleculeSelector::spin_state_index( - MucfMuonicMolecule mol, HalfSpinInt spin) const -{ - // The array index is wonky and I have to fix that... - switch (atom_) + if (atom_ == MucfMuonicAtom::tritium) { - case MucfMuonicAtom::deuterium: { - // Muonic deuterium has spin states 1/2 and 3/2, which correspond - // to indices 0 and 1, respectively - CELER_EXPECT(spin == HalfSpinInt{1} || spin == HalfSpinInt{3}); - return (spin == HalfSpinInt{1}) ? 0 : 1; - } + // DT and TT fusions are triggered by a muonic tritium + sampled_times[MMM::deuterium_deuterium] = inf; + sampled_times[MMM::deuterium_tritium] + = this->sample_exp_time(MMM::deuterium_tritium, rng); + sampled_times[MMM::tritium_tritium] + = this->sample_exp_time(MMM::tritium_tritium, rng); + } - case MucfMuonicAtom::tritium: { - // Muonic tritium has spin states 0 and 1, which correspond to - // indices 0 and 1, respectively - CELER_EXPECT(spin == HalfSpinInt{0} || spin == HalfSpinInt{1}); - return (spin == HalfSpinInt{0}) ? 0 : 1; + Result result; + for (auto mol : range(MMM::size_)) + { + if (sampled_times[mol] < result.cycle_time) + { + result.cycle_time = sampled_times[mol]; + result.molecule = mol; } - default: - CELER_ASSERT_UNREACHABLE(); } + + CELER_ENSURE(result); + return result; } //---------------------------------------------------------------------------// /*! * Sample the final exponential distribution for a given cycle rate using - * \f$ T = -\ln(r) \times \tau_\text{cycle} \f$ . If the cached rate is zero, - * the returned time is set to infinity, which effectively removes that - * molecule as a possible selection. + * \f$ T = -\ln(r) \times \tau_\text{cycle} \f$ . If the cached rate is + * zero, the returned time is set to infinity, which effectively removes + * that molecule as a possible selection. */ template CELER_FUNCTION real_type MuonicMoleculeSelector::sample_exp_time( MucfMuonicMolecule molecule, Engine& rng) { CELER_EXPECT(molecule < MucfMuonicMolecule::size_); - auto const rate = cycle_rates_[molecule][this->spin_state_index(spin_)]; + auto const rate = cycle_rates_[molecule][cycle_rate_index_]; - // Return infinity when the cached rate is zero + // Return an infinite cycle time when the cached rate is zero return (rate > 0) ? ExponentialDistribution(rate)(rng) - : std::numeric_limits::infinity(); + : std::numeric_limits::max(); } //---------------------------------------------------------------------------// diff --git a/src/celeritas/mucf/model/detail/MucfMaterialInserter.hh b/src/celeritas/mucf/model/detail/MucfMaterialInserter.hh index 076c481b7e..f5ed5cb1ec 100644 --- a/src/celeritas/mucf/model/detail/MucfMaterialInserter.hh +++ b/src/celeritas/mucf/model/detail/MucfMaterialInserter.hh @@ -58,10 +58,10 @@ class MucfMaterialInserter //// DATA //// // DTMixMucfModel host data references populated by operator() - CollectionBuilder mucfmatid_to_matid_; - CollectionBuilder + CollectionBuilder mucfmatid_to_matid_; + CollectionBuilder isotopic_fractions_; - CollectionBuilder cycle_rates_; + CollectionBuilder cycle_rates_; // Const data std::map const mass_isotope_map_{ {AtomicMassNumber{1}, MucfIsotope::protium}, diff --git a/test/celeritas/mucf/DTMixMucfModel.test.cc b/test/celeritas/mucf/DTMixMucfModel.test.cc index 19767f3d39..8701727bef 100644 --- a/test/celeritas/mucf/DTMixMucfModel.test.cc +++ b/test/celeritas/mucf/DTMixMucfModel.test.cc @@ -69,9 +69,9 @@ TEST_F(DTMixMucfModelTest, data) // Check sotopic fractions // Single material with 50/50 d and t fractions EXPECT_SOFT_EQ( - 0.5, data.isotopic_fractions[MuCfMatId{0}][MucfIsotope::deuterium]); + 0.5, data.isotopic_fractions[MucfMatId{0}][MucfIsotope::deuterium]); EXPECT_SOFT_EQ( - 0.5, data.isotopic_fractions[MuCfMatId{0}][MucfIsotope::tritium]); + 0.5, data.isotopic_fractions[MucfMatId{0}][MucfIsotope::tritium]); // Cycle times are in seconds // DD (reactivity of F = 3/2 is almost negligible, with huge cycle times) @@ -87,19 +87,19 @@ TEST_F(DTMixMucfModelTest, data) // DD cycle times EXPECT_SOFT_EQ(expected_dd_1_over_2_cycle_time, - 1 / rates[MuCfMatId{0}][Molecule::deuterium_deuterium][0]); + 1 / rates[MucfMatId{0}][Molecule::deuterium_deuterium][0]); EXPECT_SOFT_EQ(expected_dd_3_over_2_cycle_time, - 1 / rates[MuCfMatId{0}][Molecule::deuterium_deuterium][1]); + 1 / rates[MucfMatId{0}][Molecule::deuterium_deuterium][1]); // DT cycle times EXPECT_SOFT_EQ(expected_dt_0_cycle_time, - 1 / rates[MuCfMatId{0}][Molecule::deuterium_tritium][0]); + 1 / rates[MucfMatId{0}][Molecule::deuterium_tritium][0]); EXPECT_SOFT_EQ(expected_dt_1_cycle_time, - 1 / rates[MuCfMatId{0}][Molecule::deuterium_tritium][1]); + 1 / rates[MucfMatId{0}][Molecule::deuterium_tritium][1]); // TT cycle times EXPECT_SOFT_EQ(expected_tt_1_over_2_cycle_time, - 1 / rates[MuCfMatId{0}][Molecule::tritium_tritium][0]); + 1 / rates[MucfMatId{0}][Molecule::tritium_tritium][0]); EXPECT_SOFT_EQ(numeric_limits::infinity(), - 1 / rates[MuCfMatId{0}][Molecule::tritium_tritium][1]); + 1 / rates[MucfMatId{0}][Molecule::tritium_tritium][1]); } //---------------------------------------------------------------------------// diff --git a/test/celeritas/mucf/MuonicMoleculeSelector.test.cc b/test/celeritas/mucf/MuonicMoleculeSelector.test.cc index b3d2bb3b69..a86cfae799 100644 --- a/test/celeritas/mucf/MuonicMoleculeSelector.test.cc +++ b/test/celeritas/mucf/MuonicMoleculeSelector.test.cc @@ -9,6 +9,9 @@ #include "corecel/cont/Range.hh" #include "corecel/io/Logger.hh" #include "corecel/random/DiagnosticRngEngine.hh" +#include "celeritas/mucf/Types.hh" +#include "celeritas/mucf/executor/detail/MuonicAtomSelector.hh" +#include "celeritas/mucf/executor/detail/MuonicAtomSpinSelector.hh" #include "MucfInteractorHostTestBase.hh" #include "celeritas_test.hh" @@ -34,97 +37,123 @@ class MuonicMoleculeSelectorTest using MMA = MucfMuonicAtom; using HalfSpinInt = MuonicMoleculeSelector::HalfSpinInt; using CycleRatesArray = MuonicMoleculeSelector::CycleRatesArray; - using MoleculeCountArray = EnumArray; void SetUp() override { data_ = this->host_data(); rng_.reset_count(); - auto const mucfmatid = MuCfMatId{0}; - cycle_rates_ = data_.cycle_rates[mucfmatid]; - } - - MoleculeCountArray sample_molecules(MucfMuonicAtom atom, - HalfSpinInt spin, - CycleRatesArray cycle_rates, - size_type num_samples) - { - MuonicMoleculeSelector select_molecule(atom, spin, cycle_rates); - - MoleculeCountArray result; - for ([[maybe_unused]] auto i : range(num_samples)) - { - auto sampled = select_molecule(rng_); - EXPECT_LT(sampled.molecule, MucfMuonicMolecule::size_); - result[sampled.molecule]++; - } - return result; + mucfmatid_ = MucfMatId{0}; + cycle_rates_ = data_.cycle_rates[mucfmatid_]; } protected: HostCRef data_; Engine rng_; + MucfMatId mucfmatid_; CycleRatesArray cycle_rates_; }; //---------------------------------------------------------------------------// // TESTS //---------------------------------------------------------------------------// -TEST_F(MuonicMoleculeSelectorTest, deuterium_spin_selection) +TEST_F(MuonicMoleculeSelectorTest, dd_spin_selection) { - size_type const num_samples{1}; CycleRatesArray rates; - //// DD molecule ///// - { - // Make DT and TT rates zero - rates[MMM::deuterium_tritium] = {0, 0}; - rates[MMM::tritium_tritium] = {0, 0}; + // Make DT and TT rates zero + rates[MMM::deuterium_tritium] = {0, 0}; + rates[MMM::tritium_tritium] = {0, 0}; - // Test DD, F = 1/2 - rates[MMM::deuterium_deuterium] = {1e6, 0}; - auto const dd_1_result = MuonicMoleculeSelector( - MMA::deuterium, HalfSpinInt{1}, rates)(rng_); + // DD, F = 1/2 + rates[MMM::deuterium_deuterium] = {1e6, 0}; + auto const dd_1_result + = MuonicMoleculeSelector(MMA::deuterium, HalfSpinInt{1}, rates)(rng_); - EXPECT_EQ(dd_1_result.molecule, MMM::deuterium_deuterium); - EXPECT_SOFT_EQ(dd_1_result.cycle_time, 1.9989533630516e-06); + EXPECT_EQ(dd_1_result.molecule, MMM::deuterium_deuterium); + EXPECT_SOFT_EQ(dd_1_result.cycle_time, 1.9989533630516e-06); - // Test DD, F = 3/2 - rates[MMM::deuterium_deuterium] = {0, 1e6}; - auto const dd_3_result = MuonicMoleculeSelector( - MMA::deuterium, HalfSpinInt{3}, rates)(rng_); + // DD, F = 3/2 + rates[MMM::deuterium_deuterium] = {0, 1e6}; + auto const dd_3_result + = MuonicMoleculeSelector(MMA::deuterium, HalfSpinInt{3}, rates)(rng_); - EXPECT_EQ(dd_3_result.molecule, MMM::deuterium_deuterium); - EXPECT_SOFT_EQ(dd_3_result.cycle_time, 1.8031326676554e-07); - } + EXPECT_EQ(dd_3_result.molecule, MMM::deuterium_deuterium); + EXPECT_SOFT_EQ(dd_3_result.cycle_time, 1.8031326676554e-07); +} - //// DT molecule ///// - { - // Make DD and TT rates zero - rates[MMM::deuterium_deuterium] = {0, 0}; - rates[MMM::tritium_tritium] = {0, 0}; +//---------------------------------------------------------------------------// +TEST_F(MuonicMoleculeSelectorTest, dt_spin_selection) +{ + CycleRatesArray rates; - // Test DT, F = 0 - rates[MMM::deuterium_tritium] = {1e6, 0}; - auto const dt_0_result = MuonicMoleculeSelector( - MMA::deuterium, HalfSpinInt{0}, rates)(rng_); + // Make DD and TT rates zero + rates[MMM::deuterium_deuterium] = {0, 0}; + rates[MMM::tritium_tritium] = {0, 0}; - EXPECT_EQ(dt_0_result.molecule, MMM::deuterium_tritium); - EXPECT_SOFT_EQ(dt_0_result.cycle_time, 1.9989533630516e-06); + // DT, F = 0 + rates[MMM::deuterium_tritium] = {1e6, 0}; + auto const dt_0_result + = MuonicMoleculeSelector(MMA::tritium, HalfSpinInt{0}, rates)(rng_); - // Test DT, F = 1 - rates[MMM::deuterium_tritium] = {0, 1e6}; - auto const dt_1_result = MuonicMoleculeSelector( - MMA::deuterium, HalfSpinInt{1}, rates)(rng_); + EXPECT_EQ(dt_0_result.molecule, MMM::deuterium_tritium); + EXPECT_SOFT_EQ(dt_0_result.cycle_time, 1.9989533630516e-06); - EXPECT_EQ(dt_1_result.molecule, MMM::deuterium_tritium); - EXPECT_SOFT_EQ(dt_1_result.cycle_time, 1.8031326676554e-07); - } + // DT, F = 1 + rates[MMM::deuterium_tritium] = {0, 1e6}; + auto const dt_2_result + = MuonicMoleculeSelector(MMA::tritium, HalfSpinInt{2}, rates)(rng_); + + EXPECT_EQ(dt_2_result.molecule, MMM::deuterium_tritium); + EXPECT_SOFT_EQ(dt_2_result.cycle_time, 1.8031326676554e-07); +} + +//---------------------------------------------------------------------------// +TEST_F(MuonicMoleculeSelectorTest, tt_spin_selection) +{ + CycleRatesArray rates; + + // TT, F = 1/2 + rates[MMM::deuterium_deuterium] = {0, 0}; + rates[MMM::deuterium_tritium] = {0, 0}; + rates[MMM::tritium_tritium] = {1e6, 0}; + + auto const tt_1_result + = MuonicMoleculeSelector(MMA::tritium, HalfSpinInt{1}, rates)(rng_); + EXPECT_EQ(tt_1_result.molecule, MMM::tritium_tritium); + EXPECT_SOFT_EQ(tt_1_result.cycle_time, 1.9989533630516e-06); } //---------------------------------------------------------------------------// -TEST_F(MuonicMoleculeSelectorTest, tritium_spin_selection) {} +TEST_F(MuonicMoleculeSelectorTest, model_data) +{ + using MoleculeCountArray = EnumArray; + + // Mimic executor behavior: form atom, then form molecule + size_type const num_samples = 10000; + auto const& rates = data_.cycle_rates[mucfmatid_]; + + MoleculeCountArray molecule_counts; + for (auto i : range(num_samples)) + { + // Form muonic atom + auto muonic_atom = MuonicAtomSelector( + data_.isotopic_fractions[mucfmatid_][MucfIsotope::deuterium])(rng_); + auto atom_spin = MuonicAtomSpinSelector(muonic_atom)(rng_); + + // Form molecule that will call determine the interactor + auto result + = MuonicMoleculeSelector(muonic_atom, atom_spin, rates)(rng_); + molecule_counts[result.molecule]++; + } + + // From Acceleron's simulation with equivalent material: 50/50 d/t at 300K, + // the fraction of each fusion call is: DD: 20.1%, DT: 78.9%, TT: 0.9% + // The counts on this test (20%, 79%, 1%) are spot on with Acceleron. + EXPECT_EQ(2044, molecule_counts[MMM::deuterium_deuterium]); + EXPECT_EQ(7940, molecule_counts[MMM::deuterium_tritium]); + EXPECT_EQ(16, molecule_counts[MMM::tritium_tritium]); +} //---------------------------------------------------------------------------// } // namespace test From 4a9371944a6c5758f96fcf24a19814efe5eb8f79 Mon Sep 17 00:00:00 2001 From: Stefano Tognini Date: Fri, 20 Feb 2026 15:25:12 -0500 Subject: [PATCH 27/42] Fix comment --- test/celeritas/mucf/MuonicMoleculeSelector.test.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/celeritas/mucf/MuonicMoleculeSelector.test.cc b/test/celeritas/mucf/MuonicMoleculeSelector.test.cc index a86cfae799..24319aad2d 100644 --- a/test/celeritas/mucf/MuonicMoleculeSelector.test.cc +++ b/test/celeritas/mucf/MuonicMoleculeSelector.test.cc @@ -148,8 +148,8 @@ TEST_F(MuonicMoleculeSelectorTest, model_data) } // From Acceleron's simulation with equivalent material: 50/50 d/t at 300K, - // the fraction of each fusion call is: DD: 20.1%, DT: 78.9%, TT: 0.9% - // The counts on this test (20%, 79%, 1%) are spot on with Acceleron. + // the fraction of each fusion call is: DD: 20.1%, DT: 78.9%, TT: 0.9%. + // The counts on this test (20%, 79%, 0.2%) are close to Acceleron's. EXPECT_EQ(2044, molecule_counts[MMM::deuterium_deuterium]); EXPECT_EQ(7940, molecule_counts[MMM::deuterium_tritium]); EXPECT_EQ(16, molecule_counts[MMM::tritium_tritium]); From bea52c6b6cbd4e51d379a55776f1c1eeee165b2f Mon Sep 17 00:00:00 2001 From: Stefano Tognini Date: Fri, 20 Feb 2026 16:14:45 -0500 Subject: [PATCH 28/42] Update documentation --- doc/_static/zotero.bib | 82 +++++++++++++------ doc/implementation/mucf-physics.rst | 7 +- .../mucf/executor/detail/DDChannelSelector.hh | 28 ++----- .../mucf/executor/detail/DTChannelSelector.hh | 12 +-- .../mucf/executor/detail/TTChannelSelector.hh | 17 ++-- 5 files changed, 82 insertions(+), 64 deletions(-) diff --git a/doc/_static/zotero.bib b/doc/_static/zotero.bib index 36c84fa9a9..562a5539c7 100644 --- a/doc/_static/zotero.bib +++ b/doc/_static/zotero.bib @@ -3747,32 +3747,64 @@ @article{bom-experimentaldt-2005 } @article{yamashita-muonicspin-2022, - title = {Roles of resonant muonic molecule in new kinetics model and muon catalyzed fusion in compressed gas}, - author = {Yamashita, Takuma and Kino, Yasushi and Okutsu, Kenichi and Okada, Shinji and Sato, Motoyasu}, - year = 2022, - month = apr, - volume = {12}, - number = {1}, - journal = {Scientific Reports}, - shortjournal = {Sci Rep}, - issn = {2045-2322}, - doi = {10.1038/s41598-022-09487-0}, - url = {https://www.nature.com/articles/s41598-022-09487-0}, - pages = {6393} + title = {Roles of resonant muonic molecule in new kinetics model and muon catalyzed fusion in compressed gas}, + author = {Yamashita, Takuma and Kino, Yasushi and Okutsu, Kenichi and Okada, Shinji and Sato, Motoyasu}, + year = 2022, + month = apr, + volume = {12}, + number = {1}, + journal = {Scientific Reports}, + shortjournal = {Sci Rep}, + issn = {2045-2322}, + doi = {10.1038/s41598-022-09487-0}, + url = {https://www.nature.com/articles/s41598-022-09487-0}, + pages = {6393} } @article{faifman-mucfformation-1996, - title = {Quadrupole corrections to matrix elements of transitions in resonant reactions of muonic molecule formation}, - author = {Faifman, M. P. and Strizh, T. A. and Armour, E. A. G. and Harston, M. R.}, - year = 1996, - month = dec, - volume = {101-102}, - rights = {http://www.springer.com/tdm}, - issn = {0304-3834, 1572-9540}, - url = {http://link.springer.com/10.1007/BF02227621}, - doi = {10.1007/BF02227621}, - pages = {179--189}, - number = {1}, - journaltitle = {Hyperfine Interactions}, - shortjournal = {Hyperfine Interact}, + title = {Quadrupole corrections to matrix elements of transitions in resonant reactions of muonic molecule formation}, + author = {Faifman, M. P. and Strizh, T. A. and Armour, E. A. G. and Harston, M. R.}, + year = 1996, + month = dec, + volume = {101-102}, + rights = {http://www.springer.com/tdm}, + issn = {0304-3834, 1572-9540}, + url = {http://link.springer.com/10.1007/BF02227621}, + doi = {10.1007/BF02227621}, + pages = {179--189}, + number = {1}, + journaltitle = {Hyperfine Interactions}, + shortjournal = {Hyperfine Interact}, +} + +@article{balin-mucf-2011, + title = {{High precision study of muon catalyzed fusion in D2 and {HD} gas}}, + author = {Balin, D. V. and Ganzha, V. A. and Kozlov, S. M. and Maev, E. M. and Petrov, G. E. and Soroka, M. A. and Schapkin, G. N. and Semenchuk, G. G. and Trofimov, V. A. and Vasiliev, A. A. and Vorobyov, A. A. and Voropaev, N. I. and Petitjean, C. and Gartner, B. and Lauss, B. and Marton, J. and Zmeskal, J. and Case, T. and Crowe, K. M. and Kammel, P. and Hartmann, F. J. and Faifman, M. P.}, + volume = {42}, + rights = {http://www.springer.com/tdm}, + issn = {1063-7796, 1531-8559}, + url = {http://link.springer.com/10.1134/S106377961102002X}, + doi = {10.1134/S106377961102002X}, + year = 2011, + month = apr, + pages = {185--214}, + number = {2}, + journaltitle = {Physics of Particles and Nuclei}, + shortjournal = {Phys. Part. Nuclei}, +} + +@article{bogdanova-mucf-2009, + title = {Experimental investigation of muon-catalyzed t + t fusion}, + author = {Bogdanova, L. N. and Bom, V. R. and Demin, A. M. and Demin, D. L. and Van Eijk, C. W. E. and Filchagin, S. V. and Filchenkov, V. V. and Grafov, N. N. and Grishechkin, S. K. and Gritsaj, K. I. and Konin, A. D. and Kuryakin, A. V. and Medved’, S. V. and Musyaev, R. K. and Rudenko, A. I. and Tumkin, D. P. and Vinogradov, Yu. I. and Yukhimchuk, A. A. and Yukhimchuk, S. A. and Zinov, V. G. and Zlatoustovskii, S. V.}, + volume = {108}, + year = 2009, + month = feb, + rights = {http://www.springer.com/tdm}, + issn = {1063-7761, 1090-6509}, + url = {http://link.springer.com/10.1134/S1063776109020034}, + doi = {10.1134/S1063776109020034}, + pages = {216--225}, + number = {2}, + journaltitle = {Journal of Experimental and Theoretical Physics}, + shortjournal = {J. Exp. Theor. Phys.} } diff --git a/doc/implementation/mucf-physics.rst b/doc/implementation/mucf-physics.rst index 5b3008f166..390d3fcd40 100644 --- a/doc/implementation/mucf-physics.rst +++ b/doc/implementation/mucf-physics.rst @@ -133,14 +133,15 @@ host/device data via :cpp:class:`celeritas::detail::MucfMaterialInserter`. .. doxygenclass:: celeritas::detail::MucfMaterialInserter The main cycle is managed by the model's -:cpp:class:`celeritas::DTMixMucfExecutor`. The muonic atom selection is handled -by detail classes, while the Interactors are reserved for sampling final states -of the outgoing secondaries. +:cpp:class:`celeritas::DTMixMucfExecutor`. The muonic atom and molecule +formation are handled by detail classes, while the Interactors are reserved for +sampling final states of the outgoing secondaries. .. note:: Only reactive channels are implemented. .. doxygenclass:: celeritas::detail::MuonicAtomSelector .. doxygenclass:: celeritas::detail::MuonicAtomSpinSelector +.. doxygenclass:: celeritas::detail::MuonicMoleculeSelector .. doxygenclass:: celeritas::detail::DDChannelSelector .. doxygenclass:: celeritas::detail::DTChannelSelector .. doxygenclass:: celeritas::detail::TTChannelSelector diff --git a/src/celeritas/mucf/executor/detail/DDChannelSelector.hh b/src/celeritas/mucf/executor/detail/DDChannelSelector.hh index 09d20698fe..c0f3cff4ba 100644 --- a/src/celeritas/mucf/executor/detail/DDChannelSelector.hh +++ b/src/celeritas/mucf/executor/detail/DDChannelSelector.hh @@ -6,7 +6,7 @@ //---------------------------------------------------------------------------// #pragma once -#include "corecel/random/distribution/GenerateCanonical.hh" +#include "corecel/random/distribution/BernoulliDistribution.hh" #include "celeritas/mucf/interactor/DDMucfInteractor.hh" namespace celeritas @@ -21,7 +21,7 @@ namespace detail * of the outcome of the fusion ending in the \f$ ^3\text{He} \f$ channels * versus the tritium channel. If the outcome is a \f$ ^3\text{He} \f$ channel, * a constant sticking fraction of 12.2% is used to define if sticking occurs - * [ \todo https://doi.org/10.1134/S106377961102002X ]. + * \citet{balin-mucf-2011, https://doi.org/10.1134/S106377961102002X} . */ class DDChannelSelector { @@ -90,23 +90,13 @@ DDChannelSelector::operator()(Engine& rng) { Channel result{Channel::size_}; - if (generate_canonical(rng) < he3_probability_) - { - // Select between the two 3He channels - if (generate_canonical(rng) > this->sticking_fraction()) - { - result = Channel::helium3_muon_neutron; - } - else - { - result = Channel::muonichelium3_neutron; - } - } - else - { - // Select tritium channel - result = Channel::tritium_muon_proton; - } + result = BernoulliDistribution(he3_probability_)(rng) + // Select between the two 3He channels + ? (BernoulliDistribution(this->sticking_fraction())(rng) + ? Channel::muonichelium3_neutron + : Channel::helium3_muon_neutron) + // Select tritium channel + : Channel::tritium_muon_proton; CELER_ENSURE(result < Channel::size_); return result; diff --git a/src/celeritas/mucf/executor/detail/DTChannelSelector.hh b/src/celeritas/mucf/executor/detail/DTChannelSelector.hh index 76246281f4..4fcf0da6a3 100644 --- a/src/celeritas/mucf/executor/detail/DTChannelSelector.hh +++ b/src/celeritas/mucf/executor/detail/DTChannelSelector.hh @@ -6,7 +6,7 @@ //---------------------------------------------------------------------------// #pragma once -#include "corecel/random/distribution/GenerateCanonical.hh" +#include "corecel/random/distribution/BernoulliDistribution.hh" #include "celeritas/mucf/interactor/DTMucfInteractor.hh" namespace celeritas @@ -18,8 +18,8 @@ namespace detail * Select final channel for muonic dt molecules. * * The selection is based on a constant sticking fraction from - * [ \todo https://arxiv.org/abs/2112.08399 ], in which ~0.8% of the time the - * muonic alpha channel is selected. + * \citet{kamimura-mucf-2023, https://doi.org/10.1103/PhysRevC.107.034607} + * in which ~0.8% of the time the muonic alpha channel is selected. */ class DTChannelSelector { @@ -58,9 +58,9 @@ DTChannelSelector::operator()(Engine& rng) { Channel result{Channel::size_}; - generate_canonical(rng) > this->sticking_fraction() - ? result = Channel::alpha_muon_neutron - : result = Channel::muonicalpha_neutron; + result = (BernoulliDistribution(this->sticking_fraction())(rng)) + ? Channel::muonicalpha_neutron + : Channel::alpha_muon_neutron; CELER_ENSURE(result < Channel::size_); return result; diff --git a/src/celeritas/mucf/executor/detail/TTChannelSelector.hh b/src/celeritas/mucf/executor/detail/TTChannelSelector.hh index 3c797a0374..36cc215972 100644 --- a/src/celeritas/mucf/executor/detail/TTChannelSelector.hh +++ b/src/celeritas/mucf/executor/detail/TTChannelSelector.hh @@ -6,7 +6,7 @@ //---------------------------------------------------------------------------// #pragma once -#include "corecel/random/distribution/GenerateCanonical.hh" +#include "corecel/random/distribution/BernoulliDistribution.hh" #include "celeritas/mucf/interactor/TTMucfInteractor.hh" namespace celeritas @@ -17,8 +17,8 @@ namespace detail /*! * Select final channel for muonic tt molecules. * - * The selection is based on a constant sticking fraction - * from [ \todo https://link.springer.com/article/10.1134/S1063776109020034 ], + * The selection is based on a constant sticking fraction from + * \citet{bogdanova-mucf-2009, https://doi.org/10.1134/S1063776109020034} , * in which ~14% of the time the muonic alpha channel is selected. */ class TTChannelSelector @@ -55,14 +55,9 @@ TTChannelSelector::operator()(Engine& rng) { Channel result{Channel::size_}; - if (generate_canonical(rng) > sticking_fraction()) - { - result = Channel::alpha_muon_neutron_neutron; - } - else - { - result = Channel::muonicalpha_neutron_neutron; - } + result = (BernoulliDistribution(this->sticking_fraction())(rng)) + ? Channel::muonicalpha_neutron_neutron + : Channel::alpha_muon_neutron_neutron; CELER_ENSURE(result < Channel::size_); return result; From a1c9621ed026fac4ef7f0715167844d872cacb7b Mon Sep 17 00:00:00 2001 From: Stefano Tognini Date: Fri, 20 Feb 2026 16:18:39 -0500 Subject: [PATCH 29/42] Add missing CELER_FUNCTION --- src/celeritas/mucf/executor/detail/MuonicMoleculeSelector.hh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/celeritas/mucf/executor/detail/MuonicMoleculeSelector.hh b/src/celeritas/mucf/executor/detail/MuonicMoleculeSelector.hh index aaea0244fe..9efbfa5928 100644 --- a/src/celeritas/mucf/executor/detail/MuonicMoleculeSelector.hh +++ b/src/celeritas/mucf/executor/detail/MuonicMoleculeSelector.hh @@ -44,7 +44,7 @@ class MuonicMoleculeSelector real_type cycle_time{std::numeric_limits::max()}; //! Check whether the data are assigned - explicit operator bool() const + explicit CELER_FUNCTION operator bool() const { return molecule < MucfMuonicMolecule::size_ && cycle_time > 0 && cycle_time < std::numeric_limits::max(); From d818b50c2cdc77ed66e76c93abac9ced0aac0cec Mon Sep 17 00:00:00 2001 From: Stefano Tognini Date: Fri, 20 Feb 2026 16:31:30 -0500 Subject: [PATCH 30/42] Use device-compatible numeric limits --- src/celeritas/inp/MucfPhysics.cc | 3 ++- .../mucf/executor/detail/MuonicMoleculeSelector.hh | 9 +++++---- 2 files changed, 7 insertions(+), 5 deletions(-) diff --git a/src/celeritas/inp/MucfPhysics.cc b/src/celeritas/inp/MucfPhysics.cc index 22cf2a963e..d5ede1fd98 100644 --- a/src/celeritas/inp/MucfPhysics.cc +++ b/src/celeritas/inp/MucfPhysics.cc @@ -503,7 +503,8 @@ MucfCycleRate hd_1_cycle_data() /*! * Muon-catalyzed fusion cycle rate data for \em tt fusion with F = 1/2 state. * - * The tt fusion cycle rate data is currently a constant value at 2.8e6 1/s. + * The tt fusion cycle rate data is currently a constant value at 2.8e6 1/s + * citet{bogdanova-mucf-2009, https://doi.org/10.1134/S1063776109020034} . */ MucfCycleRate tt_1_over_2_cycle_data() { diff --git a/src/celeritas/mucf/executor/detail/MuonicMoleculeSelector.hh b/src/celeritas/mucf/executor/detail/MuonicMoleculeSelector.hh index 9efbfa5928..c86352aa85 100644 --- a/src/celeritas/mucf/executor/detail/MuonicMoleculeSelector.hh +++ b/src/celeritas/mucf/executor/detail/MuonicMoleculeSelector.hh @@ -8,6 +8,7 @@ #include "corecel/cont/EnumArray.hh" #include "corecel/cont/Range.hh" +#include "corecel/math/NumericLimits.hh" #include "corecel/random/distribution/ExponentialDistribution.hh" #include "celeritas/Quantities.hh" #include "celeritas/Types.hh" @@ -41,13 +42,13 @@ class MuonicMoleculeSelector struct Result { MucfMuonicMolecule molecule{MucfMuonicMolecule::size_}; - real_type cycle_time{std::numeric_limits::max()}; + real_type cycle_time{numeric_limits::max()}; //! Check whether the data are assigned explicit CELER_FUNCTION operator bool() const { return molecule < MucfMuonicMolecule::size_ && cycle_time > 0 - && cycle_time < std::numeric_limits::max(); + && cycle_time < numeric_limits::max(); } }; @@ -131,7 +132,7 @@ MuonicMoleculeSelector::operator()(Engine& rng) using CycleTimeArray = EnumArray; CycleTimeArray sampled_times; - auto const inf = std::numeric_limits::max(); + auto const inf = numeric_limits::max(); if (atom_ == MucfMuonicAtom::deuterium) { // DD fusion is only triggered by a muonic deuterium @@ -180,7 +181,7 @@ CELER_FUNCTION real_type MuonicMoleculeSelector::sample_exp_time( // Return an infinite cycle time when the cached rate is zero return (rate > 0) ? ExponentialDistribution(rate)(rng) - : std::numeric_limits::max(); + : numeric_limits::max(); } //---------------------------------------------------------------------------// From 230e139eafff2cb6a0871e651086afd3cd9ea4a3 Mon Sep 17 00:00:00 2001 From: Stefano Tognini Date: Fri, 20 Feb 2026 18:06:17 -0500 Subject: [PATCH 31/42] Use double in ChannelSelector test; Fix windows cast error. --- test/celeritas/mucf/ChannelSelector.test.cc | 147 ++++++++++++-------- 1 file changed, 87 insertions(+), 60 deletions(-) diff --git a/test/celeritas/mucf/ChannelSelector.test.cc b/test/celeritas/mucf/ChannelSelector.test.cc index b8ac7caf63..d3d55b60ca 100644 --- a/test/celeritas/mucf/ChannelSelector.test.cc +++ b/test/celeritas/mucf/ChannelSelector.test.cc @@ -33,20 +33,20 @@ class ChannelSelectorTest : public Test Engine& rng() { return rng_; } // Sticking fraction between the two dd --> 3He channels - real_type const dd_sticking_fraction() { return 0.122; } + double const dd_sticking_fraction() { return 0.122; } // Sticking fraction for dt - real_type const dt_sticking_fraction() { return 0.00857; } + double const dt_sticking_fraction() { return 0.00857; } // Sticking fraction for tt - real_type const tt_sticking_fraction() { return 0.14; } + double const tt_sticking_fraction() { return 0.14; } // Calculate dd --> 3He channel probability from the branching ratio - real_type he3_probability(real_type branching_ratio) + double he3_probability(double branching_ratio) { return branching_ratio / (branching_ratio + 1); } // Calculate sigma for the statistical tests - real_type calc_sigma(real_type num_samples, real_type success_prob) + double calc_sigma(double num_samples, double success_prob) { return std::sqrt(num_samples * success_prob * (1 - success_prob)); } @@ -97,26 +97,33 @@ TEST_F(ChannelSelectorTest, dd_channel_low_temperature) EXPECT_EQ(num_samples, helium3_count + muonichelium3_count + tritium_count); - real_type expected_tritium_count = num_samples * (1 - he3_probability); - real_type expected_helium3_count = num_samples * he3_probability - * (1 - sticking_fraction); - real_type expected_muonichelium3_count = num_samples * he3_probability - * sticking_fraction; - // 3 sigma tolerance - real_type tolerance = 3 * this->calc_sigma(num_samples, he3_probability); - - EXPECT_NEAR(expected_tritium_count, tritium_count, tolerance); - EXPECT_NEAR(expected_helium3_count, helium3_count, tolerance); - EXPECT_NEAR(expected_muonichelium3_count, muonichelium3_count, tolerance); + if constexpr (CELERITAS_REAL_TYPE == CELERITAS_REAL_TYPE_DOUBLE) + { + double const num_samples_d = static_cast(num_samples); + + double const expected_tritium_count = num_samples_d + * (1 - he3_probability); + double const expected_helium3_count = num_samples_d * he3_probability + * (1 - sticking_fraction); + double const expected_muonichelium3_count + = num_samples_d * he3_probability * sticking_fraction; + double const tolerance + = 3 * this->calc_sigma(num_samples_d, he3_probability); + + EXPECT_NEAR(expected_tritium_count, tritium_count, tolerance); + EXPECT_NEAR(expected_helium3_count, helium3_count, tolerance); + EXPECT_NEAR( + expected_muonichelium3_count, muonichelium3_count, tolerance); + } } //---------------------------------------------------------------------------// TEST_F(ChannelSelectorTest, dd_channel_mid_temperature) { // DD fusion at 50 < T < 100 K: branching_ratio = 1.0088 * (T - 50) = 25.22 - real_type const temperature = 75.0; - real_type const branching_ratio = 1.0088 * (temperature - 50); - real_type const he3_probability = this->he3_probability(branching_ratio); + double const temperature = 75.0; + double const branching_ratio = 1.0088 * (temperature - 50); + double const he3_probability = this->he3_probability(branching_ratio); DDChannelSelector select_channel(temperature); @@ -140,20 +147,24 @@ TEST_F(ChannelSelectorTest, dd_channel_mid_temperature) EXPECT_EQ(num_samples, he3_total_count + tritium_count); - real_type const expected_he3_count = num_samples * he3_probability; - real_type const tolerance - = 3 * this->calc_sigma(num_samples, he3_probability); + if constexpr (CELERITAS_REAL_TYPE == CELERITAS_REAL_TYPE_DOUBLE) + { + double const num_samples_d = static_cast(num_samples); + double const expected_he3_count = num_samples_d * he3_probability; + double const tolerance + = 3 * this->calc_sigma(num_samples_d, he3_probability); - EXPECT_NEAR(expected_he3_count, he3_total_count, tolerance); + EXPECT_NEAR(expected_he3_count, he3_total_count, tolerance); + } } //---------------------------------------------------------------------------// TEST_F(ChannelSelectorTest, dd_channel_high_temperature) { // DD fusion at T < 300 K: branching_ratio = 1.44 - real_type const temperature = 300; - real_type const branching_ratio = 1.44; - real_type const he3_probability = this->he3_probability(branching_ratio); + double const temperature = 300; + double const branching_ratio = 1.44; + double const he3_probability = this->he3_probability(branching_ratio); DDChannelSelector select_channel(temperature); @@ -177,20 +188,24 @@ TEST_F(ChannelSelectorTest, dd_channel_high_temperature) EXPECT_EQ(num_samples, he3_total_count + tritium_count); - real_type const expected_he3_count = num_samples * he3_probability; - // 3 sigma tolerance - real_type const tolerance - = 3 * this->calc_sigma(num_samples, he3_probability); + if constexpr (CELERITAS_REAL_TYPE == CELERITAS_REAL_TYPE_DOUBLE) + { + double const num_samples_d = static_cast(num_samples); + double const expected_he3_count = num_samples_d * he3_probability; + // 3 sigma tolerance + double const tolerance + = 3 * this->calc_sigma(num_samples_d, he3_probability); - EXPECT_NEAR(expected_he3_count, he3_total_count, tolerance); + EXPECT_NEAR(expected_he3_count, he3_total_count, tolerance); + } } //---------------------------------------------------------------------------// TEST_F(ChannelSelectorTest, dd_sticking_fraction_within_he3) { // Test that when He3 channel is selected, sticking fraction is 12.2% - real_type const temperature = 300; - real_type const sticking_fraction = this->dd_sticking_fraction(); + double const temperature = 300; + double const sticking_fraction = this->dd_sticking_fraction(); DDChannelSelector select_channel(temperature); @@ -215,12 +230,16 @@ TEST_F(ChannelSelectorTest, dd_sticking_fraction_within_he3) size_type total_he3 = helium3_count + muonichelium3_count; EXPECT_GT(total_he3, 0); - real_type expected_muonichelium3 = total_he3 * sticking_fraction; - // 3 sigma tolerance - real_type const tolerance - = 3 * this->calc_sigma(total_he3, sticking_fraction); + if constexpr (CELERITAS_REAL_TYPE == CELERITAS_REAL_TYPE_DOUBLE) + { + double const total_he3_d = static_cast(total_he3); + double const expected_muonichelium3 = total_he3_d * sticking_fraction; + // 3 sigma tolerance + double const tolerance + = 3 * this->calc_sigma(total_he3_d, sticking_fraction); - EXPECT_NEAR(expected_muonichelium3, muonichelium3_count, tolerance); + EXPECT_NEAR(expected_muonichelium3, muonichelium3_count, tolerance); + } } //---------------------------------------------------------------------------// @@ -252,17 +271,21 @@ TEST_F(ChannelSelectorTest, dt_channel) EXPECT_EQ(num_samples, alpha_count + muonicalpha_count); - real_type const sticking_fraction = this->dt_sticking_fraction(); - real_type const expected_muonicalpha_count = num_samples - * sticking_fraction; - real_type const expected_alpha_count = num_samples - * (1 - sticking_fraction); - // 3 sigma tolerance - real_type const tolerance - = 3 * this->calc_sigma(num_samples, sticking_fraction); - - EXPECT_NEAR(expected_muonicalpha_count, muonicalpha_count, tolerance); - EXPECT_NEAR(expected_alpha_count, alpha_count, tolerance); + if constexpr (CELERITAS_REAL_TYPE == CELERITAS_REAL_TYPE_DOUBLE) + { + double const num_samples_d = static_cast(num_samples); + double const sticking_fraction = this->dt_sticking_fraction(); + double const expected_muonicalpha_count = num_samples_d + * sticking_fraction; + double const expected_alpha_count = num_samples_d + * (1 - sticking_fraction); + // 3 sigma tolerance + double const tolerance + = 3 * this->calc_sigma(num_samples_d, sticking_fraction); + + EXPECT_NEAR(expected_muonicalpha_count, muonicalpha_count, tolerance); + EXPECT_NEAR(expected_alpha_count, alpha_count, tolerance); + } } //---------------------------------------------------------------------------// @@ -295,17 +318,21 @@ TEST_F(ChannelSelectorTest, tt_channel) EXPECT_EQ(num_samples, alpha_count + muonicalpha_count); - real_type const sticking_fraction = this->tt_sticking_fraction(); - real_type const expected_muonicalpha_count = num_samples - * sticking_fraction; - real_type const expected_alpha_count = num_samples - * (1 - sticking_fraction); - - real_type const tolerance - = 3 * this->calc_sigma(num_samples, sticking_fraction); - - EXPECT_NEAR(expected_muonicalpha_count, muonicalpha_count, tolerance); - EXPECT_NEAR(expected_alpha_count, alpha_count, tolerance); + if constexpr (CELERITAS_REAL_TYPE == CELERITAS_REAL_TYPE_DOUBLE) + { + double const num_samples_d = static_cast(num_samples); + double const sticking_fraction = this->tt_sticking_fraction(); + double const expected_muonicalpha_count = num_samples_d + * sticking_fraction; + double const expected_alpha_count = num_samples_d + * (1 - sticking_fraction); + + double const tolerance + = 3 * this->calc_sigma(num_samples_d, sticking_fraction); + + EXPECT_NEAR(expected_muonicalpha_count, muonicalpha_count, tolerance); + EXPECT_NEAR(expected_alpha_count, alpha_count, tolerance); + } } //---------------------------------------------------------------------------// From 6873adc053d2e57b1bf020d3c3b352a1b3444db8 Mon Sep 17 00:00:00 2001 From: Stefano Tognini Date: Fri, 20 Feb 2026 19:39:07 -0500 Subject: [PATCH 32/42] More casts for windows --- test/celeritas/mucf/ChannelSelector.test.cc | 29 ++++++++++++++++----- 1 file changed, 22 insertions(+), 7 deletions(-) diff --git a/test/celeritas/mucf/ChannelSelector.test.cc b/test/celeritas/mucf/ChannelSelector.test.cc index d3d55b60ca..ff927b07e9 100644 --- a/test/celeritas/mucf/ChannelSelector.test.cc +++ b/test/celeritas/mucf/ChannelSelector.test.cc @@ -150,11 +150,13 @@ TEST_F(ChannelSelectorTest, dd_channel_mid_temperature) if constexpr (CELERITAS_REAL_TYPE == CELERITAS_REAL_TYPE_DOUBLE) { double const num_samples_d = static_cast(num_samples); + double const he3_total_count_d = static_cast(he3_total_count); + double const expected_he3_count = num_samples_d * he3_probability; double const tolerance = 3 * this->calc_sigma(num_samples_d, he3_probability); - EXPECT_NEAR(expected_he3_count, he3_total_count, tolerance); + EXPECT_NEAR(expected_he3_count, he3_total_count_d, tolerance); } } @@ -191,12 +193,14 @@ TEST_F(ChannelSelectorTest, dd_channel_high_temperature) if constexpr (CELERITAS_REAL_TYPE == CELERITAS_REAL_TYPE_DOUBLE) { double const num_samples_d = static_cast(num_samples); + double const he3_total_count_d = static_cast(he3_total_count); + double const expected_he3_count = num_samples_d * he3_probability; // 3 sigma tolerance double const tolerance = 3 * this->calc_sigma(num_samples_d, he3_probability); - EXPECT_NEAR(expected_he3_count, he3_total_count, tolerance); + EXPECT_NEAR(expected_he3_count, he3_total_count_d, tolerance); } } @@ -233,12 +237,15 @@ TEST_F(ChannelSelectorTest, dd_sticking_fraction_within_he3) if constexpr (CELERITAS_REAL_TYPE == CELERITAS_REAL_TYPE_DOUBLE) { double const total_he3_d = static_cast(total_he3); + double const muonichelium3_count_d + = static_cast(muonichelium3_count); + double const expected_muonichelium3 = total_he3_d * sticking_fraction; // 3 sigma tolerance double const tolerance = 3 * this->calc_sigma(total_he3_d, sticking_fraction); - EXPECT_NEAR(expected_muonichelium3, muonichelium3_count, tolerance); + EXPECT_NEAR(expected_muonichelium3, muonichelium3_count_d, tolerance); } } @@ -274,6 +281,10 @@ TEST_F(ChannelSelectorTest, dt_channel) if constexpr (CELERITAS_REAL_TYPE == CELERITAS_REAL_TYPE_DOUBLE) { double const num_samples_d = static_cast(num_samples); + double const alpha_count_d = static_cast(alpha_count); + double const muonicalpha_count_d + = static_cast(muonicalpha_count); + double const sticking_fraction = this->dt_sticking_fraction(); double const expected_muonicalpha_count = num_samples_d * sticking_fraction; @@ -283,8 +294,8 @@ TEST_F(ChannelSelectorTest, dt_channel) double const tolerance = 3 * this->calc_sigma(num_samples_d, sticking_fraction); - EXPECT_NEAR(expected_muonicalpha_count, muonicalpha_count, tolerance); - EXPECT_NEAR(expected_alpha_count, alpha_count, tolerance); + EXPECT_NEAR(expected_muonicalpha_count, muonicalpha_count_d, tolerance); + EXPECT_NEAR(expected_alpha_count, alpha_count_d, tolerance); } } @@ -321,6 +332,10 @@ TEST_F(ChannelSelectorTest, tt_channel) if constexpr (CELERITAS_REAL_TYPE == CELERITAS_REAL_TYPE_DOUBLE) { double const num_samples_d = static_cast(num_samples); + double const alpha_count_d = static_cast(alpha_count); + double const muonicalpha_count_d + = static_cast(muonicalpha_count); + double const sticking_fraction = this->tt_sticking_fraction(); double const expected_muonicalpha_count = num_samples_d * sticking_fraction; @@ -330,8 +345,8 @@ TEST_F(ChannelSelectorTest, tt_channel) double const tolerance = 3 * this->calc_sigma(num_samples_d, sticking_fraction); - EXPECT_NEAR(expected_muonicalpha_count, muonicalpha_count, tolerance); - EXPECT_NEAR(expected_alpha_count, alpha_count, tolerance); + EXPECT_NEAR(expected_muonicalpha_count, muonicalpha_count_d, tolerance); + EXPECT_NEAR(expected_alpha_count, alpha_count_d, tolerance); } } From 9cb451b8389738a9ecb361198d7f8bea9fd826a2 Mon Sep 17 00:00:00 2001 From: Stefano Tognini Date: Fri, 20 Feb 2026 19:41:47 -0500 Subject: [PATCH 33/42] Remove trailing `,` from bib file --- doc/_static/zotero.bib | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/doc/_static/zotero.bib b/doc/_static/zotero.bib index 562a5539c7..b634ddf078 100644 --- a/doc/_static/zotero.bib +++ b/doc/_static/zotero.bib @@ -3774,7 +3774,7 @@ @article{faifman-mucfformation-1996 pages = {179--189}, number = {1}, journaltitle = {Hyperfine Interactions}, - shortjournal = {Hyperfine Interact}, + shortjournal = {Hyperfine Interact} } @article{balin-mucf-2011, @@ -3790,7 +3790,7 @@ @article{balin-mucf-2011 pages = {185--214}, number = {2}, journaltitle = {Physics of Particles and Nuclei}, - shortjournal = {Phys. Part. Nuclei}, + shortjournal = {Phys. Part. Nuclei} } @article{bogdanova-mucf-2009, From 16a7847b7344b4ccbe2f2387c0b690f078aa5c0c Mon Sep 17 00:00:00 2001 From: Stefano Tognini Date: Sun, 22 Feb 2026 19:43:53 -0500 Subject: [PATCH 34/42] Add constexpr; Improve doc --- .../mucf/executor/detail/DDChannelSelector.hh | 5 ++++- .../mucf/executor/detail/DTChannelSelector.hh | 2 +- .../mucf/executor/detail/MuonicMoleculeSelector.hh | 11 ++++++----- .../mucf/executor/detail/TTChannelSelector.hh | 5 ++++- 4 files changed, 15 insertions(+), 8 deletions(-) diff --git a/src/celeritas/mucf/executor/detail/DDChannelSelector.hh b/src/celeritas/mucf/executor/detail/DDChannelSelector.hh index c0f3cff4ba..c15b5db4ef 100644 --- a/src/celeritas/mucf/executor/detail/DDChannelSelector.hh +++ b/src/celeritas/mucf/executor/detail/DDChannelSelector.hh @@ -42,7 +42,10 @@ class DDChannelSelector real_type he3_probability_; // Constant sticking fraction between the two 3He channels - inline CELER_FUNCTION real_type sticking_fraction() const { return 0.122; } + inline CELER_FUNCTION static constexpr real_type sticking_fraction() + { + return 0.122; + } }; //---------------------------------------------------------------------------// diff --git a/src/celeritas/mucf/executor/detail/DTChannelSelector.hh b/src/celeritas/mucf/executor/detail/DTChannelSelector.hh index 4fcf0da6a3..e2c8552ca1 100644 --- a/src/celeritas/mucf/executor/detail/DTChannelSelector.hh +++ b/src/celeritas/mucf/executor/detail/DTChannelSelector.hh @@ -38,7 +38,7 @@ class DTChannelSelector private: // Constant sticking fraction of dt fusion - inline CELER_FUNCTION real_type sticking_fraction() const + inline CELER_FUNCTION static constexpr real_type sticking_fraction() { return 0.00857; } diff --git a/src/celeritas/mucf/executor/detail/MuonicMoleculeSelector.hh b/src/celeritas/mucf/executor/detail/MuonicMoleculeSelector.hh index c86352aa85..35b61c30d7 100644 --- a/src/celeritas/mucf/executor/detail/MuonicMoleculeSelector.hh +++ b/src/celeritas/mucf/executor/detail/MuonicMoleculeSelector.hh @@ -23,11 +23,12 @@ namespace detail * Select a muonic molecule and its final cycle time. * * The muonic molecule is selected by sampling the shortest fusion cycle time - * (\f$ T = -\ln(r) \times \tau_\text{cycle} \f$). + * \f[ + * T = -\ln(r) \times \tau_\text{cycle}, + * \f] * - * \note The cycle rates cached in the model data are calculated based on the - * material definition, and thus are propagated to the final cycle time sampled - * here and used in the fusion process. + * where \f$ \tau_\text{cycle} \f$ is the cycle time a given molecule + spin + * calculated from material data. */ class MuonicMoleculeSelector { @@ -141,7 +142,7 @@ MuonicMoleculeSelector::operator()(Engine& rng) sampled_times[MMM::deuterium_tritium] = inf; sampled_times[MMM::tritium_tritium] = inf; } - if (atom_ == MucfMuonicAtom::tritium) + else if (atom_ == MucfMuonicAtom::tritium) { // DT and TT fusions are triggered by a muonic tritium sampled_times[MMM::deuterium_deuterium] = inf; diff --git a/src/celeritas/mucf/executor/detail/TTChannelSelector.hh b/src/celeritas/mucf/executor/detail/TTChannelSelector.hh index 36cc215972..7916c1f797 100644 --- a/src/celeritas/mucf/executor/detail/TTChannelSelector.hh +++ b/src/celeritas/mucf/executor/detail/TTChannelSelector.hh @@ -38,7 +38,10 @@ class TTChannelSelector private: // Constant sticking fraction of tt fusion - inline CELER_FUNCTION real_type sticking_fraction() const { return 0.14; } + inline CELER_FUNCTION real_type static constexpr sticking_fraction() + { + return 0.14; + } }; //---------------------------------------------------------------------------// From 3e15238665ef9a56806fc286f175d1c7fa315d52 Mon Sep 17 00:00:00 2001 From: Stefano Tognini Date: Tue, 24 Feb 2026 09:15:19 -0600 Subject: [PATCH 35/42] Add todo notes; fix build error --- .../mucf/executor/detail/DDChannelSelector.hh | 2 ++ .../mucf/executor/detail/DTChannelSelector.hh | 2 ++ .../mucf/executor/detail/TTChannelSelector.hh | 2 ++ test/celeritas/mucf/ChannelSelector.test.cc | 14 +++++++------- 4 files changed, 13 insertions(+), 7 deletions(-) diff --git a/src/celeritas/mucf/executor/detail/DDChannelSelector.hh b/src/celeritas/mucf/executor/detail/DDChannelSelector.hh index c15b5db4ef..11b896c23c 100644 --- a/src/celeritas/mucf/executor/detail/DDChannelSelector.hh +++ b/src/celeritas/mucf/executor/detail/DDChannelSelector.hh @@ -22,6 +22,8 @@ namespace detail * versus the tritium channel. If the outcome is a \f$ ^3\text{He} \f$ channel, * a constant sticking fraction of 12.2% is used to define if sticking occurs * \citet{balin-mucf-2011, https://doi.org/10.1134/S106377961102002X} . + * + * \todo Update I/O with user-defined sticking fractions and branching ratios. */ class DDChannelSelector { diff --git a/src/celeritas/mucf/executor/detail/DTChannelSelector.hh b/src/celeritas/mucf/executor/detail/DTChannelSelector.hh index e2c8552ca1..b890d931f7 100644 --- a/src/celeritas/mucf/executor/detail/DTChannelSelector.hh +++ b/src/celeritas/mucf/executor/detail/DTChannelSelector.hh @@ -20,6 +20,8 @@ namespace detail * The selection is based on a constant sticking fraction from * \citet{kamimura-mucf-2023, https://doi.org/10.1103/PhysRevC.107.034607} * in which ~0.8% of the time the muonic alpha channel is selected. + * + * \todo Update I/O with user-defined sticking fractions. */ class DTChannelSelector { diff --git a/src/celeritas/mucf/executor/detail/TTChannelSelector.hh b/src/celeritas/mucf/executor/detail/TTChannelSelector.hh index 7916c1f797..3716d073b4 100644 --- a/src/celeritas/mucf/executor/detail/TTChannelSelector.hh +++ b/src/celeritas/mucf/executor/detail/TTChannelSelector.hh @@ -20,6 +20,8 @@ namespace detail * The selection is based on a constant sticking fraction from * \citet{bogdanova-mucf-2009, https://doi.org/10.1134/S1063776109020034} , * in which ~14% of the time the muonic alpha channel is selected. + * + * \todo Update I/O with user-defined sticking fractions. */ class TTChannelSelector { diff --git a/test/celeritas/mucf/ChannelSelector.test.cc b/test/celeritas/mucf/ChannelSelector.test.cc index ff927b07e9..442c79aadf 100644 --- a/test/celeritas/mucf/ChannelSelector.test.cc +++ b/test/celeritas/mucf/ChannelSelector.test.cc @@ -33,11 +33,11 @@ class ChannelSelectorTest : public Test Engine& rng() { return rng_; } // Sticking fraction between the two dd --> 3He channels - double const dd_sticking_fraction() { return 0.122; } + double dd_sticking_fraction() { return 0.122; } // Sticking fraction for dt - double const dt_sticking_fraction() { return 0.00857; } + double dt_sticking_fraction() { return 0.00857; } // Sticking fraction for tt - double const tt_sticking_fraction() { return 0.14; } + double tt_sticking_fraction() { return 0.14; } // Calculate dd --> 3He channel probability from the branching ratio double he3_probability(double branching_ratio) @@ -62,10 +62,10 @@ class ChannelSelectorTest : public Test TEST_F(ChannelSelectorTest, dd_channel_low_temperature) { // dd fusion at T < 50 K: branching_ratio = 1 - real_type const temperature = 30.0; - real_type const branching_ratio = 1.0; - real_type const he3_probability = this->he3_probability(branching_ratio); - real_type const sticking_fraction = this->dd_sticking_fraction(); + double const temperature = 30.0; + double const branching_ratio = 1.0; + double const he3_probability = this->he3_probability(branching_ratio); + double const sticking_fraction = this->dd_sticking_fraction(); DDChannelSelector select_channel(temperature); From 272bda42bf8f301bc15ed124a05ad2a799d4df88 Mon Sep 17 00:00:00 2001 From: Stefano Tognini Date: Wed, 25 Feb 2026 11:43:40 -0600 Subject: [PATCH 36/42] Update zotero --- doc/_static/zotero.bib | 37 +++++++++++++++++++------------------ 1 file changed, 19 insertions(+), 18 deletions(-) diff --git a/doc/_static/zotero.bib b/doc/_static/zotero.bib index b634ddf078..d66545ddc0 100644 --- a/doc/_static/zotero.bib +++ b/doc/_static/zotero.bib @@ -3778,33 +3778,34 @@ @article{faifman-mucfformation-1996 } @article{balin-mucf-2011, - title = {{High precision study of muon catalyzed fusion in D2 and {HD} gas}}, + title = {High Precision Study of Muon Catalyzed Fusion in {{D2}} and {{HD}} Gas}, author = {Balin, D. V. and Ganzha, V. A. and Kozlov, S. M. and Maev, E. M. and Petrov, G. E. and Soroka, M. A. and Schapkin, G. N. and Semenchuk, G. G. and Trofimov, V. A. and Vasiliev, A. A. and Vorobyov, A. A. and Voropaev, N. I. and Petitjean, C. and Gartner, B. and Lauss, B. and Marton, J. and Zmeskal, J. and Case, T. and Crowe, K. M. and Kammel, P. and Hartmann, F. J. and Faifman, M. P.}, + year = 2011, + month = mar, + journal = {Physics of Particles and Nuclei}, volume = {42}, - rights = {http://www.springer.com/tdm}, + number = {2}, + pages = {185--214}, issn = {1063-7796, 1531-8559}, - url = {http://link.springer.com/10.1134/S106377961102002X}, doi = {10.1134/S106377961102002X}, - year = 2011, - month = apr, - pages = {185--214}, - number = {2}, - journaltitle = {Physics of Particles and Nuclei}, - shortjournal = {Phys. Part. Nuclei} + urldate = {2026-02-20}, + copyright = {http://www.springer.com/tdm}, + langid = {english} } + @article{bogdanova-mucf-2009, - title = {Experimental investigation of muon-catalyzed t + t fusion}, - author = {Bogdanova, L. N. and Bom, V. R. and Demin, A. M. and Demin, D. L. and Van Eijk, C. W. E. and Filchagin, S. V. and Filchenkov, V. V. and Grafov, N. N. and Grishechkin, S. K. and Gritsaj, K. I. and Konin, A. D. and Kuryakin, A. V. and Medved’, S. V. and Musyaev, R. K. and Rudenko, A. I. and Tumkin, D. P. and Vinogradov, Yu. I. and Yukhimchuk, A. A. and Yukhimchuk, S. A. and Zinov, V. G. and Zlatoustovskii, S. V.}, - volume = {108}, + title = {{Experimental Investigation of Muon-Catalyzed t + t Fusion}}, + author = {Bogdanova, L. N. and Bom, V. R. and Demin, A. M. and Demin, D. L. and Van Eijk, C. W. E. and Filchagin, S. V. and Filchenkov, V. V. and Grafov, N. N. and Grishechkin, S. K. and Gritsaj, K. I. and Konin, A. D. and Kuryakin, A. V. and Medved', S. V. and Musyaev, R. K. and Rudenko, A. I. and Tumkin, D. P. and Vinogradov, {\relax Yu}. I. and Yukhimchuk, A. A. and Yukhimchuk, S. A. and Zinov, V. G. and Zlatoustovskii, S. V.}, year = 2009, month = feb, - rights = {http://www.springer.com/tdm}, + journal = {Journal of Experimental and Theoretical Physics}, + volume = {108}, + number = {2}, + pages = {216--225}, issn = {1063-7761, 1090-6509}, - url = {http://link.springer.com/10.1134/S1063776109020034}, doi = {10.1134/S1063776109020034}, - pages = {216--225}, - number = {2}, - journaltitle = {Journal of Experimental and Theoretical Physics}, - shortjournal = {J. Exp. Theor. Phys.} + urldate = {2026-02-20}, + copyright = {http://www.springer.com/tdm}, + langid = {english} } From 58767d86709cdc2b07e359b660baccf8fd308b88 Mon Sep 17 00:00:00 2001 From: Stefano Tognini Date: Wed, 25 Feb 2026 13:43:15 -0600 Subject: [PATCH 37/42] Remove trailing line --- doc/_static/zotero.bib | 1 - 1 file changed, 1 deletion(-) diff --git a/doc/_static/zotero.bib b/doc/_static/zotero.bib index d66545ddc0..bc5af07237 100644 --- a/doc/_static/zotero.bib +++ b/doc/_static/zotero.bib @@ -3793,7 +3793,6 @@ @article{balin-mucf-2011 langid = {english} } - @article{bogdanova-mucf-2009, title = {{Experimental Investigation of Muon-Catalyzed t + t Fusion}}, author = {Bogdanova, L. N. and Bom, V. R. and Demin, A. M. and Demin, D. L. and Van Eijk, C. W. E. and Filchagin, S. V. and Filchenkov, V. V. and Grafov, N. N. and Grishechkin, S. K. and Gritsaj, K. I. and Konin, A. D. and Kuryakin, A. V. and Medved', S. V. and Musyaev, R. K. and Rudenko, A. I. and Tumkin, D. P. and Vinogradov, {\relax Yu}. I. and Yukhimchuk, A. A. and Yukhimchuk, S. A. and Zinov, V. G. and Zlatoustovskii, S. V.}, From 8997e260ab07cc7c5ba9420c0ac346dddd54f6aa Mon Sep 17 00:00:00 2001 From: Stefano Tognini Date: Fri, 27 Feb 2026 10:20:18 -0500 Subject: [PATCH 38/42] Add [[maybe_unused]] to test --- test/celeritas/mucf/MuonicMoleculeSelector.test.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/celeritas/mucf/MuonicMoleculeSelector.test.cc b/test/celeritas/mucf/MuonicMoleculeSelector.test.cc index 24319aad2d..e03e4c2c42 100644 --- a/test/celeritas/mucf/MuonicMoleculeSelector.test.cc +++ b/test/celeritas/mucf/MuonicMoleculeSelector.test.cc @@ -134,7 +134,7 @@ TEST_F(MuonicMoleculeSelectorTest, model_data) auto const& rates = data_.cycle_rates[mucfmatid_]; MoleculeCountArray molecule_counts; - for (auto i : range(num_samples)) + for ([[maybe_unused]] auto i : range(num_samples)) { // Form muonic atom auto muonic_atom = MuonicAtomSelector( From aa6833d4ac13423fbeef298412b1e95d939d7979 Mon Sep 17 00:00:00 2001 From: Stefano Tognini Date: Fri, 27 Feb 2026 11:04:04 -0500 Subject: [PATCH 39/42] Add constexpr for reactive states --- src/celeritas/mucf/Types.hh | 29 +++++++++++++++++++ .../executor/detail/MuonicMoleculeSelector.hh | 11 +++---- 2 files changed, 35 insertions(+), 5 deletions(-) diff --git a/src/celeritas/mucf/Types.hh b/src/celeritas/mucf/Types.hh index c986fcf900..cbfe46d169 100644 --- a/src/celeritas/mucf/Types.hh +++ b/src/celeritas/mucf/Types.hh @@ -7,6 +7,7 @@ #pragma once #include "corecel/OpaqueId.hh" +#include "celeritas/Quantities.hh" namespace celeritas { @@ -57,5 +58,33 @@ enum class MucfMuonicMolecule //! Opaque index of a muCF material component using MucfMatId = OpaqueId; +//---------------------------------------------------------------------------// +// REACTIVE SPIN STATES +//---------------------------------------------------------------------------// + +//! State F = 1/2 (dd and tt) +inline constexpr units::HalfSpinInt spin_one_half() +{ + return units::HalfSpinInt{1}; +} + +//! State F = 3/2 (dd) +inline constexpr units::HalfSpinInt spin_three_halves() +{ + return units::HalfSpinInt{3}; +} + +//! State F = 0 (dt) +inline constexpr units::HalfSpinInt spin_zero() +{ + return units::HalfSpinInt{0}; +} + +//! State F = 1 (dt) +inline constexpr units::HalfSpinInt spin_one() +{ + return units::HalfSpinInt{2}; +} + //---------------------------------------------------------------------------// } // namespace celeritas diff --git a/src/celeritas/mucf/executor/detail/MuonicMoleculeSelector.hh b/src/celeritas/mucf/executor/detail/MuonicMoleculeSelector.hh index 35b61c30d7..a79e457d6e 100644 --- a/src/celeritas/mucf/executor/detail/MuonicMoleculeSelector.hh +++ b/src/celeritas/mucf/executor/detail/MuonicMoleculeSelector.hh @@ -101,19 +101,20 @@ MuonicMoleculeSelector::MuonicMoleculeSelector(MucfMuonicAtom atom, { case MucfMuonicAtom::deuterium: { // Related to DD states F = 1/2 and F = 3/2 - CELER_EXPECT(spin == HalfSpinInt{1} || spin == HalfSpinInt{3}); + CELER_EXPECT(spin == spin_one_half() + || spin == spin_three_halves()); // F = 1/2 and F = 3/2 correspond to indices 0 and 1, respectively - cycle_rate_index_ = (spin == HalfSpinInt{1}) ? 0 : 1; + cycle_rate_index_ = (spin == spin_one_half()) ? 0 : 1; break; } case MucfMuonicAtom::tritium: { // Related to DT states F = 0 and F = 1, and TT state F = 1/2 - CELER_EXPECT(spin == HalfSpinInt{0} || spin == HalfSpinInt{2} - || spin == HalfSpinInt{1}); + CELER_EXPECT(spin == spin_zero() || spin == spin_one() + || spin == spin_one_half()); // DT: F = 0 and F = 1 correspond to indices 0 and 1, respectively // TT: F = 1/2 corresponds to index 0 cycle_rate_index_ - = (spin == HalfSpinInt{0} || spin == HalfSpinInt{1}) ? 0 : 1; + = (spin == spin_zero() || spin == spin_one_half()) ? 0 : 1; break; } default: From 2070ad4f8a53b61e92a2fad78bd36c73939a5e77 Mon Sep 17 00:00:00 2001 From: Stefano Tognini Date: Fri, 27 Feb 2026 11:11:26 -0500 Subject: [PATCH 40/42] Fix if/else statement --- src/celeritas/mucf/executor/detail/MuonicMoleculeSelector.hh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/celeritas/mucf/executor/detail/MuonicMoleculeSelector.hh b/src/celeritas/mucf/executor/detail/MuonicMoleculeSelector.hh index a79e457d6e..796dcf2624 100644 --- a/src/celeritas/mucf/executor/detail/MuonicMoleculeSelector.hh +++ b/src/celeritas/mucf/executor/detail/MuonicMoleculeSelector.hh @@ -143,7 +143,7 @@ MuonicMoleculeSelector::operator()(Engine& rng) sampled_times[MMM::deuterium_tritium] = inf; sampled_times[MMM::tritium_tritium] = inf; } - else if (atom_ == MucfMuonicAtom::tritium) + else { // DT and TT fusions are triggered by a muonic tritium sampled_times[MMM::deuterium_deuterium] = inf; From ac85854efa8ec8e4035c4b66db6d46a37c33124e Mon Sep 17 00:00:00 2001 From: Stefano Tognini Date: Fri, 27 Feb 2026 12:38:34 -0500 Subject: [PATCH 41/42] Make spin states variables --- src/celeritas/mucf/Types.hh | 20 ++++--------------- .../executor/detail/MuonicMoleculeSelector.hh | 14 ++++++------- 2 files changed, 10 insertions(+), 24 deletions(-) diff --git a/src/celeritas/mucf/Types.hh b/src/celeritas/mucf/Types.hh index cbfe46d169..f23f7e3b80 100644 --- a/src/celeritas/mucf/Types.hh +++ b/src/celeritas/mucf/Types.hh @@ -63,28 +63,16 @@ using MucfMatId = OpaqueId; //---------------------------------------------------------------------------// //! State F = 1/2 (dd and tt) -inline constexpr units::HalfSpinInt spin_one_half() -{ - return units::HalfSpinInt{1}; -} +inline constexpr units::HalfSpinInt spin_one_half{1}; //! State F = 3/2 (dd) -inline constexpr units::HalfSpinInt spin_three_halves() -{ - return units::HalfSpinInt{3}; -} +inline constexpr units::HalfSpinInt spin_three_halves{3}; //! State F = 0 (dt) -inline constexpr units::HalfSpinInt spin_zero() -{ - return units::HalfSpinInt{0}; -} +inline constexpr units::HalfSpinInt spin_zero{0}; //! State F = 1 (dt) -inline constexpr units::HalfSpinInt spin_one() -{ - return units::HalfSpinInt{2}; -} +inline constexpr units::HalfSpinInt spin_one{2}; //---------------------------------------------------------------------------// } // namespace celeritas diff --git a/src/celeritas/mucf/executor/detail/MuonicMoleculeSelector.hh b/src/celeritas/mucf/executor/detail/MuonicMoleculeSelector.hh index 796dcf2624..c5e70045fc 100644 --- a/src/celeritas/mucf/executor/detail/MuonicMoleculeSelector.hh +++ b/src/celeritas/mucf/executor/detail/MuonicMoleculeSelector.hh @@ -100,21 +100,19 @@ MuonicMoleculeSelector::MuonicMoleculeSelector(MucfMuonicAtom atom, switch (atom_) { case MucfMuonicAtom::deuterium: { - // Related to DD states F = 1/2 and F = 3/2 - CELER_EXPECT(spin == spin_one_half() - || spin == spin_three_halves()); + CELER_EXPECT(spin == spin_one_half || spin == spin_three_halves); + // DD fusion // F = 1/2 and F = 3/2 correspond to indices 0 and 1, respectively - cycle_rate_index_ = (spin == spin_one_half()) ? 0 : 1; + cycle_rate_index_ = (spin == spin_one_half) ? 0 : 1; break; } case MucfMuonicAtom::tritium: { - // Related to DT states F = 0 and F = 1, and TT state F = 1/2 - CELER_EXPECT(spin == spin_zero() || spin == spin_one() - || spin == spin_one_half()); + CELER_EXPECT(spin == spin_zero || spin == spin_one + || spin == spin_one_half); // DT: F = 0 and F = 1 correspond to indices 0 and 1, respectively // TT: F = 1/2 corresponds to index 0 cycle_rate_index_ - = (spin == spin_zero() || spin == spin_one_half()) ? 0 : 1; + = (spin == spin_zero || spin == spin_one_half) ? 0 : 1; break; } default: From 0bee8dc9e6f6970ab5115afb1c602c796e58d619 Mon Sep 17 00:00:00 2001 From: Stefano Tognini Date: Thu, 5 Mar 2026 09:42:13 -0500 Subject: [PATCH 42/42] Enable tests only when double precision is available --- test/celeritas/mucf/MuonicMoleculeSelector.test.cc | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/celeritas/mucf/MuonicMoleculeSelector.test.cc b/test/celeritas/mucf/MuonicMoleculeSelector.test.cc index e03e4c2c42..256b872ffa 100644 --- a/test/celeritas/mucf/MuonicMoleculeSelector.test.cc +++ b/test/celeritas/mucf/MuonicMoleculeSelector.test.cc @@ -57,7 +57,7 @@ class MuonicMoleculeSelectorTest //---------------------------------------------------------------------------// // TESTS //---------------------------------------------------------------------------// -TEST_F(MuonicMoleculeSelectorTest, dd_spin_selection) +TEST_F(MuonicMoleculeSelectorTest, TEST_IF_CELERITAS_DOUBLE(dd_spin_selection)) { CycleRatesArray rates; @@ -83,7 +83,7 @@ TEST_F(MuonicMoleculeSelectorTest, dd_spin_selection) } //---------------------------------------------------------------------------// -TEST_F(MuonicMoleculeSelectorTest, dt_spin_selection) +TEST_F(MuonicMoleculeSelectorTest, TEST_IF_CELERITAS_DOUBLE(dt_spin_selection)) { CycleRatesArray rates; @@ -109,7 +109,7 @@ TEST_F(MuonicMoleculeSelectorTest, dt_spin_selection) } //---------------------------------------------------------------------------// -TEST_F(MuonicMoleculeSelectorTest, tt_spin_selection) +TEST_F(MuonicMoleculeSelectorTest, TEST_IF_CELERITAS_DOUBLE(tt_spin_selection)) { CycleRatesArray rates;