Support decay in stepping loop

src/accel/SharedParams.cc

@@ -68,7 +68,6 @@
 #include "celeritas/phys/ParticleParams.hh"
 #include "celeritas/phys/PhysicsParams.hh"
 #include "celeritas/phys/Process.hh"
-#include "celeritas/phys/ProcessBuilder.hh"
 #include "celeritas/setup/FrameworkInput.hh"
 #include "celeritas/track/SimParams.hh"
 #include "celeritas/track/TrackInitParams.hh"

src/celeritas/CMakeLists.txt

@@ -36,6 +36,7 @@ list(APPEND PUBLIC_DEPS ${_core_geo_deps})
 list(APPEND SOURCES
   Types.cc
   TypesIO.json.cc
+  decay/DecayProcess.cc
   em/detail/Utils.cc
   em/params/AtomicRelaxationParams.cc
   em/params/FluctuationParams.cc
@@ -118,6 +119,7 @@ list(APPEND SOURCES
   optical/surface/SurfacePhysicsParams.cc
   optical/surface/SurfaceSteppingAction.cc
   phys/CutoffParams.cc
+  phys/detail/DecayTableInserter.cc
   phys/detail/EnergyMaxXsCalculator.cc
   phys/GeneratorInterface.cc
   phys/GeneratorRegistry.cc
@@ -328,6 +330,7 @@ celeritas_polysource(alongstep/AlongStepNeutralAction)
 celeritas_polysource(alongstep/AlongStepUniformMscAction)
 celeritas_polysource(alongstep/AlongStepRZMapFieldMscAction)
 celeritas_polysource(alongstep/AlongStepCylMapFieldMscAction)
+celeritas_polysource(decay/channel/MuDecayChannel)
 celeritas_polysource(em/model/BetheHeitlerModel)
 celeritas_polysource(em/model/BetheBlochModel)
 celeritas_polysource(em/model/BraggModel)

src/celeritas/Types.cc

@@ -90,6 +90,18 @@ char const* to_cstring(NuclearFormFactorType value)
     return to_cstring_impl(value);
 }
 
+//---------------------------------------------------------------------------//
+/*!
+ * Get a string corresponding to the decay channel.
+ */
+char const* to_cstring(DecayChannelType value)
+{
+    static EnumStringMapper<DecayChannelType> const to_cstring_impl{
+        "muon",
+    };
+    return to_cstring_impl(value);
+}
+
 //---------------------------------------------------------------------------//
 /*!
  * Get a string corresponding to the wavelength shifting time model selection.

src/celeritas/Types.hh

@@ -68,6 +68,9 @@ using TrackId = OpaqueId<struct Track_>;
 //! Opaque index of particle-nucleon cascade channel
 using ChannelId = OpaqueId<struct Channel_>;
 
+//! Opaque index of decay channel
+using DecayChannelId = OpaqueId<struct DecayChannelData>;
+
 //! Opaque index to one elemental component datum in a particular material
 using ElementComponentId = OpaqueId<struct MatElementComponent>;
 
@@ -202,6 +205,14 @@ enum class NuclearFormFactorType
     size_
 };
 
+//---------------------------------------------------------------------------//
+//! Decay channel type
+enum class DecayChannelType
+{
+    muon,
+    size_
+};
+
 //---------------------------------------------------------------------------//
 //! Optical photon wavelength shifting time model
 enum class WlsTimeProfile
@@ -277,6 +288,9 @@ char const* to_cstring(MscStepLimitAlgorithm value);
 // Get a string corresponding to the nuclear form factor model
 char const* to_cstring(NuclearFormFactorType value);
 
+// Get a string corresponding to the decay channel
+char const* to_cstring(DecayChannelType value);
+
 // Get a string corresponding to the interpolation method
 char const* to_cstring(InterpolationType value);
 

src/celeritas/decay/DecayData.hh

@@ -0,0 +1,54 @@
+//------------------------------- -*- C++ -*- -------------------------------//
+// Copyright Celeritas contributors: see top-level COPYRIGHT file for details
+// SPDX-License-Identifier: (Apache-2.0 OR MIT)
+//---------------------------------------------------------------------------//
+//! \file celeritas/decay/DecayData.hh
+//---------------------------------------------------------------------------//
+#pragma once
+
+#include "corecel/Macros.hh"
+#include "corecel/Types.hh"
+#include "corecel/data/Collection.hh"
+#include "celeritas/Types.hh"
+
+namespace celeritas
+{
+//---------------------------------------------------------------------------//
+/*!
+ * Data for a decay interactor.
+ *
+ * This stores the particle IDs of the daughters for a specific decay channel.
+ */
+struct DecayChannelData
+{
+    //! Daughter particle IDs
+    ItemRange<ParticleId> daughters;
+
+    //! Whether the data is assigned
+    explicit CELER_FUNCTION operator bool() const
+    {
+        return !daughters.empty();
+    }
+};
+
+//---------------------------------------------------------------------------//
+/*!
+ * Decay channels for a particle type.
+ */
+struct DecayTableData
+{
+    //! Decay channels
+    ItemRange<DecayChannelId> channel_ids;
+    //! Branching ratio of each decay channel
+    ItemRange<real_type> branching_ratios;
+
+    //! Whether the data is assigned
+    explicit CELER_FUNCTION operator bool() const
+    {
+        return !channel_ids.empty()
+               && branching_ratios.size() == channel_ids.size();
+    }
+};
+
+//---------------------------------------------------------------------------//
+}  // namespace celeritas

