Skip to content

Commit

Permalink
Use consistent energy intervals for model applicability and selection (…
Browse files Browse the repository at this point in the history 
…celeritas-project#1015)

* Use consistent energy intervals for model applicability and selection

* Add check for SB upper limit and update limit in combined brems interactor

---------

Co-authored-by: Seth R. Johnson <johnsonsr@ornl.gov>
  • Loading branch information
@amandalund @sethrj
amandalund and sethrj committed Nov 20, 2023
1 parent bad52d7 commit 88237d1
Show file tree
Hide file tree
Showing 10 changed files with 52 additions and 74 deletions.
2 changes: 1 addition & 1 deletion src/celeritas/em/interactor/BetheHeitlerInteractor.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -150,7 +150,7 @@ CELER_FUNCTION BetheHeitlerInteractor::BetheHeitlerInteractor(
Energy{inc_energy_})
{
CELER_EXPECT(particle.particle_id() == shared_.ids.gamma);
CELER_EXPECT(inc_energy_ > 2 * value_as<Mass>(shared_.electron_mass));
CELER_EXPECT(inc_energy_ >= 2 * value_as<Mass>(shared_.electron_mass));

epsilon0_ = value_as<Mass>(shared_.electron_mass) / inc_energy_;
}
Expand Down
12 changes: 3 additions & 9 deletions src/celeritas/em/interactor/CombinedBremInteractor.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -119,7 +119,8 @@ CombinedBremInteractor::CombinedBremInteractor(
{
CELER_EXPECT(is_electron_
|| particle.particle_id() == shared.rb_data.ids.positron);
CELER_EXPECT(gamma_cutoff_.value() > 0);
CELER_EXPECT(gamma_cutoff_ > zero_quantity());
CELER_EXPECT(inc_energy_ > gamma_cutoff_);
}

//---------------------------------------------------------------------------//
Expand All @@ -129,13 +130,6 @@ CombinedBremInteractor::CombinedBremInteractor(
template<class Engine>
CELER_FUNCTION Interaction CombinedBremInteractor::operator()(Engine& rng)
{
// TODO: reject this interaction before executing the kernel by using
// correct material-dependent lower bounds for the interaction
if (gamma_cutoff_ > inc_energy_)
{
return Interaction::from_unchanged(inc_energy_, inc_direction_);
}

// Allocate space for the brems photon
Secondary* secondaries = allocate_(1);
if (secondaries == nullptr)
Expand All @@ -146,7 +140,7 @@ CELER_FUNCTION Interaction CombinedBremInteractor::operator()(Engine& rng)

// Sample the bremsstrahlung photon energy
Energy gamma_energy;
if (inc_energy_ > detail::seltzer_berger_limit())
if (inc_energy_ >= detail::seltzer_berger_limit())
{
detail::RBEnergySampler sample_energy{
shared_.rb_data, inc_energy_, cutoffs_, material_, elcomp_id_};
Expand Down
10 changes: 2 additions & 8 deletions src/celeritas/em/interactor/MollerBhabhaInteractor.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -104,6 +104,8 @@ CELER_FUNCTION MollerBhabhaInteractor::MollerBhabhaInteractor(
{
CELER_EXPECT(particle.particle_id() == shared_.ids.electron
|| particle.particle_id() == shared_.ids.positron);
CELER_EXPECT(inc_energy_
> (inc_particle_is_electron_ ? 2 : 1) * electron_cutoff_);
}

//---------------------------------------------------------------------------//
Expand All @@ -116,14 +118,6 @@ CELER_FUNCTION MollerBhabhaInteractor::MollerBhabhaInteractor(
template<class Engine>
CELER_FUNCTION Interaction MollerBhabhaInteractor::operator()(Engine& rng)
{
if (inc_energy_ <= (inc_particle_is_electron_ ? 2 : 1) * electron_cutoff_)
{
// The secondary should not be emitted. This interaction cannot
// happen and the incident particle must undergo an energy loss
// process.
return Interaction::from_unchanged(Energy{inc_energy_}, inc_direction_);
}

// Allocate memory for the produced electron
Secondary* electron_secondary = allocate_(1);

Expand Down
2 changes: 1 addition & 1 deletion src/celeritas/em/interactor/RelativisticBremInteractor.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -113,7 +113,7 @@ RelativisticBremInteractor::RelativisticBremInteractor(
|| particle.particle_id() == shared_.ids.positron);

// Valid energy region of the relativistic e-/e+ Bremsstrahlung model
CELER_EXPECT(inc_energy_ > detail::seltzer_berger_limit());
CELER_EXPECT(inc_energy_ >= detail::seltzer_berger_limit());
}

//---------------------------------------------------------------------------//
Expand Down
12 changes: 2 additions & 10 deletions src/celeritas/em/interactor/SeltzerBergerInteractor.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -128,6 +128,8 @@ CELER_FUNCTION SeltzerBergerInteractor::SeltzerBergerInteractor(
CELER_EXPECT(particle.particle_id() == shared_.ids.electron
|| particle.particle_id() == shared_.ids.positron);
CELER_EXPECT(gamma_cutoff_ > zero_quantity());
CELER_EXPECT(inc_energy_ > gamma_cutoff_
&& inc_energy_ < detail::seltzer_berger_limit());
}

//---------------------------------------------------------------------------//
Expand All @@ -139,16 +141,6 @@ CELER_FUNCTION SeltzerBergerInteractor::SeltzerBergerInteractor(
template<class Engine>
CELER_FUNCTION Interaction SeltzerBergerInteractor::operator()(Engine& rng)
{
// Check if secondary can be produced. If not, this interaction cannot
// happen and the incident particle must undergo an energy loss process
// instead.
// TODO: reject this interaction before executing the kernel by using
// correct material-dependent lower bounds for the interaction
if (gamma_cutoff_ > inc_energy_)
{
return Interaction::from_unchanged(inc_energy_, inc_direction_);
}

// Allocate space for the brems photon
Secondary* secondaries = allocate_(1);
if (secondaries == nullptr)
Expand Down
6 changes: 4 additions & 2 deletions src/celeritas/em/model/BetheHeitlerModel.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -57,10 +57,12 @@ BetheHeitlerModel::BetheHeitlerModel(ActionId id,
*/
auto BetheHeitlerModel::applicability() const -> SetApplicability
{
using Energy = units::MevEnergy;

Applicability photon_applic;
photon_applic.particle = data_.ids.gamma;
photon_applic.lower = zero_quantity();
photon_applic.upper = units::MevEnergy{1e8};
photon_applic.lower = Energy{2 * this->host_ref().electron_mass.value()};
photon_applic.upper = Energy{1e8};

return {photon_applic};
}
Expand Down
2 changes: 1 addition & 1 deletion src/celeritas/em/model/EPlusGGModel.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -46,7 +46,7 @@ auto EPlusGGModel::applicability() const -> SetApplicability
{
Applicability applic;
applic.particle = data_.ids.positron;
applic.lower = neg_max_quantity(); // Valid at rest
applic.lower = zero_quantity(); // Valid at rest
applic.upper = units::MevEnergy{1e8}; // 100 TeV

return {applic};
Expand Down
18 changes: 3 additions & 15 deletions src/celeritas/phys/Applicability.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,10 +23,9 @@ namespace celeritas
* Range where a model and/or process is valid.
*
* This class is used during setup for specifying the ranges of applicability
* for a physics model or process. The interval is *open* on the lower energy
* range and *closed* on the upper energy. So a threshold reaction should have
* the lower energy set to the threshold. Models valid to zero energy but have
* special "at rest" models should set upper to zero.
* for a physics model or process. The interval is *closed* on the lower energy
* range and *open* on the upper energy. So a threshold reaction should have
* the lower energy set to the threshold.
*
* An unset value for "material" means it applies to all materials; however,
* the particle ID should always be set.
Expand All @@ -41,17 +40,6 @@ struct Applicability
Energy lower = zero_quantity();
Energy upper = max_quantity();

//! Range for a particle at rest
static inline Applicability at_rest(ParticleId id)
{
CELER_EXPECT(id);
Applicability result;
result.particle = id;
result.lower = neg_max_quantity();
result.upper = zero_quantity();
return result;
}

//! Whether applicability is in a valid state
inline explicit operator bool() const
{
Expand Down
9 changes: 9 additions & 0 deletions test/celeritas/em/BetheHeitler.test.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -274,6 +274,15 @@ TEST_F(BetheHeitlerInteractorTest, distributions)
= {1209, 1073, 911, 912, 844, 881, 903, 992, 1066, 1209};
EXPECT_VEC_EQ(expected_eps_dist, eps_dist);
}

// Interaction threshold energy
{
std::vector<int> eps_dist
= bin_epsilon(2 * data_.electron_mass.value());
static int const expected_eps_dist[]
= {0, 0, 0, 0, 0, 10000, 0, 0, 0, 0};
EXPECT_VEC_EQ(expected_eps_dist, eps_dist);
}
}

//---------------------------------------------------------------------------//
Expand Down
53 changes: 26 additions & 27 deletions test/celeritas/em/CombinedBrem.test.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -271,7 +271,6 @@ TEST_F(CombinedBremTest, stress_test_combined)
for (real_type inc_e : test_energy)
{
SCOPED_TRACE("Incident energy: " + std::to_string(inc_e));
// this->set_inc_particle(particle, MevEnergy{inc_e});

RandomEngine& rng_engine = this->rng();
RandomEngine::size_type num_particles_sampled = 0;
Expand Down Expand Up @@ -320,37 +319,37 @@ TEST_F(CombinedBremTest, stress_test_combined)
12.9641,
12.5832,
12.4988,
12.3433,
12.4378,
13.2556,
15.3633,
14.2262,
13.262,
12.9294,
12.5754,
12.508,
12.3334,
12.4193,
13.293,
15.3784};
12.31,
12.4381,
13.2552,
15.3604,
14.2257,
13.2616,
12.9286,
12.5763,
12.5076,
12.3059,
12.4207,
13.2906,
15.3809};
static double const expected_avg_energy_samples[] = {0.20338654094171,
0.53173619503507,
0.99638562846318,
4.4359411867158,
8.7590072534526,
85.185116736899,
905.94487251514,
10719.081816783,
149600.77957549,
0.18914626656986,
0.52230134540886,
0.98770452529095,
4.4238993615396,
8.4950149725315,
85.418339892001,
917.61799096706,
10758.713294023,
146932.68621334};
85.769352932541,
905.74010590236,
10724.127172904,
149587.47778726,
0.18927693789863,
0.52259561993542,
0.98783539434815,
4.4286338859014,
8.495313663667,
85.853825010643,
917.26619467326,
10760.512214274,
146925.51055711};

EXPECT_VEC_SOFT_EQ(expected_avg_engine_samples, avg_engine_samples);
EXPECT_VEC_SOFT_EQ(expected_avg_energy_samples, avg_energy_samples);
Expand Down

0 comments on commit 88237d1

Please sign in to comment.