Add neutron capture and inelastic processes

src/celeritas/CMakeLists.txt

@@ -62,7 +62,9 @@ list(APPEND SOURCES
   mat/MaterialParams.cc
   mat/MaterialParamsOutput.cc
   mat/detail/Utils.cc
+  neutron/process/NeutronCaptureProcess.cc
   neutron/process/NeutronElasticProcess.cc
+  neutron/process/NeutronInelasticProcess.cc
   optical/CerenkovParams.cc
   optical/ScintillationParams.cc
   optical/OpticalPropertyParams.cc
@@ -254,7 +256,9 @@ celeritas_polysource(global/alongstep/AlongStepGeneralLinearAction)
 celeritas_polysource(global/alongstep/AlongStepNeutralAction)
 celeritas_polysource(global/alongstep/AlongStepUniformMscAction)
 celeritas_polysource(global/alongstep/AlongStepRZMapFieldMscAction)
-celeritas_polysource(neutron/model/ChipsNeutronElasticModel)
+celeritas_polysource(neutron/model/NeutronCaptureModel)
+celeritas_polysource(neutron/model/NeutronElasticModel)
+celeritas_polysource(neutron/model/NeutronInelasticModel)
 celeritas_polysource(phys/detail/DiscreteSelectAction)
 celeritas_polysource(phys/detail/PreStepAction)
 celeritas_polysource(random/RngReseed)

src/celeritas/io/ImportData.hh

@@ -58,7 +58,9 @@ struct ImportData
     using ImportSBMap = std::map<ZInt, ImportSBTable>;
     using ImportLivermorePEMap = std::map<ZInt, ImportLivermorePE>;
     using ImportAtomicRelaxationMap = std::map<ZInt, ImportAtomicRelaxation>;
+    using ImportNeutronCaptureMap = std::map<ZInt, ImportPhysicsVector>;
     using ImportNeutronElasticMap = std::map<ZInt, ImportPhysicsVector>;
+    using ImportNeutronInelasticMap = std::map<ZInt, ImportPhysicsVector>;
     using ImportOpticalMap = std::map<MatIdx, ImportOpticalMaterial>;
     //!@}
 
@@ -75,7 +77,9 @@ struct ImportData
     ImportOpticalParameters optical_params;
     ImportSBMap sb_data;
     ImportLivermorePEMap livermore_pe_data;
+    ImportNeutronCaptureMap neutron_capture_data;
     ImportNeutronElasticMap neutron_elastic_data;
+    ImportNeutronInelasticMap neutron_inelastic_data;
     ImportAtomicRelaxationMap atomic_relaxation_data;
 
     std::string units;  //!< "cgs", "clhep", or "si"

src/celeritas/io/ImportProcess.cc

@@ -63,7 +63,9 @@ char const* to_cstring(ImportProcessClass value)
         "mu_brems",
         "mu_pair_prod",
         "gamma_general",
+        "neutron_capture",
         "neutron_elastic",
+        "neutron_inelastic",
     };
     return to_cstring_impl(value);
 }
@@ -93,7 +95,9 @@ char const* to_geant_name(ImportProcessClass value)
         "muBrems",  // mu_brems,
         "muPairProd",  // mu_pair_prod,
         "GammaGeneralProc",  // gamma_general,
+        "neutronCaptureProc",  // neutron_capture,
         "neutronElasticProc",  // neutron_elastic,
+        "neutronInelasticProc",  // neutron_inelastic,
     };
     return to_name_impl(value);
 }

src/celeritas/io/ImportProcess.hh

@@ -72,7 +72,9 @@ enum class ImportProcessClass
     mu_pair_prod,
     gamma_general,  // Will be decomposed into other processes
     // Neutron
+    neutron_capture,
     neutron_elastic,
+    neutron_inelastic,
     size_
 };
 

src/celeritas/io/NeutronXsReader.cc

@@ -12,8 +12,8 @@
 
 #include "corecel/Assert.hh"
 #include "corecel/Types.hh"
+#include "corecel/io/EnumStringMapper.hh"
 #include "corecel/io/Logger.hh"
-#include "corecel/math/Algorithms.hh"
 #include "corecel/sys/Environment.hh"
 #include "celeritas/Quantities.hh"
 #include "celeritas/Units.hh"
@@ -26,20 +26,21 @@ namespace celeritas
  * Construct the reader using the G4PARTICLEXSDATA environment variable to get
  * the path to the data.
  */
-NeutronXsReader::NeutronXsReader()
+NeutronXsReader::NeutronXsReader(NeutronXsType type) : type_(type)
 {
     std::string const& dir = celeritas::getenv("G4PARTICLEXSDATA");
     CELER_VALIDATE(!dir.empty(),
                    << "environment variable G4PARTICLEXSDATA is not defined "
-                      "(needed to locate neutron elastic cross section data)");
+                      "(needed to locate neutron cross section data)");
     path_ = dir + "/neutron";
 }
 
 //---------------------------------------------------------------------------//
 /*!
  * Construct the reader with the path to the directory containing the data.
  */