src/celeritas/decay/DecayProcess.cc

@@ -0,0 +1,90 @@
+//------------------------------- -*- C++ -*- -------------------------------//
+// Copyright Celeritas contributors: see top-level COPYRIGHT file for details
+// SPDX-License-Identifier: (Apache-2.0 OR MIT)
+//---------------------------------------------------------------------------//
+//! \file celeritas/decay/DecayProcess.cc
+//---------------------------------------------------------------------------//
+#include "DecayProcess.hh"
+
+#include <set>
+
+#include "corecel/Assert.hh"
+#include "celeritas/decay/DecayData.hh"
+#include "celeritas/decay/channel/MuDecayChannel.hh"
+#include "celeritas/phys/PDGNumber.hh"
+#include "celeritas/phys/ParticleParams.hh"
+
+namespace celeritas
+{
+//---------------------------------------------------------------------------//
+/*!
+ * Construct from particles and input data.
+ */
+DecayProcess::DecayProcess(SPConstParticles particles,
+                           inp::DecayPhysics const& input)
+    : particles_(std::move(particles)), input_(input)
+{
+    CELER_EXPECT(particles_);
+    CELER_VALIDATE(input_, << "no decay tables are present");
+}
+
+//---------------------------------------------------------------------------//
+/*!
+ * Construct the decay channels.
+ */
+auto DecayProcess::build_channels(ActionIdIter start_id) const -> VecChannel
+{
+    VecChannel result;
+    std::set<DecayChannelType> types;
+    for (auto const& [pdg, table] : input_.tables)
+    {
+        for (auto const& channel : table)
+        {
+            // Identify the unique channel types
+            auto [iter, inserted] = types.insert(channel.type);
+            if (inserted)
+            {
+                // Construct an action for each channel
+                switch (*iter)
+                {
+                    case DecayChannelType::muon:
+                        result.push_back(std::make_shared<MuDecayChannel>(
+                            *start_id++, input_));
+                        break;
+                    default:
+                        CELER_NOT_IMPLEMENTED("Decay channel type");
+                }
+            }
+        }
+    }
+    CELER_VALIDATE(!result.empty(),
+                   << "no supported channels for decay process");
+    return result;
+}
+
+//---------------------------------------------------------------------------//
+/*!
+ * Get the decay table for a particle.
+ */
+auto DecayProcess::decay_table(ParticleId pid) const -> DecayTable
+{
+    auto iter = input_.tables.find(particles_->id_to_pdg(pid));
+    if (iter == input_.tables.end())
+    {
+        return {};
+    }
+    return iter->second;
+}
+
+//---------------------------------------------------------------------------//
+/*!
+ * Whether the decay process applies to the particle.
+ */
+bool DecayProcess::is_applicable(ParticleId pid) const
+{
+    return input_.tables.find(particles_->id_to_pdg(pid))
+           != input_.tables.end();
+}
+
+//---------------------------------------------------------------------------//
+}  // namespace celeritas

