diff --git a/doc/_static/zotero.bib b/doc/_static/zotero.bib index f8605f67e7..025b9ff1c6 100644 --- a/doc/_static/zotero.bib +++ b/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} +} diff --git a/src/celeritas/CMakeLists.txt b/src/celeritas/CMakeLists.txt index 24232686dc..d0492a9e6b 100644 --- a/src/celeritas/CMakeLists.txt +++ b/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 diff --git a/src/celeritas/em/data/GammaNuclearData.hh b/src/celeritas/em/data/GammaNuclearData.hh index efd4cfefdb..66b02c7e48 100644 --- a/src/celeritas/em/data/GammaNuclearData.hh +++ b/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 micro_xs; - - // Backend data + // Microscopic cross sections using G4PARTICLEXS/gamma nuclear (IAEA) data + ElementItems xs_iaea; Items reals; + // Microscopic cross sections using parameterized CHIPS data at high energy + ElementItems 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; diff --git a/src/celeritas/em/model/GammaNuclearModel.cc b/src/celeritas/em/model/GammaNuclearModel.cc index b294572bef..a6f4375977 100644 --- a/src/celeritas/em/model/GammaNuclearModel.cc +++ b/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 data; + helper_ = std::make_shared(); + // 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{std::move(data)}; @@ -98,5 +111,49 @@ void GammaNuclearModel::step(CoreParams const&, CoreStateDevice&) const } #endif +//---------------------------------------------------------------------------// +/*! + * Build CHIPS gamma-nuclear element cross sections using G4GammaNuclearXS. + */ +inp::Grid +GammaNuclearModel::calc_chips_xs(AtomicNumber z, double emin, double emax) const +{ + CELER_EXPECT(z); + + 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 = 5e+4; + 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 + EmExtraPhysicsHelper::MmSqXs xs; + for (size_type i = 0; i < nbin_total; ++i) + { + xs = helper_->calc_gamma_nuclear_xs(z, MevEnergy{result.x[i]}); + result.y[i] + = native_value_to(native_value_from(xs)).value(); + } + + return result; +} + //---------------------------------------------------------------------------// } // namespace celeritas diff --git a/src/celeritas/em/model/GammaNuclearModel.hh b/src/celeritas/em/model/GammaNuclearModel.hh index 385d98f066..664f1f4a99 100644 --- a/src/celeritas/em/model/GammaNuclearModel.hh +++ b/src/celeritas/em/model/GammaNuclearModel.hh @@ -7,6 +7,7 @@ #pragma once #include +#include #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,6 +77,7 @@ class GammaNuclearModel final : public Model, public StaticConcreteAction private: //// DATA //// + std::shared_ptr helper_; // Host/device storage and reference CollectionMirror data_; @@ -69,6 +85,11 @@ class GammaNuclearModel final : public Model, public StaticConcreteAction //// TYPES //// using HostXsData = HostVal; + + //// HELPER FUNCTIONS //// + + inp::Grid + calc_chips_xs(AtomicNumber atomic_number, double emin, double emax) const; }; //---------------------------------------------------------------------------// diff --git a/src/celeritas/em/xs/GammaNuclearMicroXsCalculator.hh b/src/celeritas/em/xs/GammaNuclearMicroXsCalculator.hh index ac09d7986f..b3651c080f 100644 --- a/src/celeritas/em/xs/GammaNuclearMicroXsCalculator.hh +++ b/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}; diff --git a/src/celeritas/g4/EmExtraPhysicsHelper.cc b/src/celeritas/g4/EmExtraPhysicsHelper.cc new file mode 100644 index 0000000000..30ae31d30d --- /dev/null +++ b/src/celeritas/g4/EmExtraPhysicsHelper.cc @@ -0,0 +1,36 @@ +//------------------------------- -*- 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 + +namespace celeritas +{ +//---------------------------------------------------------------------------// +/*! + * Construct with Geant4 provided cross section classes + */ +EmExtraPhysicsHelper::EmExtraPhysicsHelper() +{ + gn_xs_ = std::make_shared(); +} + +//---------------------------------------------------------------------------// +/*! + * Calculate the gamma-nuclear element cross section using G4GammaNuclearXS. + */ +auto EmExtraPhysicsHelper::calc_gamma_nuclear_xs(AtomicNumber z, + MevEnergy energy) const + -> MmSqXs +{ + MmSqXs xs; + xs.value() = gn_xs_->ElementCrossSection(energy.value(), z.get()); + return xs; +} + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/celeritas/g4/EmExtraPhysicsHelper.hh b/src/celeritas/g4/EmExtraPhysicsHelper.hh new file mode 100644 index 0000000000..1cbcde97bb --- /dev/null +++ b/src/celeritas/g4/EmExtraPhysicsHelper.hh @@ -0,0 +1,74 @@ +//------------------------------- -*- 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 + +#include "corecel/Config.hh" + +#include "corecel/math/Quantity.hh" +#include "corecel/math/UnitUtils.hh" +#include "celeritas/Quantities.hh" +#include "celeritas/UnitTypes.hh" +#include "celeritas/phys/AtomicNumber.hh" + +class G4GammaNuclearXS; + +namespace celeritas +{ + +//---------------------------------------------------------------------------// +/*! + * A helper class for interfacing with Geant4 cross section calculations and + * other properties. + * + * This class primarily severs as a wrapper around Geant4 cross section + * calculation methods, which are not directly accessible from Celeritas EM + * physics models. Use of this class requires CELERITAS_USE_GEANT4 to be + * enabled. + */ +class EmExtraPhysicsHelper +{ + public: + //!@{ + using MevEnergy = units::MevEnergy; + using MmSqXs + = Quantity, double>; + //!@} + + public: + // Construct EM extra physics helper + EmExtraPhysicsHelper(); + + // Calculate gamma-nuclear element cross section + MmSqXs calc_gamma_nuclear_xs(AtomicNumber z, MevEnergy energy) const; + + private: + //// DATA //// + std::shared_ptr gn_xs_; +}; + +//---------------------------------------------------------------------------// +// INLINE DEFINITIONS +//---------------------------------------------------------------------------// + +#if !CELERITAS_USE_GEANT4 +inline EmExtraPhysicsHelper::EmExtraPhysicsHelper() +{ + CELER_NOT_CONFIGURED("Geant4"); +} + +inline EmExtraPhysicsHelper::MmSqXs +EmExtraPhysicsHelper::calc_gamma_nuclear_xs(AtomicNumber, MevEnergy) const +{ + CELER_ASSERT_UNREACHABLE(); +} + +#endif + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/test/celeritas/em/GammaNuclear.test.cc b/test/celeritas/em/GammaNuclear.test.cc index d986beb8b2..56cf15e109 100644 --- a/test/celeritas/em/GammaNuclear.test.cc +++ b/test/celeritas/em/GammaNuclear.test.cc @@ -84,7 +84,8 @@ TEST_F(GammaNuclearTest, micro_xs) // Check the size of the element cross section data (G4PARTICLEXS4.1) HostCRef shared = model_->host_ref(); - NonuniformGridRecord grid = shared.micro_xs[el_id]; + + NonuniformGridRecord grid = shared.xs_iaea[el_id]; EXPECT_EQ(grid.grid.size(), 260); // Microscopic cross section (units::BarnXs) in [0.5:100.5] (MeV) @@ -107,6 +108,30 @@ TEST_F(GammaNuclearTest, micro_xs) // Check the gamma-nuclear element cross section at the upper bound XsCalculator calc_upper_xs(shared, MevEnergy{130}); EXPECT_SOFT_EQ(calc_upper_xs(el_id).value(), 0.010895100000000003); + + // Calculate the gamma-nuclear cross section at the high energy region + // using parameterizated data + + NonuniformGridRecord grid_high = shared.xs_chips[el_id]; + EXPECT_EQ(grid_high.grid.size(), 300); + + // Expected microscopic cross section (units::BarnXs) in [130:1e+8] (MeV) + std::vector> const energy_xs + = {{130, 0.010895100000000003}, + {140, 0.016056145231123135}, + {150, 0.02121368730750826}, + {200, 0.041931723222538624}, + {1e+3, 0.032829279254133224}, + {5e+3, 0.018822644663262746}, + {5e+4, 0.01448519295151751}, + {1e+6, 0.017122133350371736}, + {1e+8, 0.027254443598797456}}; + + for (auto i : range(energy_xs.size())) + { + XsCalculator calc_micro_xs(shared, MevEnergy{energy_xs[i].first}); + EXPECT_SOFT_EQ(calc_micro_xs(el_id).value(), energy_xs[i].second); + } } TEST_F(GammaNuclearTest, macro_xs) @@ -116,28 +141,31 @@ TEST_F(GammaNuclearTest, macro_xs) auto calc_xs = MacroXsCalculator( model_->host_ref(), material); - // Macroscopic cross section (\f$ cm^{-1} \f$) in [0.5:100.5] (MeV) - std::vector const expected_macro_xs = {0, - 0.67518515551801506, - 0.27924724815369489, - 0.30744953728122743, - 0.32743928018832685, - 0.31806243520165606}; - - real_type energy = 0.5; - real_type const factor = 2e+1; - for (auto i : range(expected_macro_xs.size())) + // Expected macroscopic cross section (\f$ cm^{-1} \f$)} in [0.5:1e+8](MeV) + std::vector> const energy_xs + = {{0.5, 0}, + {20.5, 0.67518515551801506}, + {40.5, 0.27924724815369489}, + {60.5, 0.30744953728122743}, + {80.5, 0.32743928018832685}, + {100.5, 0.31806243520165606}, + {130, 0.27716766602987458}, + {140, 0.067383743323079975}, + {150, 0.086034497899296999}, + {200, 0.17200535827285135}, + {1e+3, 0.1353591424632776}, + {5e+3, 0.077905738172584824}, + {5e+4, 0.060230059626849054}, + {1e+6, 0.071194572007666074}, + {1e+8, 0.11332515683749959}}; + + for (auto i : range(energy_xs.size())) { - EXPECT_SOFT_EQ( - native_value_to(calc_xs(MevEnergy{energy})).value(), - expected_macro_xs[i]); - energy += factor; + EXPECT_SOFT_EQ(native_value_to( + calc_xs(MevEnergy{energy_xs[i].first})) + .value(), + energy_xs[i].second); } - - // Check the gamma-nuclear interaction cross section at the upper bound - EXPECT_SOFT_EQ( - native_value_to(calc_xs(MevEnergy{130})).value(), - 0.27716766602987458); } //---------------------------------------------------------------------------//