-NeutronXsReader::NeutronXsReader(char const* path) : path_(path)
+NeutronXsReader::NeutronXsReader(char const* path, NeutronXsType type)
+    : path_(path), type_(type)
 {
     CELER_EXPECT(!path_.empty());
     if (path_.back() == '/')
@@ -58,13 +59,14 @@ NeutronXsReader::operator()(AtomicNumber atomic_number) const
     CELER_EXPECT(atomic_number);
 
     std::string z_str = std::to_string(atomic_number.unchecked_get());
-    CELER_LOG(debug) << "Reading neutron elastic xs data for Z=" << z_str;
+    CELER_LOG(debug) << "Reading neutron xs data for " << to_cstring(type_)
+                     << " Z=" << z_str;
 
     result_type result;
 
     // Read neutron elastic cross section data for the given atomic_number
     {
-        std::string filename = path_ + "/el" + z_str;
+        std::string filename = path_ + "/" + to_cstring(type_) + z_str;
         std::ifstream infile(filename);
         CELER_VALIDATE(infile,
                        << "failed to open '" << filename
@@ -99,5 +101,20 @@ NeutronXsReader::operator()(AtomicNumber atomic_number) const
     return result;
 }
 
+//---------------------------------------------------------------------------//
+// FREE FUNCTIONS
+//---------------------------------------------------------------------------//
+
+//---------------------------------------------------------------------------//
+/*!
+ * Get the string value for a neutron cross section data type.
+ */
+char const* to_cstring(NeutronXsType value)
+{
+    static EnumStringMapper<NeutronXsType> const to_cstring_impl{
+        "cap", "el", "inel"};
+    return to_cstring_impl(value);
+}
+
 //---------------------------------------------------------------------------//
 }  // namespace celeritas

src/celeritas/io/NeutronXsReader.hh

@@ -21,7 +21,20 @@ namespace celeritas
 {
 //---------------------------------------------------------------------------//
 /*!
- * Load the neutron elastic cross section (G4PARTICLEXSDATA/neutron/elZ) data.
+ * Types of microscopic cross sections in G4PARTICLEXSDATA/neutron data.
+ */
+enum class NeutronXsType
+{
+    cap,  //!< Capture cross section
+    el,  //!< Elastic cross section
+    inel,  //!< Inelastic cross section
+    size_
+};
+
+//---------------------------------------------------------------------------//
+/*!
+ * Load the neutron cross section (G4PARTICLEXSDATA/neutron) data by the
+ * interaction type (capture, elastic, and inelastic).
  */
 class NeutronXsReader
 {
@@ -35,18 +48,29 @@ class NeutronXsReader
 
   public:
     // Construct the reader and locate the data using the environment variable
-    NeutronXsReader();
+    explicit NeutronXsReader(NeutronXsType type);
 
-    // Construct the reader from the path to the data directory
-    explicit NeutronXsReader(char const* path);
+    // Construct the reader from the path to the data directory and the type
+    explicit NeutronXsReader(char const* path, NeutronXsType type);
 
     // Read the data for the given element
     result_type operator()(AtomicNumber atomic_number) const;
 
+  private:
+    // Get the string value for the cross section data type
+
   private:
     // Directory containing the neutron elastic cross section data
     std::string path_;
+    NeutronXsType type_;
 };
 
+//---------------------------------------------------------------------------//
+// FREE FUNCTIONS
+//---------------------------------------------------------------------------//
+
+// Get the string value for a neutron cross section type
+char const* to_cstring(NeutronXsType value);
+
 //---------------------------------------------------------------------------//
 }  // namespace celeritas

src/celeritas/phys/PhysicsData.hh

@@ -18,7 +18,9 @@
 #include "celeritas/em/data/LivermorePEData.hh"
 #include "celeritas/grid/ValueGridType.hh"
 #include "celeritas/grid/XsGridData.hh"
+#include "celeritas/neutron/data/NeutronCaptureData.hh"
 #include "celeritas/neutron/data/NeutronElasticData.hh"
+#include "celeritas/neutron/data/NeutronInelasticData.hh"
 
 #include "Interaction.hh"
 #include "Secondary.hh"
@@ -174,11 +176,21 @@ struct HardwiredModels
     ModelId eplusgg;
     EPlusGGData eplusgg_data;
 
+    // Neutron capture
+    ProcessId neutron_capture;
+    ModelId rad_capture;
+    NeutronCaptureData<W, M> capture_data;
+
     // Neutron elastic
     ProcessId neutron_elastic;
     ModelId chips;
     NeutronElasticData<W, M> chips_data;
 
+    // Neutron inelastic
+    ProcessId neutron_inelastic;
+    ModelId bertini;
+    NeutronInelasticData<W, M> bertini_data;
+
     //// MEMBER FUNCTIONS ////
 
     //! Assign from another set of hardwired models
@@ -199,13 +211,29 @@ struct HardwiredModels
         eplusgg = other.eplusgg;
         eplusgg_data = other.eplusgg_data;
 
+        neutron_capture = other.neutron_capture;
+        if (neutron_capture)
+        {
+            // Only assign neutron_cpature data if that process is present
+            rad_capture = other.rad_capture;
+            capture_data = other.capture_data;
+        }
         neutron_elastic = other.neutron_elastic;
+
         if (neutron_elastic)
         {
             // Only assign neutron_elastic data if that process is present
             chips = other.chips;
             chips_data = other.chips_data;
         }
+        neutron_inelastic = other.neutron_inelastic;
+
+        if (neutron_inelastic)
+        {
+            // Only assign neutron_inelastic data if that process is present
+            bertini = other.bertini;
+            bertini_data = other.bertini_data;
+        }
 
         return *this;
     }

