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); 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 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); diff --git a/src/celeritas/ext/detail/GeantPhysicsList.cc b/src/celeritas/ext/detail/GeantPhysicsList.cc index 93e50fe9d1..7d10573e30 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 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"; } 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); 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; }; //---------------------------------------------------------------------------// 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); } }