Support Geant4 11.1+

app/demo-geant-integration/demo-geant-integration.cc

@@ -69,7 +69,15 @@ void run(std::string const& macro_filename)
     run_manager->SetUserInitialization(new FTFP_BERT{/* verbosity = */ 0});
     run_manager->SetUserInitialization(new demo_geant::ActionInitialization());
 
-    demo_geant::GlobalSetup::Instance()->SetIgnoreProcesses({"CoulombScat"});
+    std::vector<std::string> ignore_processes = {"CoulombScat"};
+    if (G4VERSION_NUMBER >= 1110)
+    {
+        CELER_LOG(warning) << "Default Rayleigh scattering 'MinKinEnergyPrim' "
+                              "is not compatible between Celeritas and "
+                              "Geant4@11.1: disabling Rayleigh scattering";
+        ignore_processes.push_back("Rayl");
+    }
+    demo_geant::GlobalSetup::Instance()->SetIgnoreProcesses(ignore_processes);
 
     G4UImanager* ui = G4UImanager::GetUIpointer();
     CELER_ASSERT(ui);

src/celeritas/ext/GeantImporter.cc

@@ -30,6 +30,7 @@
 #include <G4ParticleTable.hh>
 #include <G4ProcessManager.hh>
 #include <G4ProcessType.hh>
+#include <G4PropagatorInField.hh>
 #include <G4ProcessVector.hh>
 #include <G4ProductionCuts.hh>
 #include <G4ProductionCutsTable.hh>

src/celeritas/ext/detail/GeantLoggerAdapter.cc

@@ -9,7 +9,12 @@
 
 #include <cctype>
 #include <memory>
-#include <G4strstreambuf.hh>
+#include <G4Version.hh>
+#if CELER_G4SSBUF
+#    include <G4strstreambuf.hh>
+#else
+#    include <G4ios.hh>
+#endif
 
 #include "corecel/io/Logger.hh"
 #include "corecel/io/StringUtils.hh"
@@ -23,19 +28,27 @@ namespace detail
  * Redirect geant4's stdout/cerr on construction.
  */
 GeantLoggerAdapter::GeantLoggerAdapter()
+#if CELER_G4SSBUF
     : saved_cout_(G4coutbuf.GetDestination())
     , saved_cerr_(G4cerrbuf.GetDestination())
 {
     G4coutbuf.SetDestination(this);
     G4cerrbuf.SetDestination(this);
 }
+#else
+{
+    // See Geant4 change global-V11-01-01
+    G4iosSetDestination(this);
+}
+#endif
 
 //---------------------------------------------------------------------------//
 /*!
  * Restore iostream buffers on destruction.
  */
 GeantLoggerAdapter::~GeantLoggerAdapter()
 {
+#if CELER_G4SSBUF
     if (G4coutbuf.GetDestination() == this)
     {
         G4coutbuf.SetDestination(saved_cout_);
@@ -44,6 +57,9 @@ GeantLoggerAdapter::~GeantLoggerAdapter()
     {
         G4cerrbuf.SetDestination(saved_cerr_);
     }
+#else
+    G4iosSetDestination(nullptr);
+#endif
 }
 
 //---------------------------------------------------------------------------//

src/celeritas/ext/detail/GeantLoggerAdapter.hh

@@ -9,7 +9,14 @@
 
 #include <G4String.hh>
 #include <G4Types.hh>
+#include <G4Version.hh>
 #include <G4coutDestination.hh>
+#if G4VERSION_NUMBER > 1111
+// No ability to include G4strstreambuf
+#    define CELER_G4SSBUF 0
+#else
+#    define CELER_G4SSBUF 1
+#endif
 
 #include "corecel/io/LoggerTypes.hh"
 
@@ -35,8 +42,10 @@ class GeantLoggerAdapter : public G4coutDestination
   private:
     //// DATA ////
 
+#if CELER_G4SSBUF
     G4coutDestination* saved_cout_;
     G4coutDestination* saved_cerr_;
+#endif
 
     //// IMPLEMENTATION ////
     G4int log_impl(G4String const& str, LogLevel level);

src/celeritas/ext/detail/GeantPhysicsList.cc

@@ -28,6 +28,7 @@
 #include <G4Proton.hh>
 #include <G4RayleighScattering.hh>
 #include <G4UrbanMscModel.hh>
+#include <G4Version.hh>
 #include <G4WentzelVIModel.hh>
 #include <G4eCoulombScatteringModel.hh>
 #include <G4eIonisation.hh>
@@ -188,7 +189,14 @@ void GeantPhysicsList::add_gamma_processes()
     if (options_.rayleigh_scattering)
     {
         // Rayleigh: G4LivermoreRayleighModel
-        add_process(new G4RayleighScattering());
+        auto rayl = std::make_unique<G4RayleighScattering>();
+        if (G4VERSION_NUMBER >= 1110)
+        {
+            // Revert MinKinEnergyPrim (150 keV) to its Geant 11.0 value so that the lower
+            // and energy grids have the same log spacing
+            rayl->SetMinKinEnergyPrim(100 * CLHEP::keV);
+        }
+        add_process(rayl.release());
         CELER_LOG(debug) << "Loaded Rayleigh scattering with "
                             "G4LivermoreRayleighModel";
     }

src/celeritas/grid/ValueGridBuilder.cc

@@ -51,11 +51,6 @@ bool has_log_spacing(SpanConstReal vec)
     return true;
 }
 
