From 7cf690238bf41f75348b95463c36e1c818f9cd0c Mon Sep 17 00:00:00 2001 From: Amanda Lund Date: Thu, 5 Sep 2024 16:05:31 +0000 Subject: [PATCH 1/6] Add failing Wentzel distribution test Fails for positrons with `src/celeritas/random/distribution/RejectionSampler.hh:90: celeritas: precondition failed: fmax_ >= f_"` --- test/celeritas/em/CoulombScattering.test.cc | 106 ++++++++++++++------ 1 file changed, 74 insertions(+), 32 deletions(-) diff --git a/test/celeritas/em/CoulombScattering.test.cc b/test/celeritas/em/CoulombScattering.test.cc index f40d04e163..9327a01f4a 100644 --- a/test/celeritas/em/CoulombScattering.test.cc +++ b/test/celeritas/em/CoulombScattering.test.cc @@ -438,46 +438,88 @@ TEST_F(CoulombScatteringTest, distribution) = this->cutoff_params()->get(mat_id_).energy(ParticleId{0}); std::vector avg_angles; + std::vector avg_engine_samples; - for (real_type energy : {1, 50, 100, 200, 1000, 13000}) - { - this->set_inc_particle(pdg::electron(), MevEnergy{energy}); + auto sample = [&](PDGNumber pdg) { + avg_angles.clear(); + avg_engine_samples.clear(); - WentzelHelper helper(this->particle_track(), - material, - isotope.atomic_number(), - wentzel_->host_ref(), - model_->host_ref().ids, - cutoff); - WentzelDistribution sample_angle(wentzel_->host_ref(), - helper, - this->particle_track(), - isotope, - el_id_, - helper.cos_thetamax_nuclear(), - model_->host_ref().cos_thetamax()); + for (real_type energy : {1, 50, 100, 200, 1000, 13000}) + { + this->set_inc_particle(pdg, MevEnergy{energy}); + + WentzelHelper helper(this->particle_track(), + material, + isotope.atomic_number(), + wentzel_->host_ref(), + model_->host_ref().ids, + cutoff); + WentzelDistribution sample_angle(wentzel_->host_ref(), + helper, + this->particle_track(), + isotope, + el_id_, + helper.cos_thetamax_nuclear(), + model_->host_ref().cos_thetamax()); - RandomEngine& rng_engine = this->rng(); + RandomEngine& rng = this->rng(); - real_type avg_angle = 0; + real_type avg_angle = 0; - int const num_samples = 4096; - for ([[maybe_unused]] int i : range(num_samples)) - { - avg_angle += sample_angle(rng_engine); + int const num_samples = 4096; + for ([[maybe_unused]] int i : range(num_samples)) + { + avg_angle += sample_angle(rng); + } + + avg_angle /= num_samples; + avg_angles.push_back(avg_angle); + avg_engine_samples.push_back(real_type(rng.count()) / num_samples); } + }; - avg_angle /= num_samples; - avg_angles.push_back(avg_angle); + { + sample(pdg::electron()); + static double const expected_avg_angles[] = { + 0.99957853627426, + 0.99999954645904, + 0.99999989882947, + 0.99999996985799, + 0.99999999945722, + 0.99999999999487, + }; + static double const expected_avg_engine_samples[] = { + 5.9287109375, + 5.93359375, + 5.9287109375, + 5.943359375, + 5.927734375, + 5.927734375, + }; + EXPECT_VEC_SOFT_EQ(expected_avg_angles, avg_angles); + EXPECT_VEC_SOFT_EQ(expected_avg_engine_samples, avg_engine_samples); + } + { + sample(pdg::positron()); + static double const expected_avg_angles[] = { + 0.99960786019912, + 0.99999969317473, + 0.99999989582094, + 0.99999998024112, + 0.99999999932915, + 0.99999999996876, + }; + static double const expected_avg_engine_samples[] = { + 5.93701171875, + 5.9306640625, + 5.9345703125, + 5.92431640625, + 5.9345703125, + 5.9267578125, + }; + EXPECT_VEC_SOFT_EQ(expected_avg_angles, avg_angles); + EXPECT_VEC_SOFT_EQ(expected_avg_engine_samples, avg_engine_samples); } - - static double const expected_avg_angles[] = {0.99957853627426, - 0.99999954645904, - 0.99999989882947, - 0.99999996985799, - 0.99999999945722, - 0.99999999999487}; - EXPECT_VEC_SOFT_EQ(expected_avg_angles, avg_angles); } //---------------------------------------------------------------------------// From 7b655aa312bd71a924f9ca99ea9c7e1c042393c6 Mon Sep 17 00:00:00 2001 From: Amanda Lund Date: Thu, 5 Sep 2024 16:14:20 +0000 Subject: [PATCH 2/6] Don't use rejection sampler in Wentzel distribution --- src/celeritas/em/distribution/WentzelDistribution.hh | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/celeritas/em/distribution/WentzelDistribution.hh b/src/celeritas/em/distribution/WentzelDistribution.hh index f7090fdec6..8a83afdb4e 100644 --- a/src/celeritas/em/distribution/WentzelDistribution.hh +++ b/src/celeritas/em/distribution/WentzelDistribution.hh @@ -18,7 +18,6 @@ #include "celeritas/mat/IsotopeView.hh" #include "celeritas/phys/ParticleTrackView.hh" #include "celeritas/random/distribution/BernoulliDistribution.hh" -#include "celeritas/random/distribution/RejectionSampler.hh" namespace celeritas { @@ -178,9 +177,10 @@ CELER_FUNCTION real_type WentzelDistribution::operator()(Engine& rng) const std::sqrt(particle_.beta_sq())); real_type xs = mott_xsec(cos_theta) * ipow<2>(this->calculate_form_factor(cos_theta)); - if (RejectionSampler(xs, helper_.mott_factor())(rng)) + if (xs < helper_.mott_factor() * generate_canonical(rng)) { - // Reject scattering event: no change in direction + // Reject scattering event: no change in direction. Note that the + // cross section can be larger than the Mott factor cos_theta = 1; } } From e35e9b6ece3e1f94f4a1813382600ac887f3353f Mon Sep 17 00:00:00 2001 From: Amanda Lund Date: Thu, 5 Sep 2024 16:27:10 +0000 Subject: [PATCH 3/6] Fix Coulomb scattering construction in celer em physics list - Construct proton when Coulomb scattering is enabled: `Geant4 error: G4ParticleDefinition should be created in PreInit state G4ParticleDefintion::G4ParticleDefinition(): 'PART101' failed` - Always construct Coulomb scattering model with `isCombined=true` (cross sections are zero otherwise) - Set polar angle limit between SS and MSC to zero if single Coulomb scattering is used --- src/celeritas/em/model/CoulombScatteringModel.cc | 4 ++++ src/celeritas/ext/detail/CelerEmStandardPhysics.cc | 9 ++++++--- 2 files changed, 10 insertions(+), 3 deletions(-) diff --git a/src/celeritas/em/model/CoulombScatteringModel.cc b/src/celeritas/em/model/CoulombScatteringModel.cc index 544454c190..e58a4e1624 100644 --- a/src/celeritas/em/model/CoulombScatteringModel.cc +++ b/src/celeritas/em/model/CoulombScatteringModel.cc @@ -60,6 +60,10 @@ CoulombScatteringModel::CoulombScatteringModel(ActionId id, = imported_.energy_grid_bounds(data_.ids.electron, MaterialId{0}); // Check that the bounds are the same for all particles/materials + // \todo This is only expected when using Coulomb scattering with the + // Wentzel VI model above the MSC energy limit. When the MSC energy limit + // is not set, the model energy grid bounds are material dependent and + // require material-dependent applicability for (auto pid : {data_.ids.electron, data_.ids.positron}) { for (auto mid : range(MaterialId{materials.num_materials()})) diff --git a/src/celeritas/ext/detail/CelerEmStandardPhysics.cc b/src/celeritas/ext/detail/CelerEmStandardPhysics.cc index 5038ebd257..e6359eb6eb 100644 --- a/src/celeritas/ext/detail/CelerEmStandardPhysics.cc +++ b/src/celeritas/ext/detail/CelerEmStandardPhysics.cc @@ -155,7 +155,7 @@ void CelerEmStandardPhysics::ConstructParticle() G4Gamma::GammaDefinition(); G4Electron::ElectronDefinition(); G4Positron::PositronDefinition(); - if (options_.msc != MscModelSelection::none) + if (options_.msc != MscModelSelection::none || options_.coulomb_scattering) { G4Proton::ProtonDefinition(); } @@ -371,14 +371,17 @@ void CelerEmStandardPhysics::add_e_processes(G4ParticleDefinition* p) else { auto process = std::make_unique(); - auto model = std::make_unique( - /* isCombined = */ options_.msc != MMS::none); + auto model = std::make_unique(); if (set_energy_limit) { process->SetMinKinEnergy(msc_energy_limit); model->SetLowEnergyLimit(msc_energy_limit); model->SetActivationLowEnergyLimit(msc_energy_limit); } + if (options_.msc == MMS::none) + { + G4EmParameters::Instance()->SetMscThetaLimit(0); + } CELER_LOG(debug) << "Loaded single Coulomb scattering with " "G4eCoulombScatteringModel from " From 8e0fea7c3bfffd15d0ec93a00fbdac8d6cf13a9e Mon Sep 17 00:00:00 2001 From: Amanda Lund Date: Sun, 8 Sep 2024 18:46:35 +0000 Subject: [PATCH 4/6] Revert "Add failing Wentzel distribution test" This reverts commit 7cf690238bf41f75348b95463c36e1c818f9cd0c. --- test/celeritas/em/CoulombScattering.test.cc | 106 ++++++-------------- 1 file changed, 32 insertions(+), 74 deletions(-) diff --git a/test/celeritas/em/CoulombScattering.test.cc b/test/celeritas/em/CoulombScattering.test.cc index 9327a01f4a..f40d04e163 100644 --- a/test/celeritas/em/CoulombScattering.test.cc +++ b/test/celeritas/em/CoulombScattering.test.cc @@ -438,88 +438,46 @@ TEST_F(CoulombScatteringTest, distribution) = this->cutoff_params()->get(mat_id_).energy(ParticleId{0}); std::vector avg_angles; - std::vector avg_engine_samples; - auto sample = [&](PDGNumber pdg) { - avg_angles.clear(); - avg_engine_samples.clear(); - - for (real_type energy : {1, 50, 100, 200, 1000, 13000}) - { - this->set_inc_particle(pdg, MevEnergy{energy}); - - WentzelHelper helper(this->particle_track(), - material, - isotope.atomic_number(), - wentzel_->host_ref(), - model_->host_ref().ids, - cutoff); - WentzelDistribution sample_angle(wentzel_->host_ref(), - helper, - this->particle_track(), - isotope, - el_id_, - helper.cos_thetamax_nuclear(), - model_->host_ref().cos_thetamax()); + for (real_type energy : {1, 50, 100, 200, 1000, 13000}) + { + this->set_inc_particle(pdg::electron(), MevEnergy{energy}); - RandomEngine& rng = this->rng(); + WentzelHelper helper(this->particle_track(), + material, + isotope.atomic_number(), + wentzel_->host_ref(), + model_->host_ref().ids, + cutoff); + WentzelDistribution sample_angle(wentzel_->host_ref(), + helper, + this->particle_track(), + isotope, + el_id_, + helper.cos_thetamax_nuclear(), + model_->host_ref().cos_thetamax()); - real_type avg_angle = 0; + RandomEngine& rng_engine = this->rng(); - int const num_samples = 4096; - for ([[maybe_unused]] int i : range(num_samples)) - { - avg_angle += sample_angle(rng); - } + real_type avg_angle = 0; - avg_angle /= num_samples; - avg_angles.push_back(avg_angle); - avg_engine_samples.push_back(real_type(rng.count()) / num_samples); + int const num_samples = 4096; + for ([[maybe_unused]] int i : range(num_samples)) + { + avg_angle += sample_angle(rng_engine); } - }; - { - sample(pdg::electron()); - static double const expected_avg_angles[] = { - 0.99957853627426, - 0.99999954645904, - 0.99999989882947, - 0.99999996985799, - 0.99999999945722, - 0.99999999999487, - }; - static double const expected_avg_engine_samples[] = { - 5.9287109375, - 5.93359375, - 5.9287109375, - 5.943359375, - 5.927734375, - 5.927734375, - }; - EXPECT_VEC_SOFT_EQ(expected_avg_angles, avg_angles); - EXPECT_VEC_SOFT_EQ(expected_avg_engine_samples, avg_engine_samples); - } - { - sample(pdg::positron()); - static double const expected_avg_angles[] = { - 0.99960786019912, - 0.99999969317473, - 0.99999989582094, - 0.99999998024112, - 0.99999999932915, - 0.99999999996876, - }; - static double const expected_avg_engine_samples[] = { - 5.93701171875, - 5.9306640625, - 5.9345703125, - 5.92431640625, - 5.9345703125, - 5.9267578125, - }; - EXPECT_VEC_SOFT_EQ(expected_avg_angles, avg_angles); - EXPECT_VEC_SOFT_EQ(expected_avg_engine_samples, avg_engine_samples); + avg_angle /= num_samples; + avg_angles.push_back(avg_angle); } + + static double const expected_avg_angles[] = {0.99957853627426, + 0.99999954645904, + 0.99999989882947, + 0.99999996985799, + 0.99999999945722, + 0.99999999999487}; + EXPECT_VEC_SOFT_EQ(expected_avg_angles, avg_angles); } //---------------------------------------------------------------------------// From 8e75382e4fda468ac34125b873a2fbc2726d235b Mon Sep 17 00:00:00 2001 From: Amanda Lund Date: Sun, 8 Sep 2024 18:46:55 +0000 Subject: [PATCH 5/6] Revert "Don't use rejection sampler in Wentzel distribution" This reverts commit 7b655aa312bd71a924f9ca99ea9c7e1c042393c6. --- src/celeritas/em/distribution/WentzelDistribution.hh | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/celeritas/em/distribution/WentzelDistribution.hh b/src/celeritas/em/distribution/WentzelDistribution.hh index 8a83afdb4e..f7090fdec6 100644 --- a/src/celeritas/em/distribution/WentzelDistribution.hh +++ b/src/celeritas/em/distribution/WentzelDistribution.hh @@ -18,6 +18,7 @@ #include "celeritas/mat/IsotopeView.hh" #include "celeritas/phys/ParticleTrackView.hh" #include "celeritas/random/distribution/BernoulliDistribution.hh" +#include "celeritas/random/distribution/RejectionSampler.hh" namespace celeritas { @@ -177,10 +178,9 @@ CELER_FUNCTION real_type WentzelDistribution::operator()(Engine& rng) const std::sqrt(particle_.beta_sq())); real_type xs = mott_xsec(cos_theta) * ipow<2>(this->calculate_form_factor(cos_theta)); - if (xs < helper_.mott_factor() * generate_canonical(rng)) + if (RejectionSampler(xs, helper_.mott_factor())(rng)) { - // Reject scattering event: no change in direction. Note that the - // cross section can be larger than the Mott factor + // Reject scattering event: no change in direction cos_theta = 1; } } From 8997fbed539ec15abcd02b5eb88c53d4830dc2d4 Mon Sep 17 00:00:00 2001 From: Amanda Lund Date: Sun, 8 Sep 2024 18:48:40 +0000 Subject: [PATCH 6/6] Address feedback --- src/celeritas/em/model/CoulombScatteringModel.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/celeritas/em/model/CoulombScatteringModel.cc b/src/celeritas/em/model/CoulombScatteringModel.cc index e58a4e1624..d36d2003a4 100644 --- a/src/celeritas/em/model/CoulombScatteringModel.cc +++ b/src/celeritas/em/model/CoulombScatteringModel.cc @@ -60,7 +60,7 @@ CoulombScatteringModel::CoulombScatteringModel(ActionId id, = imported_.energy_grid_bounds(data_.ids.electron, MaterialId{0}); // Check that the bounds are the same for all particles/materials - // \todo This is only expected when using Coulomb scattering with the + // TODO: This is only expected when using Coulomb scattering with the // Wentzel VI model above the MSC energy limit. When the MSC energy limit // is not set, the model energy grid bounds are material dependent and // require material-dependent applicability