Add MuCF muonic molecule selector and fusion channel selectors

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.},
+  year = 2011,
+  month = mar,
+  journal = {Physics of Particles and Nuclei},
+  volume = {42},
+  number = {2},
+  pages = {185--214},
+  issn = {1063-7796, 1531-8559},
+  doi = {10.1134/S106377961102002X},
+  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, {\relax Yu}. I. and Yukhimchuk, A. A. and Yukhimchuk, S. A. and Zinov, V. G. and Zlatoustovskii, S. V.},
+  year = 2009,
+  month = feb,
+  journal = {Journal of Experimental and Theoretical Physics},
+  volume = {108},
+  number = {2},
+  pages = {216--225},
+  issn = {1063-7761, 1090-6509},
+  doi = {10.1134/S1063776109020034},
+  urldate = {2026-02-20},
+  copyright = {http://www.springer.com/tdm},
+  langid = {english}
 }

doc/implementation/mucf-physics.rst

@@ -133,14 +133,18 @@ 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
 .. doxygenclass:: celeritas::DDMucfInteractor
 .. doxygenclass:: celeritas::DTMucfInteractor
 .. doxygenclass:: celeritas::TTMucfInteractor

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()
 {

src/celeritas/mucf/Types.hh

@@ -7,6 +7,7 @@
 #pragma once
 
 #include "corecel/OpaqueId.hh"
+#include "celeritas/Quantities.hh"
 
 namespace celeritas
 {
@@ -55,7 +56,23 @@ enum class MucfMuonicMolecule
 //---------------------------------------------------------------------------//
 
 //! Opaque index of a muCF material component
-using MuCfMatId = OpaqueId<struct MuCfMaterialComponent_>;
+using MucfMatId = OpaqueId<struct MucfMaterialComponent_>;
+
+//---------------------------------------------------------------------------//
+// REACTIVE SPIN STATES
+//---------------------------------------------------------------------------//
+
+//! State F = 1/2 (dd and tt)
+inline constexpr units::HalfSpinInt spin_one_half{1};
+
+//! State F = 3/2 (dd)
+inline constexpr units::HalfSpinInt spin_three_halves{3};
+
+//! State F = 0 (dt)
+inline constexpr units::HalfSpinInt spin_zero{0};
+
+//! State F = 1 (dt)
+inline constexpr units::HalfSpinInt spin_one{2};
 
 //---------------------------------------------------------------------------//
 }  // namespace celeritas

src/celeritas/mucf/data/DTMixMucfData.hh

@@ -106,9 +106,9 @@ struct DTMixMucfData
     template<class T>
     using Items = Collection<T, W, M>;
     template<class T>
-    using MaterialItems = Collection<T, W, M, MuCfMatId>;
+    using MaterialItems = Collection<T, W, M, MucfMatId>;
     using GridRecord = NonuniformGridRecord;
-    using CycleTimesArray = EnumArray<MucfMuonicMolecule, Array<real_type, 2>>;
+    using CycleRatesArray = EnumArray<MucfMuonicMolecule, Array<real_type, 2>>;
     using MaterialFractionsArray = EnumArray<MucfIsotope, real_type>;
 
     //! Particles
@@ -122,12 +122,12 @@ struct DTMixMucfData
 
     //!@{
     //! Material-dependent data calculated at model construction
-    //! \c PhysMatId indexed by \c MuCfMatId
+    //! \c PhysMatId indexed by \c MucfMatId
     MaterialItems<PhysMatId> mucfmatid_to_matid;
     //! Isotopic fractions per material: [mat_comp_id][isotope]
     MaterialItems<MaterialFractionsArray> isotopic_fractions;
-    //! Cycle times per material: [mat_comp_id][muonic_molecule][spin_index]
-    MaterialItems<CycleTimesArray> cycle_times;  //!< In [s]
+    //! Cycle rates per material: [mat_comp_id][muonic_molecule][spin_index]
+    MaterialItems<CycleRatesArray> 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;
     }

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
@@ -44,39 +43,36 @@ 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();
-
-    // Muon decay may compete against other "actions" in this executor
-    real_type const decay_len{};  //! \todo Set muon decay interaction length
+    //! \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 {
+    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);
 
-    // 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
@@ -85,39 +81,23 @@ 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::MuonicMoleculeSelector(
+        muonic_atom, atom_spin, 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)
     {
         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;
         }