Add CHIPS gamma-nuclear cross sections

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,19 @@ struct GammaNuclearData
     // Scalar data
     GammaNuclearScalars scalars;
 
-    // Microscopic (element) cross section data (G4PARTICLEXS/gamma/inelZ)
-    ElementItems<NonuniformGridRecord> micro_xs;
+    // Microscopic cross sections using G4PARTICLEXS/gamma nuclear (IAEA) data
+    ElementItems<NonuniformGridRecord> micro_xs_iaea;
+    Items<real_type> reals_iaea;
 
-    // Backend data
-    Items<real_type> reals;
+    // Microscopic cross sections using parameterized CHIPS data at high energy
+    ElementItems<NonuniformGridRecord> micro_xs_chips;
+    Items<real_type> reals_chips;
 
     //! Whether the data are assigned
     explicit CELER_FUNCTION operator bool() const
     {
-        return scalars && !micro_xs.empty() && !reals.empty();
+        return scalars && !micro_xs_iaea.empty() && !reals_iaea.empty()
+               && !micro_xs_chips.empty() && !reals_chips.empty();
     }
 
     //! Assign from another set of data
@@ -79,8 +76,10 @@ struct GammaNuclearData
     {
         CELER_EXPECT(other);
         scalars = other.scalars;
-        micro_xs = other.micro_xs;
-        reals = other.reals;
+        micro_xs_iaea = other.micro_xs_iaea;
+        reals_iaea = other.reals_iaea;
+        micro_xs_chips = other.micro_xs_chips;
+        reals_chips = other.reals_chips;
 
         return *this;
     }

src/celeritas/em/model/GammaNuclearModel.cc

@@ -7,6 +7,7 @@
 #include "GammaNuclearModel.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 +33,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 +46,24 @@ 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_micro_xs_iaea{&data.reals_iaea,
+                                                &data.micro_xs_iaea};
+    NonuniformGridInserter insert_micro_xs_chips{&data.reals_chips,
+                                                 &data.micro_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_micro_xs_iaea(load_data(z));
+
+        // Build element cross sections above the upper bound of G4PARTICLEXS
+        double emin = data.reals_iaea[data.micro_xs_iaea[el_id].grid.back()];
+        insert_micro_xs_chips(this->calc_chips_xs(z, emin, emax));
     }
-    CELER_ASSERT(data.micro_xs.size() == materials.num_elements());
+    CELER_ASSERT(data.micro_xs_iaea.size() == materials.num_elements());
+    CELER_ASSERT(data.micro_xs_iaea.size() == data.micro_xs_chips.size());
 
     // Move to mirrored data, copying to device
     data_ = CollectionMirror<GammaNuclearData>{std::move(data)};
@@ -98,5 +112,64 @@ void GammaNuclearModel::step(CoreParams const&, CoreStateDevice&) const
 }
 #endif
 
+//---------------------------------------------------------------------------//
+/*!
+ * Build CHIPS gamma-nuclear element cross sections using G4GammaNuclearXS
+ * above the upper energy limit of G4PARTICLEXS/gamma (IAEA) data. The cross
+ * sections are derived from the parameterization developed by M. V. Kossov
+ * (CERN/ITEP Moscow) in the high energy region (106 MeV < E < 50 GeV) and
+ * from a Reggeon-based parameterization in the ultra-high-energy region
+ * (E > 50 GeV). G4GammaNuclearXS uses CHIPS (G4PhotoNuclearCrossSection)
+ * above 150 MeV and performs linear interpolation between the upper energy
+ * limit of G4PARTICLEXS/gamma (IAEA) data and 150 MeV.
+ */
+GammaNuclearModel::result_type
+GammaNuclearModel::calc_chips_xs(AtomicNumber atomic_number,
+                                 double emin,
+                                 double emax) const
+{
+    CELER_EXPECT(atomic_number);
+
+    result_type 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_high = 250;
+    size_type nbin_ultra = 50;
+
+    result.x.resize(nbin_high + nbin_ultra);
+    result.y.resize(nbin_high + nbin_ultra);
+
+    // Bin sizes in logarithmic energy scale
+    double const de_high = (std::log(emid) - std::log(emin)) / (nbin_high - 1);
+    double const de_ultra = (std::log(emax) - std::log(emid))
+                            / (nbin_ultra - 1);
+
+    // Tabulate the cross section from emin to emax
+    MmSqMicroXs xs;
+    for (size_type i = 0; i < nbin_high + nbin_ultra; ++i)
+    {
+        double energy
+            = (i < nbin_high)
+                  ? std::exp(std::log(emin) + de_high * i)
+                  : std::exp(std::log(emid) + de_ultra * (i - nbin_high + 1));
+
+        result.x[i] = energy;
+        xs.value()
+            = helper_->GammaNuclearElementXS(energy, atomic_number.get());
+        result.y[i] = native_value_to<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,6 +22,8 @@ namespace celeritas
 class MaterialParams;
 class ParticleParams;
 
+class EmExtraPhysicsHelper;
+
 //---------------------------------------------------------------------------//
 /*!
  * Set up and launch the gamma-nuclear model interaction.
@@ -31,6 +34,10 @@ class GammaNuclearModel final : public Model, public StaticConcreteAction
     //!@{
     using MevEnergy = units::MevEnergy;
     using ReadData = std::function<inp::Grid(AtomicNumber)>;
+    using result_type = inp::Grid;
+    using BarnXs = units::BarnXs;
+    using MmSqMicroXs
+        = Quantity<UnitProduct<units::Millimeter, units::Millimeter>, double>;
     using HostRef = HostCRef<GammaNuclearData>;
     using DeviceRef = DeviceCRef<GammaNuclearData>;
     //!@}
@@ -62,13 +69,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 ////
+
+    result_type
+    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,39 @@ 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;
+    NonuniformGridCalculator::Values values;
 
-    // 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_iaea[data_.micro_xs_iaea[el_id].grid.back()])
+    {
+        CELER_EXPECT(el_id < data_.micro_xs_iaea.size());
+        grid = data_.micro_xs_iaea[el_id];
+        values = data_.reals_iaea;
+    }
+    else
+    {
+        CELER_EXPECT(el_id < data_.micro_xs_chips.size());
+        grid = data_.micro_xs_chips[el_id];
+        values = data_.reals_chips;
+    }
 
     // Calculate micro cross section at the given energy
-    NonuniformGridCalculator calc_xs(grid, shared_.reals);
+    NonuniformGridCalculator calc_xs(grid, values);
     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