From 75f5e5cf9a4ff16fdd184093a566392ab7ff6e18 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Mon, 22 May 2023 10:06:04 -0400 Subject: [PATCH 01/12] Add header for forward-declared class --- src/celeritas/ext/GeantImporter.cc | 1 + 1 file changed, 1 insertion(+) diff --git a/src/celeritas/ext/GeantImporter.cc b/src/celeritas/ext/GeantImporter.cc index 8908e9a695..5270ab3902 100644 --- a/src/celeritas/ext/GeantImporter.cc +++ b/src/celeritas/ext/GeantImporter.cc @@ -30,6 +30,7 @@ #include #include #include +#include #include #include #include From 3d6dbaaf129d05e316bf8727d12feb37a708afd4 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Mon, 22 May 2023 10:13:03 -0400 Subject: [PATCH 02/12] Fudge lower cross section value if below energy cutoff --- src/celeritas/grid/ValueGridBuilder.cc | 31 ++++++++++++++++++++------ 1 file changed, 24 insertions(+), 7 deletions(-) diff --git a/src/celeritas/grid/ValueGridBuilder.cc b/src/celeritas/grid/ValueGridBuilder.cc index a149a2ac6c..2af1361c28 100644 --- a/src/celeritas/grid/ValueGridBuilder.cc +++ b/src/celeritas/grid/ValueGridBuilder.cc @@ -12,6 +12,7 @@ #include #include "corecel/Types.hh" +#include "corecel/io/Logger.hh" #include "corecel/cont/Range.hh" #include "corecel/grid/UniformGrid.hh" #include "corecel/grid/UniformGridData.hh" @@ -51,11 +52,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,7 +105,6 @@ 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(), @@ -123,8 +118,30 @@ ValueGridXsBuilder::from_geant(SpanConstReal lambda_energy, dst = std::copy(lambda_prim.begin(), lambda_prim.end(), dst); CELER_ASSERT(dst == xs.end()); + + real_type const emin = [&] { + real_type const log_delta_lo = calc_log_delta(lambda_energy); + real_type const log_delta_hi = calc_log_delta(lambda_prim_energy); + + if (soft_equal(log_delta_lo, log_delta_hi)) + { + // Always the case for Geant4 < 11.1 + return lambda_energy.front(); + } + real_type const adjusted + = lambda_energy.back() + / std::pow(log_delta_hi, lambda_energy.size() - 1); + CELER_LOG(warning) + << "Adjusting cross section energy lower limit from " + << lambda_energy.front() << " to " << adjusted + << " [MeV] due to unexpected grid spacing"; + + // Reset lower energy point based on high energy spacing + return adjusted; + }(); + // Construct the grid - return std::make_unique(lambda_energy.front(), + return std::make_unique(emin, lambda_prim_energy.front(), lambda_prim_energy.back(), VecReal(std::move(xs))); From ce3ad5e0e1f04f3344c2c1b2fdcaa3a9b461b8c2 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Mon, 22 May 2023 10:14:02 -0400 Subject: [PATCH 03/12] Loosen tolerance presumably due to slightly different cross sections --- test/celeritas/global/AlongStep.test.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/celeritas/global/AlongStep.test.cc b/test/celeritas/global/AlongStep.test.cc index f732710e4f..c3820852b0 100644 --- a/test/celeritas/global/AlongStep.test.cc +++ b/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); } } From 722ee3b3b7bf9df5b8d63bce34e156b41b68f7b7 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Mon, 22 May 2023 10:16:33 -0400 Subject: [PATCH 04/12] Adjust logger to work with newer Geant4 interface (global-V11-01-01) --- src/celeritas/ext/detail/GeantLoggerAdapter.cc | 18 +++++++++++++++++- src/celeritas/ext/detail/GeantLoggerAdapter.hh | 9 +++++++++ 2 files changed, 26 insertions(+), 1 deletion(-) diff --git a/src/celeritas/ext/detail/GeantLoggerAdapter.cc b/src/celeritas/ext/detail/GeantLoggerAdapter.cc index 8900355703..3d2d898bf7 100644 --- a/src/celeritas/ext/detail/GeantLoggerAdapter.cc +++ b/src/celeritas/ext/detail/GeantLoggerAdapter.cc @@ -9,7 +9,12 @@ #include #include -#include +#include +#if CELER_G4SSBUF +# include +#else +# include +#endif #include "corecel/io/Logger.hh" #include "corecel/io/StringUtils.hh" @@ -23,12 +28,19 @@ 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 //---------------------------------------------------------------------------// /*! @@ -36,6 +48,7 @@ GeantLoggerAdapter::GeantLoggerAdapter() */ 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 } //---------------------------------------------------------------------------// diff --git a/src/celeritas/ext/detail/GeantLoggerAdapter.hh b/src/celeritas/ext/detail/GeantLoggerAdapter.hh index 0e0d36cff8..4e60e6b020 100644 --- a/src/celeritas/ext/detail/GeantLoggerAdapter.hh +++ b/src/celeritas/ext/detail/GeantLoggerAdapter.hh @@ -9,7 +9,14 @@ #include #include +#include #include +#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); From 09ffeb5606ce465572d306263fdd28be92bdd66a Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Mon, 22 May 2023 11:51:27 -0400 Subject: [PATCH 05/12] REVERTME: work with geant4-dev --- src/celeritas/ext/detail/GeantLoggerAdapter.hh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/celeritas/ext/detail/GeantLoggerAdapter.hh b/src/celeritas/ext/detail/GeantLoggerAdapter.hh index 4e60e6b020..5d8f8ec82e 100644 --- a/src/celeritas/ext/detail/GeantLoggerAdapter.hh +++ b/src/celeritas/ext/detail/GeantLoggerAdapter.hh @@ -11,7 +11,7 @@ #include #include #include -#if G4VERSION_NUMBER > 1111 +#if G4VERSION_NUMBER >= 1111 // No ability to include G4strstreambuf # define CELER_G4SSBUF 0 #else From 007186bdaea875a3a009575d57c7e1c8b0867cda Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Mon, 22 May 2023 17:32:17 -0400 Subject: [PATCH 06/12] Revert "Fudge lower cross section value if below energy cutoff" This reverts commit 3d6dbaaf129d05e316bf8727d12feb37a708afd4. --- src/celeritas/grid/ValueGridBuilder.cc | 31 ++++++-------------------- 1 file changed, 7 insertions(+), 24 deletions(-) diff --git a/src/celeritas/grid/ValueGridBuilder.cc b/src/celeritas/grid/ValueGridBuilder.cc index 2af1361c28..a149a2ac6c 100644 --- a/src/celeritas/grid/ValueGridBuilder.cc +++ b/src/celeritas/grid/ValueGridBuilder.cc @@ -12,7 +12,6 @@ #include #include "corecel/Types.hh" -#include "corecel/io/Logger.hh" #include "corecel/cont/Range.hh" #include "corecel/grid/UniformGrid.hh" #include "corecel/grid/UniformGridData.hh" @@ -52,6 +51,11 @@ 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( @@ -105,6 +109,7 @@ 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(), @@ -118,30 +123,8 @@ ValueGridXsBuilder::from_geant(SpanConstReal lambda_energy, dst = std::copy(lambda_prim.begin(), lambda_prim.end(), dst); CELER_ASSERT(dst == xs.end()); - - real_type const emin = [&] { - real_type const log_delta_lo = calc_log_delta(lambda_energy); - real_type const log_delta_hi = calc_log_delta(lambda_prim_energy); - - if (soft_equal(log_delta_lo, log_delta_hi)) - { - // Always the case for Geant4 < 11.1 - return lambda_energy.front(); - } - real_type const adjusted - = lambda_energy.back() - / std::pow(log_delta_hi, lambda_energy.size() - 1); - CELER_LOG(warning) - << "Adjusting cross section energy lower limit from " - << lambda_energy.front() << " to " << adjusted - << " [MeV] due to unexpected grid spacing"; - - // Reset lower energy point based on high energy spacing - return adjusted; - }(); - // Construct the grid - return std::make_unique(emin, + return std::make_unique(lambda_energy.front(), lambda_prim_energy.front(), lambda_prim_energy.back(), VecReal(std::move(xs))); From a910dc0ce0295c417c5da9f5ae729ed9b114ee11 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Mon, 22 May 2023 18:01:23 -0400 Subject: [PATCH 07/12] Override new default Rayleigh 'MinKinEnergyPrim' when we create the process --- src/celeritas/ext/detail/GeantPhysicsList.cc | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/src/celeritas/ext/detail/GeantPhysicsList.cc b/src/celeritas/ext/detail/GeantPhysicsList.cc index 93e50fe9d1..6f252846b0 100644 --- a/src/celeritas/ext/detail/GeantPhysicsList.cc +++ b/src/celeritas/ext/detail/GeantPhysicsList.cc @@ -28,6 +28,7 @@ #include #include #include +#include #include #include #include @@ -188,7 +189,14 @@ void GeantPhysicsList::add_gamma_processes() if (options_.rayleigh_scattering) { // Rayleigh: G4LivermoreRayleighModel - add_process(new G4RayleighScattering()); + auto rayl = std::make_unique(); + if (G4VERSION_NUMBER >= 1110) + { + // Revert MinKinEnergy 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"; } From fa2189b21f18b424f992d93d284a3672b369b189 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Mon, 22 May 2023 18:01:47 -0400 Subject: [PATCH 08/12] Add Rayleigh to the ignored list for demo-geant-integration for now --- app/demo-geant-integration/demo-geant-integration.cc | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/app/demo-geant-integration/demo-geant-integration.cc b/app/demo-geant-integration/demo-geant-integration.cc index e4ad7dff58..b2c78875fc 100644 --- a/app/demo-geant-integration/demo-geant-integration.cc +++ b/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 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); From 665f1e260f3269fa00b31cf33657fb403f7b9479 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Mon, 22 May 2023 22:00:12 -0400 Subject: [PATCH 09/12] Always validate the grid spacing --- src/celeritas/grid/ValueGridBuilder.cc | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/src/celeritas/grid/ValueGridBuilder.cc b/src/celeritas/grid/ValueGridBuilder.cc index a149a2ac6c..b38c6906a8 100644 --- a/src/celeritas/grid/ValueGridBuilder.cc +++ b/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); From fa3bc3135727b2bdd174292cf75a70bf98a15e8a Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Mon, 22 May 2023 22:00:55 -0400 Subject: [PATCH 10/12] Add a context exception to describe failures during import ``` ExceptionConverter.cc:157: critical: The following error is from: Particle ID=0, process 'rayleigh, material ID=0 -------- EEEE ------- G4Exception-START -------- EEEE ------- *** G4Exception : celer0001 issued by : celeritas/grid/ValueGridBuilder.cc:119 Lower and upper energy grids have inconsistent spacing: log delta E for lower grid is 1.39434 log(MeV) per bin but upper is 1.38778 *** Fatal Exception *** core dump *** ```` --- src/celeritas/phys/ImportedProcessAdapter.cc | 59 +++++++++++++++++++- src/celeritas/phys/ImportedProcessAdapter.hh | 6 +- 2 files changed, 62 insertions(+), 3 deletions(-) diff --git a/src/celeritas/phys/ImportedProcessAdapter.cc b/src/celeritas/phys/ImportedProcessAdapter.cc index 0f54cfb9be..9db907a571 100644 --- a/src/celeritas/phys/ImportedProcessAdapter.cc +++ b/src/celeritas/phys/ImportedProcessAdapter.cc @@ -8,6 +8,7 @@ #include "ImportedProcessAdapter.hh" #include +#include #include #include @@ -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)); diff --git a/src/celeritas/phys/ImportedProcessAdapter.hh b/src/celeritas/phys/ImportedProcessAdapter.hh index 5689f4d091..132e73722d 100644 --- a/src/celeritas/phys/ImportedProcessAdapter.hh +++ b/src/celeritas/phys/ImportedProcessAdapter.hh @@ -95,7 +95,7 @@ class ImportedProcessAdapter std::initializer_list 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 ids_; + + // Construct step limits from the given particle/material type + StepLimitBuilders step_limits_impl(Applicability const& applic) const; }; //---------------------------------------------------------------------------// From 060fa54f29b59bf33b297e7bf943c45bdd624ec8 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Mon, 22 May 2023 22:02:08 -0400 Subject: [PATCH 11/12] Revert "REVERTME: work with geant4-dev" This reverts commit 09ffeb5606ce465572d306263fdd28be92bdd66a. --- src/celeritas/ext/detail/GeantLoggerAdapter.hh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/celeritas/ext/detail/GeantLoggerAdapter.hh b/src/celeritas/ext/detail/GeantLoggerAdapter.hh index 5d8f8ec82e..4e60e6b020 100644 --- a/src/celeritas/ext/detail/GeantLoggerAdapter.hh +++ b/src/celeritas/ext/detail/GeantLoggerAdapter.hh @@ -11,7 +11,7 @@ #include #include #include -#if G4VERSION_NUMBER >= 1111 +#if G4VERSION_NUMBER > 1111 // No ability to include G4strstreambuf # define CELER_G4SSBUF 0 #else From ddcac154d62e283a5687863efbf1233c233f04d1 Mon Sep 17 00:00:00 2001 From: "Seth R. Johnson" Date: Tue, 23 May 2023 11:02:52 -0400 Subject: [PATCH 12/12] Update src/celeritas/ext/detail/GeantPhysicsList.cc Co-authored-by: Soon Yung Jun --- src/celeritas/ext/detail/GeantPhysicsList.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/celeritas/ext/detail/GeantPhysicsList.cc b/src/celeritas/ext/detail/GeantPhysicsList.cc index 6f252846b0..7d10573e30 100644 --- a/src/celeritas/ext/detail/GeantPhysicsList.cc +++ b/src/celeritas/ext/detail/GeantPhysicsList.cc @@ -192,7 +192,7 @@ void GeantPhysicsList::add_gamma_processes() auto rayl = std::make_unique(); if (G4VERSION_NUMBER >= 1110) { - // Revert MinKinEnergy to its Geant 11.0 value so that the lower + // 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); }