Add CHIPS gamma-nuclear cross sections

doc/_static/zotero.bib

@@ -3718,3 +3718,15 @@ @article{wright-geant4bertini-2015
   issn = {01689002},
   doi = {10.1016/j.nima.2015.09.058}
 }
+
+@article{chips-gamma-nuclear-xs-2000,
+  title = {The Chiral invariant phase space event generator, iii: modeling of
+real and virtual photon interactionswith nuclei below pion production threshold.},
+  author = {Degtyarenko P.V., Kossov M.V. and Wellisch, H.-P.},
+  year = 2000,
+  month = dec,
+  journal = {European Physical Journal A},
+  volume = {9},
+  pages = {421--424},
+  doi = {10.1007/s100500070026}
+}

src/celeritas/CMakeLists.txt

@@ -201,6 +201,7 @@ endif()
 if(CELERITAS_USE_Geant4)
   set(_cg4_sources
     g4/DetectorConstruction.cc
+    g4/EmExtraPhysicsHelper.cc
     g4/SupportedEmStandardPhysics.cc
     g4/SupportedOpticalPhysics.cc
     g4/detail/GeantBremsstrahlungProcess.cc

src/celeritas/em/data/GammaNuclearData.hh

@@ -28,13 +28,7 @@ struct GammaNuclearScalars
     //! Model's maximum energy limit [MeV]
     static CELER_CONSTEXPR_FUNCTION units::MevEnergy max_valid_energy()
     {
-        return units::MevEnergy{5e+4};
-    }
-
-    //! Maximum energy limit for PARTICLEXS/gamma cross sections [MeV]
-    static CELER_CONSTEXPR_FUNCTION units::MevEnergy max_low_energy()
-    {
-        return units::MevEnergy{130};
+        return units::MevEnergy{1e+8};
     }
 
     //! Whether data are assigned
@@ -61,16 +55,18 @@ struct GammaNuclearData
     // Scalar data
     GammaNuclearScalars scalars;
 
-    // Microscopic (element) cross section data (G4PARTICLEXS/gamma/inelZ)
-    ElementItems<NonuniformGridRecord> micro_xs;
-
-    // Backend data
+    // Microscopic cross sections using G4PARTICLEXS/gamma nuclear (IAEA) data
+    ElementItems<NonuniformGridRecord> xs_iaea;
     Items<real_type> reals;
 
+    // Microscopic cross sections using parameterized CHIPS data at high energy
+    ElementItems<NonuniformGridRecord> xs_chips;
+
     //! Whether the data are assigned
     explicit CELER_FUNCTION operator bool() const
     {
-        return scalars && !micro_xs.empty() && !reals.empty();
+        return scalars && !xs_iaea.empty() && !reals.empty()
+               && !xs_chips.empty();
     }
 
     //! Assign from another set of data
@@ -79,7 +75,8 @@ struct GammaNuclearData
     {
         CELER_EXPECT(other);
         scalars = other.scalars;
-        micro_xs = other.micro_xs;
+        xs_iaea = other.xs_iaea;
+        xs_chips = other.xs_chips;
         reals = other.reals;
 
         return *this;

src/celeritas/em/model/GammaNuclearModel.cc

@@ -6,7 +6,9 @@
 //---------------------------------------------------------------------------//
 #include "GammaNuclearModel.hh"
 
+#include "corecel/grid/VectorUtils.hh"
 #include "corecel/math/Quantity.hh"
+#include "celeritas/g4/EmExtraPhysicsHelper.hh"
 #include "celeritas/global/CoreParams.hh"
 #include "celeritas/global/CoreState.hh"
 #include "celeritas/grid/NonuniformGridInserter.hh"
@@ -32,6 +34,8 @@ GammaNuclearModel::GammaNuclearModel(ActionId id,
 
     HostVal<GammaNuclearData> data;
 
+    helper_ = std::make_shared<EmExtraPhysicsHelper>();
+
     // Save IDs
     data.scalars.gamma_id = particles.find(pdg::gamma());
 
@@ -43,13 +47,22 @@ GammaNuclearModel::GammaNuclearModel(ActionId id,
     CELER_EXPECT(data.scalars);
 
     // Load gamma-nuclear element cross section data
-    NonuniformGridInserter insert_micro_xs{&data.reals, &data.micro_xs};
+    NonuniformGridInserter insert_xs_iaea{&data.reals, &data.xs_iaea};
+    NonuniformGridInserter insert_xs_chips{&data.reals, &data.xs_chips};
+
+    double const emax = data.scalars.max_valid_energy().value();
     for (auto el_id : range(ElementId{materials.num_elements()}))
     {
         AtomicNumber z = materials.get(el_id).atomic_number();
-        insert_micro_xs(load_data(z));
+        // Build element cross sections from G4PARTICLEXS
+        insert_xs_iaea(load_data(z));
+
+        // Build element cross sections above the upper bound of G4PARTICLEXS
+        double emin = data.reals[data.xs_iaea[el_id].grid.back()];
+        insert_xs_chips(this->calc_chips_xs(z, emin, emax));
     }
-    CELER_ASSERT(data.micro_xs.size() == materials.num_elements());
+    CELER_ASSERT(data.xs_iaea.size() == materials.num_elements());
+    CELER_ASSERT(data.xs_iaea.size() == data.xs_chips.size());
 
     // Move to mirrored data, copying to device
     data_ = CollectionMirror<GammaNuclearData>{std::move(data)};
@@ -98,5 +111,51 @@ void GammaNuclearModel::step(CoreParams const&, CoreStateDevice&) const
 }
 #endif
 
+//---------------------------------------------------------------------------//
+/*!
+ * Build CHIPS gamma-nuclear element cross sections using G4GammaNuclearXS.
+ */
+inp::Grid GammaNuclearModel::calc_chips_xs(AtomicNumber atomic_number,
+                                           double emin,
+                                           double emax) const
+{
+    CELER_EXPECT(atomic_number);
+
+    inp::Grid result;
+
+    // Tabulate cross sections using separate parameterizations for the high
+    // energy region (emin < E < 50 GeV) and the ultra high energy region up
+    // to the maximum valid energy (emax). The numbers of bins are chosen to
+    // adequately capture both the parameterized points (224 bins from 106 MeV
+    // to 50 GeV) and the calculations used in G4PhotoNuclearCrossSection,
+    // which can be made configurable if needed (TODO). Note that the linear
+    // interpolation between the upper limit of the IAEA cross-section data
+    // and 150 MeV, as used in G4GammaNuclearXS, is also included in this
+    // tabulation.
+    double const emid = helper_->max_high_energy();
+    size_type nbin_total = 300;
+    size_type nbin_ultra = 50;
+
+    result.x.resize(nbin_total);
+    result.y.resize(nbin_total);
+
+    result.x = geomspace(emin, emid, nbin_total - nbin_ultra);
+    result.x.pop_back();
+    auto ultra = geomspace(emid, emax, nbin_ultra + 1);
+    result.x.insert(result.x.end(), ultra.begin(), ultra.end());
+
+    // Tabulate the cross section from emin to emax
+    Quantity<UnitProduct<units::Millimeter, units::Millimeter>, double> xs;
+    for (size_type i = 0; i < nbin_total; ++i)
+    {
+        xs.value()
+            = helper_->GammaNuclearElementXS(result.x[i], atomic_number.get());
+        result.y[i]
+            = native_value_to<units::BarnXs>(native_value_from(xs)).value();
+    }
+
+    return result;
+}
+
 //---------------------------------------------------------------------------//
 }  // namespace celeritas

src/celeritas/em/model/GammaNuclearModel.hh

@@ -7,6 +7,7 @@
 #pragma once
 
 #include <functional>
+#include <memory>
 
 #include "corecel/data/CollectionMirror.hh"
 #include "corecel/inp/Grid.hh"
@@ -21,9 +22,23 @@ namespace celeritas
 class MaterialParams;
 class ParticleParams;
 
+class EmExtraPhysicsHelper;
+
 //---------------------------------------------------------------------------//
 /*!
  * Set up and launch the gamma-nuclear model interaction.
+ *
+ * The class also builds element cross-section tables using G4PARTICLEXS/gamma
+ * (IAEA) data at low energies and CHIPS gamma–nuclear cross sections using
+ * G4GammaNuclearXS above the IAEA upper energy limit (~130 MeV). The CHIPS
+ * cross sections are based on the parameterization developed by M. V. Kossov
+ * (CERN/ITEP Moscow) for the high energy region (106 MeV < E < 50 GeV) and on
+ * a Reggeon-based parameterization for the ultra high energy region
+ * (E > 50 GeV), as described in
+ * \citet{chips-gamma-nuclear-xs-2000, https://doi.org/10.1007/s100500070026}.
+ * G4GammaNuclearXS uses CHIPS (G4PhotoNuclearCrossSection) above 150 MeV and
+ * performs linear interpolation between the upper limit of the G4PARTICLEXS
+ * gamma-nuclear (IAEA) data and 150 MeV.
  */
 class GammaNuclearModel final : public Model, public StaticConcreteAction
 {
@@ -62,13 +77,19 @@ class GammaNuclearModel final : public Model, public StaticConcreteAction
 
   private:
     //// DATA ////
+    std::shared_ptr<EmExtraPhysicsHelper> helper_;
 
     // Host/device storage and reference
     CollectionMirror<GammaNuclearData> data_;
 
     //// TYPES ////
 
     using HostXsData = HostVal<GammaNuclearData>;
+
+    //// HELPER FUNCTIONS ////
+
+    inp::Grid
+    calc_chips_xs(AtomicNumber atomic_number, double emin, double emax) const;
 };
 
 //---------------------------------------------------------------------------//

src/celeritas/em/xs/GammaNuclearMicroXsCalculator.hh

@@ -38,8 +38,8 @@ class GammaNuclearMicroXsCalculator
     inline CELER_FUNCTION BarnXs operator()(ElementId el_id) const;
 
   private:
-    // Shared constant physics properties
-    ParamsRef const& shared_;
+    // Shared cross section data
+    ParamsRef const& data_;
     // Incident photon energy
     real_type const inc_energy_;
 };
@@ -52,25 +52,36 @@ class GammaNuclearMicroXsCalculator
  */
 CELER_FUNCTION
 GammaNuclearMicroXsCalculator::GammaNuclearMicroXsCalculator(
-    ParamsRef const& shared, Energy energy)
-    : shared_(shared), inc_energy_(energy.value())
+    ParamsRef const& data, Energy energy)
+    : data_(data), inc_energy_(energy.value())
 {
 }
 
 //---------------------------------------------------------------------------//
 /*!
- * Compute microscopic (element) cross section
+ * Compute microscopic gamma-nuclear cross section at the given gamma energy.
  */
 CELER_FUNCTION
 auto GammaNuclearMicroXsCalculator::operator()(ElementId el_id) const -> BarnXs
 {
-    CELER_EXPECT(el_id < shared_.micro_xs.size());
+    NonuniformGridRecord grid;
 
-    // Get element cross section data
-    NonuniformGridRecord grid = shared_.micro_xs[el_id];
+    // Use G4PARTICLEXS gamma-nuclear cross sections at low energy and CHIPS
+    // parameterized cross sections at high energy.
+
+    if (inc_energy_ <= data_.reals[data_.xs_iaea[el_id].grid.back()])
+    {
+        CELER_EXPECT(el_id < data_.xs_iaea.size());
+        grid = data_.xs_iaea[el_id];
+    }
+    else
+    {
+        CELER_EXPECT(el_id < data_.xs_chips.size());
+        grid = data_.xs_chips[el_id];
+    }
 
     // Calculate micro cross section at the given energy
-    NonuniformGridCalculator calc_xs(grid, shared_.reals);
+    NonuniformGridCalculator calc_xs(grid, data_.reals);
     real_type result = calc_xs(inc_energy_);
 
     return BarnXs{result};

src/celeritas/g4/EmExtraPhysicsHelper.cc

@@ -0,0 +1,33 @@
+//------------------------------- -*- C++ -*- -------------------------------//
+// Copyright Celeritas contributors: see top-level COPYRIGHT file for details
+// SPDX-License-Identifier: (Apache-2.0 OR MIT)
+//---------------------------------------------------------------------------//
+//! \file celeritas/g4/EmExtraPhysicsHelper.cc
+//---------------------------------------------------------------------------//
+#include "EmExtraPhysicsHelper.hh"
+
+#include <G4GammaNuclearXS.hh>
+
+namespace celeritas
+{
+//---------------------------------------------------------------------------//
+/*!
+ * Construct with Geant4 provided cross section classes
+ */
+EmExtraPhysicsHelper::EmExtraPhysicsHelper()
+{
+    gn_xs_ = std::make_shared<G4GammaNuclearXS>();
+}
+
+//---------------------------------------------------------------------------//
+/*!
+ * Calculate the gamma-nuclear element cross section using G4GammaNuclearXS
+ * in the native Geant4 unit [mb]
+ */
+double EmExtraPhysicsHelper::GammaNuclearElementXS(double energy, int z)
+{
+    return gn_xs_->ElementCrossSection(energy, z);
+}
+
+//---------------------------------------------------------------------------//
+}  // namespace celeritas

src/celeritas/g4/EmExtraPhysicsHelper.hh

@@ -0,0 +1,41 @@
+//------------------------------- -*- C++ -*- -------------------------------//
+// Copyright Celeritas contributors: see top-level COPYRIGHT file for details
+// SPDX-License-Identifier: (Apache-2.0 OR MIT)
+//---------------------------------------------------------------------------//
+//! \file celeritas/g4/EmExtraPhysicsHelper.hh
+//---------------------------------------------------------------------------//
+#pragma once
+
+#include <memory>
+
+class G4GammaNuclearXS;
+
+namespace celeritas
+{
+
+//---------------------------------------------------------------------------//
+/*!
+ * A helper class to interface with Geant4 cross sections.
+ */
+class EmExtraPhysicsHelper
+{
+  public:
+    // Construct EM extra physics helper
+    EmExtraPhysicsHelper();
+
+    // Calculate gamma-nuclear element cross section
+    double GammaNuclearElementXS(double energy, int z);
+
+    // The maximum high energy of G4PhotoNuclearCrossSection
+    static constexpr double max_high_energy()
+    {
+        return 5e+4;  // clhep::MeV
+    }
+
+  private:
+    //// DATA ////
+    std::shared_ptr<G4GammaNuclearXS> gn_xs_;
+};
+
+//---------------------------------------------------------------------------//
+}  // namespace celeritas