src/celeritas/phys/PhysicsParams.cc

@@ -43,8 +43,12 @@
 #include "celeritas/mat/MaterialData.hh"
 #include "celeritas/mat/MaterialParams.hh"
 #include "celeritas/mat/MaterialView.hh"
+#include "celeritas/neutron/data/NeutronCaptureData.hh"
 #include "celeritas/neutron/data/NeutronElasticData.hh"
-#include "celeritas/neutron/model/ChipsNeutronElasticModel.hh"
+#include "celeritas/neutron/data/NeutronInelasticData.hh"
+#include "celeritas/neutron/model/NeutronCaptureModel.hh"
+#include "celeritas/neutron/model/NeutronElasticModel.hh"
+#include "celeritas/neutron/model/NeutronInelasticModel.hh"
 
 #include "Model.hh"
 #include "ParticleParams.hh"
@@ -419,12 +423,26 @@ void PhysicsParams::build_ids(ParticleParams const& particles,
             data->hardwired.eplusgg = ModelId{model_idx};
             data->hardwired.eplusgg_data = epgg_model->device_ref();
         }
-        else if (auto* ne_model
-                 = dynamic_cast<ChipsNeutronElasticModel const*>(&model))
+        else if (auto* ncap_model
+                 = dynamic_cast<NeutronCaptureModel const*>(&model))
+        {
+            data->hardwired.neutron_capture = process_id;
+            data->hardwired.rad_capture = ModelId{model_idx};
+            data->hardwired.capture_data = ncap_model->device_ref();
+        }
+        else if (auto* nel_model
+                 = dynamic_cast<NeutronElasticModel const*>(&model))
         {
             data->hardwired.neutron_elastic = process_id;
             data->hardwired.chips = ModelId{model_idx};
-            data->hardwired.chips_data = ne_model->device_ref();
+            data->hardwired.chips_data = nel_model->device_ref();
+        }
+        else if (auto* ninel_model
+                 = dynamic_cast<NeutronInelasticModel const*>(&model))
+        {
+            data->hardwired.neutron_inelastic = process_id;
+            data->hardwired.bertini = ModelId{model_idx};
+            data->hardwired.bertini_data = ninel_model->device_ref();
         }
     }
 

src/celeritas/phys/PhysicsTrackView.hh

@@ -19,7 +19,9 @@
 #include "celeritas/grid/XsCalculator.hh"
 #include "celeritas/mat/MaterialView.hh"
 #include "celeritas/mat/TabulatedElementSelector.hh"
+#include "celeritas/neutron/xs/NeutronCaptureMicroXsCalculator.hh"
 #include "celeritas/neutron/xs/NeutronElasticMicroXsCalculator.hh"
+#include "celeritas/neutron/xs/NeutronInelasticMicroXsCalculator.hh"
 #include "celeritas/phys/MacroXsCalculator.hh"
 
 #include "PhysicsData.hh"
@@ -418,12 +420,24 @@ CELER_FUNCTION real_type PhysicsTrackView::calc_xs(ParticleProcessId ppid,
                 params_.hardwired.eplusgg_data, material);
             result = calc_xs(energy);
         }
+        else if (model_id == params_.hardwired.rad_capture)
+        {
+            auto calc_xs = MacroXsCalculator<NeutronCaptureMicroXsCalculator>(
+                params_.hardwired.capture_data, material);
+            result = calc_xs(energy);
+        }
         else if (model_id == params_.hardwired.chips)
         {
             auto calc_xs = MacroXsCalculator<NeutronElasticMicroXsCalculator>(
                 params_.hardwired.chips_data, material);
             result = calc_xs(energy);
         }
+        else if (model_id == params_.hardwired.bertini)
+        {
+            auto calc_xs = MacroXsCalculator<NeutronInelasticMicroXsCalculator>(
+                params_.hardwired.bertini_data, material);
+            result = calc_xs(energy);
+        }
     }
     else if (auto grid_id = this->value_grid(ValueGridType::macro_xs, ppid))
     {