Add MuCF interactor for muonic dd molecules

src/celeritas/mucf/data/DTMixMucfData.hh

@@ -30,8 +30,9 @@ struct MucfParticleIds
 
     //!@{
     //! Elementary particles and nuclei
-    ParticleId neutron;
     ParticleId proton;
+    ParticleId triton;
+    ParticleId neutron;
     ParticleId alpha;
     ParticleId he3;
     //!@}
@@ -48,8 +49,9 @@ struct MucfParticleIds
     //! Check whether all particles are assigned
     CELER_FUNCTION explicit operator bool() const
     {
-        return mu_minus && neutron && proton && alpha && he3 && muonic_hydrogen
-               && muonic_alpha && muonic_triton && muonic_he3;
+        return mu_minus && proton && triton && neutron && alpha && he3
+               && muonic_hydrogen && muonic_alpha && muonic_triton
+               && muonic_he3;
     }
 };
 
@@ -64,8 +66,9 @@ struct MucfParticleMasses
 
     //!@{
     //! Elementary particles and nuclei
-    units::MevMass neutron;
     units::MevMass proton;
+    units::MevMass triton;
+    units::MevMass neutron;
     units::MevMass alpha;
     units::MevMass he3;
     //!@}
@@ -82,9 +85,10 @@ struct MucfParticleMasses
     //! Check whether all data are assigned
     CELER_FUNCTION explicit operator bool() const
     {
-        return mu_minus > zero_quantity() && neutron > zero_quantity()
-               && proton > zero_quantity() && alpha > zero_quantity()
-               && he3 > zero_quantity() && muonic_hydrogen > zero_quantity()
+        return mu_minus > zero_quantity() && proton > zero_quantity()
+               && triton > zero_quantity() && neutron > zero_quantity()
+               && alpha > zero_quantity() && he3 > zero_quantity()
+               && muonic_hydrogen > zero_quantity()
                && muonic_alpha > zero_quantity()
                && muonic_triton > zero_quantity()
                && muonic_he3 > zero_quantity();

src/celeritas/mucf/interactor/DDMucfInteractor.hh

@@ -8,13 +8,12 @@
 
 #include "corecel/Macros.hh"
 #include "corecel/data/StackAllocator.hh"
-#include "celeritas/mat/ElementView.hh"
-#include "celeritas/mat/MaterialView.hh"
 #include "celeritas/mucf/data/DTMixMucfData.hh"
 #include "celeritas/phys/Interaction.hh"
-#include "celeritas/phys/ParticleTrackView.hh"
 #include "celeritas/phys/Secondary.hh"
 
+#include "detail/MucfInteractorUtils.hh"
+
 namespace celeritas
 {
 //---------------------------------------------------------------------------//
@@ -24,17 +23,25 @@ namespace celeritas
  * Fusion channels:
  * - \f$ ^3\text{He} + \mu + n \f$
  * - \f$ (^3\text{He})_\mu + n \f$
- * - \f$ ^3\text{H} + \mu + p \f$
+ * - \f$ t + \mu + p \f$
+ *
+ * The mass ratios between \f$ ^3\text{He} \f$ and the neutron, and between
+ * tritium and the proton, are both about 3:1. This leads to the neutron and
+ * proton kinetic energies being 3/4 of the total kinetic energy available in
+ * their respective channels.
+ *
+ * \note This interactor has a similar implementation as \c DTMucfInteractor ,
+ * where energy is not fully conserved. See its documentation for details.
  */
 class DDMucfInteractor
 {
   public:
-    //! \todo Implement muonichydrogen3_proton (\f$ (^3\text{H})_\mu + p \f$)
+    //! \todo Implement muonictriton_proton (\f$ (t)_\mu + p \f$)
     enum class Channel
     {
         helium3_muon_neutron,  //!< \f$ ^3\text{He} + \mu + n \f$
         muonichelium3_neutron,  //!< \f$ (^3\text{He})_\mu + n \f$
-        hydrogen3_muon_proton,  //!< \f$ ^3\text{H} + \mu + p \f$
+        tritium_muon_proton,  //!< \f$ t + \mu + p \f$
         size_
     };
 
@@ -52,20 +59,39 @@ class DDMucfInteractor
     // Shared constant physics properties
     NativeCRef<DTMixMucfData> const& data_;
     // Selected fusion channel
-    Channel channel_{Channel::size_};
+    Channel channel_;
     // Allocate space for secondary particles
     StackAllocator<Secondary>& allocate_;
     // Number of secondaries per channel
     EnumArray<Channel, size_type> num_secondaries_{
         3,  // helium3_muon_neutron
         2,  // muonichelium3_neutron
-        3  // hydrogen3_muon_proton
+        3  // tritium_muon_proton
     };
 
-    // Sample Interaction secondaries
-    template<class Engine>
-    inline CELER_FUNCTION Span<Secondary>
-    sample_secondaries(Secondary* secondaries /*, other args */, Engine&);
+    // Neutron channels total kinetic energy
+    inline CELER_FUNCTION units::MevEnergy total_energy_neutron_channels() const
+    {
+        return units::MevEnergy{3.3};
+    }
+
+    // Proton channel total kinetic energy
+    inline CELER_FUNCTION units::MevEnergy total_energy_proton_channel() const
+    {
+        return units::MevEnergy{4.03};
+    }
+
+    // Outgoing neutron kinetic energy
+    inline CELER_FUNCTION units::MevEnergy neutron_kinetic_energy() const
+    {
+        return 0.75 * this->total_energy_neutron_channels();
+    }
+
+    // Outgoing proton kinetic energy
+    inline CELER_FUNCTION units::MevEnergy proton_kinetic_energy() const
+    {
+        return 0.75 * this->total_energy_proton_channel();
+    }
 };
 
 //---------------------------------------------------------------------------//
@@ -76,7 +102,7 @@ class DDMucfInteractor
  */
 CELER_FUNCTION
 DDMucfInteractor::DDMucfInteractor(NativeCRef<DTMixMucfData> const& data,
-                                   Channel const channel,
+                                   Channel channel,
                                    StackAllocator<Secondary>& allocate)
     : data_(data), channel_(channel), allocate_(allocate)
 {
@@ -92,36 +118,88 @@ template<class Engine>
 CELER_FUNCTION Interaction DDMucfInteractor::operator()(Engine& rng)
 {
     // Allocate space for the final fusion channel
-    Secondary* secondaries = allocate_(num_secondaries_[channel_]);
-    if (secondaries == nullptr)
+    Secondary* sec = allocate_(num_secondaries_[channel_]);
+    if (sec == nullptr)
     {
         // Failed to allocate space for secondaries
         return Interaction::from_failure();
     }
 
-    // Kill primary and generate secondaries
-    Interaction result = Interaction::from_absorption();
-    result.secondaries = this->sample_secondaries(secondaries, rng);
+    // Short aliases for secondary indices
+    size_type const neutron_idx{0}, proton_idx{0};
+    size_type const muon_idx{1}, muonic_he3{1};
+    size_type const he3_idx{2}, tritium_idx{2};
 
-    return result;
-}
+    IsotropicDistribution sample_isotropic;
 
-//---------------------------------------------------------------------------//
-/*!
- * Sample the secondaries of the selected channel.
- *
- * Since secondaries come from an at rest interaction, their final state is
- * a simple combination of random direction + momentum conservation
- */
-template<class Engine>
-CELER_FUNCTION Span<Secondary>
-DDMucfInteractor::sample_secondaries(Secondary* secondaries /*, other args */,
-                                     Engine&)
-{
-    // TODO: switch on channel_
-    CELER_ASSERT_UNREACHABLE();
+    switch (channel_)
+    {
+        case Channel::helium3_muon_neutron: {
+            // Neutron: random direction with known energy
+            sec[neutron_idx] = detail::sample_mucf_secondary(
+                data_.particle_ids.neutron, this->neutron_kinetic_energy(), rng);
+
+            // Muon: random direction with energy sampled from its CDF
+            sec[muon_idx] = detail::sample_mucf_muon(
+                data_.particle_ids.mu_minus,
+                NonuniformGridCalculator{data_.muon_energy_cdf, data_.reals},
+                rng);
+
+            // Helium-3: momentum conservation
+            sec[he3_idx]
+                = detail::calc_third_secondary(sec[neutron_idx],
+                                               data_.particle_masses.neutron,
+                                               sec[muon_idx],
+                                               data_.particle_masses.mu_minus,
+                                               data_.particle_ids.he3,
+                                               data_.particle_masses.he3);
+            break;
+        }
+
+        case Channel::muonichelium3_neutron: {
+            // Neutron: random direction with known energy
+            sec[neutron_idx] = detail::sample_mucf_secondary(
+                data_.particle_ids.neutron, this->neutron_kinetic_energy(), rng);
 
-    return Span<Secondary>{secondaries, num_secondaries_[channel_]};
+            // Muonic helium-3: momentum conservation
+            sec[muonic_he3].particle_id = data_.particle_ids.muonic_he3;
+            sec[muonic_he3].energy = this->total_energy_neutron_channels()
+                                     - this->neutron_kinetic_energy();
+            sec[muonic_he3].direction
+                = detail::opposite(sec[neutron_idx].direction);
+            break;
+        }
+
+        case Channel::tritium_muon_proton: {
+            // Proton: random direction with known energy
+            sec[proton_idx] = detail::sample_mucf_secondary(
+                data_.particle_ids.proton, this->proton_kinetic_energy(), rng);
+
+            // Muon: random direction with energy sampled from its CDF
+            sec[muon_idx] = detail::sample_mucf_muon(
+                data_.particle_ids.mu_minus,
+                NonuniformGridCalculator{data_.muon_energy_cdf, data_.reals},
+                rng);
+
+            // tritium: momentum conservation
+            sec[tritium_idx]
+                = detail::calc_third_secondary(sec[proton_idx],
+                                               data_.particle_masses.proton,
+                                               sec[muon_idx],
+                                               data_.particle_masses.mu_minus,
+                                               data_.particle_ids.triton,
+                                               data_.particle_masses.triton);
+            break;
+        }
+
+        default:
+            CELER_ASSERT_UNREACHABLE();
+    }
+
+    // Kill primary and generate secondaries
+    Interaction result = Interaction::from_absorption();
+    result.secondaries = {sec, num_secondaries_[channel_]};
+    return result;
 }
 
 //---------------------------------------------------------------------------//