-bool has_same_log_spacing(SpanConstReal first, SpanConstReal second)
-{
-    return soft_equal(calc_log_delta(first), calc_log_delta(second));
-}
-
 bool is_nonnegative(SpanConstReal vec)
 {
     return std::all_of(
@@ -109,13 +104,20 @@ ValueGridXsBuilder::from_geant(SpanConstReal lambda_energy,
     CELER_EXPECT(is_contiguous_increasing(lambda_energy, lambda_prim_energy));
     CELER_EXPECT(has_log_spacing(lambda_energy)
                  && has_log_spacing(lambda_prim_energy));
-    CELER_EXPECT(has_same_log_spacing(lambda_energy, lambda_prim_energy));
     CELER_EXPECT(lambda.size() == lambda_energy.size());
     CELER_EXPECT(lambda_prim.size() == lambda_prim_energy.size());
     CELER_EXPECT(soft_equal(lambda.back(),
                             lambda_prim.front() / lambda_prim_energy.front()));
     CELER_EXPECT(is_nonnegative(lambda) && is_nonnegative(lambda_prim));
 
+    real_type const log_delta_lo = calc_log_delta(lambda_energy);
+    real_type const log_delta_hi = calc_log_delta(lambda_prim_energy);
+    CELER_VALIDATE(
+        soft_equal(log_delta_lo, log_delta_hi),
+        << "Lower and upper energy grids have inconsistent spacing: "
+           "log delta E for lower grid is "
+        << log_delta_lo << " log(MeV) per bin but upper is " << log_delta_hi);
+
     // Concatenate the two XS vectors: insert the scaled (lambda_prim) value at
     // the coincident point.
     VecReal xs(lambda.size() + lambda_prim.size() - 1);

src/celeritas/phys/ImportedProcessAdapter.cc

@@ -8,6 +8,7 @@
 #include "ImportedProcessAdapter.hh"
 
 #include <algorithm>
+#include <exception>
 #include <tuple>
 #include <type_traits>
 
@@ -26,6 +27,42 @@
 
 namespace celeritas
 {
+namespace
+{
+//---------------------------------------------------------------------------//
+//! Rough helper class to hopefully help a little with debugging errors
+class IPAContextException : public RichContextException
+{
+  public:
+    IPAContextException(ParticleId id, ImportProcessClass ipc, MaterialId mid);
+
+    //! This class type
+    char const* type() const final { return "ImportProcessAdapterContext"; }
+
+    // Save context to a JSON object
+    void output(JsonPimpl*) const final {}
+
+    //! Get an explanatory message
+    char const* what() const noexcept final { return what_.c_str(); }
+
+  private:
+    std::string what_;
+};
+
+//---------------------------------------------------------------------------//
+IPAContextException::IPAContextException(ParticleId id,
+                                         ImportProcessClass ipc,
+                                         MaterialId mid)
+{
+    std::stringstream os;
+    os << "Particle ID=" << id.unchecked_get() << ", process '"
+       << to_cstring(ipc) << ", material ID=" << mid.unchecked_get();
+    what_ = os.str();
+}
+
+//---------------------------------------------------------------------------//
+}  // namespace
+
 //---------------------------------------------------------------------------//
 /*!
  * Construct with imported data.
@@ -106,7 +143,7 @@ ImportedProcessAdapter::ImportedProcessAdapter(SPConstImported imported,
                                                SPConstParticles const& particles,
                                                ImportProcessClass process_class,
                                                SpanConstPDG pdg_numbers)
-    : imported_(std::move(imported))
+    : imported_(std::move(imported)), process_class_(process_class)
 {
     CELER_EXPECT(particles);
     CELER_EXPECT(!pdg_numbers.empty());
@@ -186,7 +223,25 @@ ImportedProcessAdapter::ImportedProcessAdapter(
 /*!
  * Get the interaction cross sections for the given material and particle.
  */
-auto ImportedProcessAdapter::step_limits(Applicability applic) const
+auto ImportedProcessAdapter::step_limits(Applicability const& applic) const
+    -> StepLimitBuilders
+{
+    try
+    {
+        return this->step_limits_impl(std::move(applic));
+    }
+    catch (...)
+    {
+        std::throw_with_nested(IPAContextException(
+            applic.particle, process_class_, applic.material));
+    }
+}
+
+//---------------------------------------------------------------------------//
+/*!
+ * Get the interaction cross sections for the given material and particle.
+ */
+auto ImportedProcessAdapter::step_limits_impl(Applicability const& applic) const
     -> StepLimitBuilders
 {
     CELER_EXPECT(ids_.count(applic.particle));

src/celeritas/phys/ImportedProcessAdapter.hh

@@ -95,7 +95,7 @@ class ImportedProcessAdapter
                            std::initializer_list<PDGNumber> pdg_numbers);
 
     // Construct step limits from the given particle/material type
-    StepLimitBuilders step_limits(Applicability applic) const;
+    StepLimitBuilders step_limits(Applicability const& applic) const;
 
     // Get the lambda table for the given particle ID
     inline ImportPhysicsTable const& get_lambda(ParticleId id) const;
@@ -117,7 +117,11 @@ class ImportedProcessAdapter
     };
 
     SPConstImported imported_;
+    ImportProcessClass process_class_;
     std::map<ParticleId, ParticleProcessIds> ids_;
+
+    // Construct step limits from the given particle/material type
+    StepLimitBuilders step_limits_impl(Applicability const& applic) const;
 };
 
 //---------------------------------------------------------------------------//

test/celeritas/global/AlongStep.test.cc

@@ -488,7 +488,7 @@ TEST_F(SimpleCmsRZFieldAlongStepTest, msc_rzfield)
 
         auto result = this->run(inp, num_tracks);
         EXPECT_SOFT_EQ(4.1632771293464517, result.displacement);
-        EXPECT_SOFT_EQ(-0.59445466152831616, result.angle);
+        EXPECT_SOFT_NEAR(-0.59445466152831616, result.angle, 2e-12);
     }
 }