src/celeritas/decay/DecayProcess.hh

@@ -0,0 +1,64 @@
+//------------------------------- -*- C++ -*- -------------------------------//
+// Copyright Celeritas contributors: see top-level COPYRIGHT file for details
+// SPDX-License-Identifier: (Apache-2.0 OR MIT)
+//---------------------------------------------------------------------------//
+//! \file celeritas/decay/DecayProcess.hh
+//---------------------------------------------------------------------------//
+#pragma once
+
+#include <memory>
+
+#include "celeritas/Types.hh"
+#include "celeritas/decay/channel/DecayChannel.hh"
+#include "celeritas/inp/Physics.hh"
+#include "celeritas/phys/ImportedProcessAdapter.hh"
+#include "celeritas/phys/Process.hh"
+
+namespace celeritas
+{
+class ParticleParams;
+
+//---------------------------------------------------------------------------//
+/*!
+ * Process for decay.
+ */
+class DecayProcess final : public Process
+{
+  public:
+    //!@{
+    //! \name Type aliases
+    using SPConstChannel = std::shared_ptr<DecayChannel const>;
+    using SPConstParticles = std::shared_ptr<ParticleParams const>;
+    using VecChannel = std::vector<SPConstChannel>;
+    using DecayTable = inp::DecayPhysics::DecayTable;
+    //!@}
+
+  public:
+    // Construct from particles and input data
+    DecayProcess(SPConstParticles, inp::DecayPhysics const&);
+
+    // Construct the decay channels
+    VecChannel build_channels(ActionIdIter start_id) const;
+
+    // Get the decay table for a particle
+    DecayTable decay_table(ParticleId) const;
+
+    // Whether the decay process applies to the particle
+    bool is_applicable(ParticleId) const;
+
+    //!@{
+    //! \name Process interface
+
+    //! Whether the process applies when the particle is stopped
+    bool applies_at_rest() const final { return true; }
+    //! Name of the process
+    std::string_view label() const final { return "Decay"; }
+    //!@}
+
+  private:
+    SPConstParticles particles_;
+    inp::DecayPhysics const& input_;
+};
+
+//---------------------------------------------------------------------------//
+}  // namespace celeritas

src/celeritas/decay/channel/DecayChannel.hh

@@ -0,0 +1,30 @@
+//------------------------------- -*- C++ -*- -------------------------------//
+// Copyright Celeritas contributors: see top-level COPYRIGHT file for details
+// SPDX-License-Identifier: (Apache-2.0 OR MIT)
+//---------------------------------------------------------------------------//
+//! \file celeritas/decay/channel/DecayChannel.hh
+//---------------------------------------------------------------------------//
+#pragma once
+
+#include "corecel/Types.hh"
+#include "celeritas/Types.hh"
+#include "celeritas/global/ActionInterface.hh"
+
+namespace celeritas
+{
+//---------------------------------------------------------------------------//
+/*!
+ * Abstract base class representing a decay channel.
+ *
+ * This class is responsible for validating the input decay table data and
+ * launching a kernel for the decay interactor.
+ */
+class DecayChannel : public CoreStepActionInterface
+{
+  public:
+    //! Dependency ordering of the action
+    StepActionOrder order() const final { return StepActionOrder::post; }
+};
+
+//---------------------------------------------------------------------------//
+}  // namespace celeritas