From 843c0309f3e349df2853c347837c7233cbb3441b Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Mon, 25 Oct 2021 22:20:48 +0200
Subject: [PATCH] more unit testing
---
.../framework/process/ProcessSequence.inl | 10 +-
.../process/SwitchProcessSequence.inl | 86 +++++------
examples/mars.cpp | 2 +-
tests/framework/testProcessSequence.cpp | 135 +++++++++++++-----
4 files changed, 151 insertions(+), 82 deletions(-)
diff --git a/corsika/detail/framework/process/ProcessSequence.inl b/corsika/detail/framework/process/ProcessSequence.inl
index 024c82809..5711ff249 100644
--- a/corsika/detail/framework/process/ProcessSequence.inl
+++ b/corsika/detail/framework/process/ProcessSequence.inl
@@ -312,10 +312,10 @@ namespace corsika {
int IndexProcess2>
template <typename TParticle>
inline CrossSectionType
- ProcessSequence<TProcess1, TProcess2, IndexStart, IndexProcess1,
- IndexProcess2>::getCrossSection(TParticle const& projectile,
- Code const targetId,
- FourMomentum const& targetP4) const {
+ ProcessSequence<TProcess1, TProcess2, IndexStart, IndexProcess1, IndexProcess2>::
+ getCrossSection([[maybe_unused]] TParticle const& projectile,
+ [[maybe_unused]] Code const targetId,
+ [[maybe_unused]] FourMomentum const& targetP4) const {
CrossSectionType tot = CrossSectionType::zero();
@@ -338,8 +338,8 @@ namespace corsika {
} else if constexpr (process2_type::is_process_sequence) {
tot += B_.getCrossSection(projectile, targetId, targetP4);
}
- return tot;
}
+ return tot;
}
template <typename TProcess1, typename TProcess2, int IndexStart, int IndexProcess1,
diff --git a/corsika/detail/framework/process/SwitchProcessSequence.inl b/corsika/detail/framework/process/SwitchProcessSequence.inl
index 02e197ec0..fdf1e685e 100644
--- a/corsika/detail/framework/process/SwitchProcessSequence.inl
+++ b/corsika/detail/framework/process/SwitchProcessSequence.inl
@@ -256,17 +256,18 @@ namespace corsika {
auto const& projectile = view.parent();
Code const projectileId = projectile.getPID();
- unsigned int const projectileA = projectile.getNuclearA();
// get cross section vector for all material components
+ // for selected process A
static_assert(
has_method_getCrossSection_v<TSequence, // process object
CrossSectionType, // return type
- Code, // parameters
- Code, FourMomentum const&, FourMomentum const&>,
+ Code, Code, // parameters
+ FourMomentum const&, FourMomentum const&>,
"TSequence has no method with correct signature \"CrossSectionType "
- "getCrossSection(Code, Code, FourMomentum const&, FourMomntum const&)\" "
- "required by InteractionProcess<TSequence>. ");
+ "getCrossSection(Code, Code, FourMomentum const&, FourMomentum "
+ "const&)\" required by "
+ "InteractionProcess<TSequence>. ");
std::vector<CrossSectionType> const weightedCrossSections =
composition.getWeighted([=](Code const targetId) -> CrossSectionType {
@@ -289,21 +290,21 @@ namespace corsika {
{0_GeV, 0_GeV, 0_GeV}));
// interface checking on TProcess1
- static_assert(has_method_doInteract_v<TSequence, // process object
- void, // return type
- TSecondaryView, // template argument
- TSecondaryView&, // method parameters
- Code, Code, FourMomentum const&,
- FourMomentum const&>,
- "USequence has no method with correct signature \"void "
- "doInteraction<TSecondaryView>(TSecondaryView&, Code, "
- "Code, FourMomentum const&, FourMomntum const&)\" required for "
- "InteractionProcess<USequence>. ");
+ static_assert(
+ has_method_doInteract_v<TSequence, // process object
+ void, // return type
+ TSecondaryView, // template argument
+ TSecondaryView&, // method parameters
+ Code, Code, FourMomentum const&,
+ FourMomentum const&>,
+ "TSequence has no method with correct signature \"void "
+ "doInteraction<TSecondaryView>(TSecondaryView&, "
+ "Code, Code, FourMomentum const&, FourMomentum const&)\" required for "
+ "InteractionProcess<TSequence>. ");
A_.template doInteraction(view, projectileId, targetId, projectileP4, targetP4);
return ProcessReturn::Interacted;
-
} // end collision branch A
}
@@ -317,17 +318,16 @@ namespace corsika {
auto const& projectile = view.parent();
Code const projectileId = projectile.getPID();
- unsigned int const projectileA = projectile.getNuclearA();
- // get cross section vector for all material components
- static_assert(
- has_method_getCrossSection_v<USequence, // process object
- CrossSectionType, // return type
- Code, // parameters
- Code, FourMomentum const&, FourMomentum const&>,
- "USequence has no method with correct signature \"CrossSectionType "
- "getCrossSection(Code, Code, FourMomentum const&, FourMomentum const&)\" "
- "required by InteractionProcess<USequence>. ");
+ // get cross section vector for all material components, for selected process B
+ static_assert(has_method_getCrossSection_v<USequence, // process object
+ CrossSectionType, // return type
+ Code, Code, FourMomentum const&,
+ FourMomentum const&>, // parameters
+ "USequence has no method with correct signature \"CrossSectionType "
+ "getCrossSection(Code, Code, FourMomentum const&, FourMomentum "
+ "const&)\" required by "
+ "InteractionProcess<USequence>. ");
std::vector<CrossSectionType> const weightedCrossSections =
composition.getWeighted([=](Code const targetId) -> CrossSectionType {
@@ -338,10 +338,11 @@ namespace corsika {
return B_.getCrossSection(projectileId, targetId, projectileP4, targetP4);
});
- cx_sum += std::accumulate(weightedCrossSections.cbegin(),
- weightedCrossSections.cend(), CrossSectionType::zero());
+ cx_sum += std::accumulate(weightedCrossSections.begin(),
+ weightedCrossSections.end(), CrossSectionType::zero());
- if (cx_select < cx_sum) {
+ // check if we should execute THIS process and then EXIT
+ if (cx_select <= cx_sum) {
// now also sample targetId from weighted cross sections
Code const targetId = composition.sampleTarget(weightedCrossSections, rng);
@@ -350,19 +351,20 @@ namespace corsika {
MomentumVector(projectile.getMomentum().getCoordinateSystem(),
{0_GeV, 0_GeV, 0_GeV}));
- // interface checking on TProcess1
- static_assert(has_method_doInteract_v<USequence, // process object
- void, // return type
- TSecondaryView, // template argument
- TSecondaryView&, // method parameters
- Code, Code, FourMomentum const&,
- FourMomentum const&>,
- "USequence has no method with correct signature \"void "
- "doInteraction<TSecondaryView>(TSecondaryView&, Code, "
- "Code, FourMomentum const&, FourMomentum const&)\" required for "
- "InteractionProcess<USequence>. ");
-
- B_.template doInteraction(view, projectileId, targetId, projectileP4, targetP4);
+ // interface checking on TProcess2
+ static_assert(
+ has_method_doInteract_v<USequence, // process object
+ void, // return type
+ TSecondaryView, // template argument
+ TSecondaryView&, // method parameters
+ Code, Code, FourMomentum const&,
+ FourMomentum const&>,
+ "USequence has no method with correct signature \"void "
+ "doInteraction<TSecondaryView>(TSecondaryView&, "
+ "Code, Code, FourMomentum const&, FourMomentum const&)\" required for "
+ "InteractionProcess<USequence>. ");
+
+ B_.doInteraction(view, projectileId, targetId, projectileP4, targetP4);
return ProcessReturn::Interacted;
} // end collision in branch B
diff --git a/examples/mars.cpp b/examples/mars.cpp
index e5deda391..60b4f87d4 100644
--- a/examples/mars.cpp
+++ b/examples/mars.cpp
@@ -360,7 +360,7 @@ int main(int argc, char** argv) {
HEPEnergyType cutE_;
EnergySwitch(HEPEnergyType cutE)
: cutE_(cutE) {}
- bool operator()(const Particle& p) { return (p.getKineticEnergy() < cutE_); }
+ bool operator()(Particle const& p) const { return (p.getKineticEnergy() < cutE_); }
};
auto hadronSequence = make_select(EnergySwitch(63.1_GeV), urqmdCounted, heModelCounted);
auto decaySequence = make_sequence(decayPythia, decaySibyll);
diff --git a/tests/framework/testProcessSequence.cpp b/tests/framework/testProcessSequence.cpp
index ac4c6392d..7cf2c55c4 100644
--- a/tests/framework/testProcessSequence.cpp
+++ b/tests/framework/testProcessSequence.cpp
@@ -63,7 +63,6 @@ struct DummyData {
return MomentumVector{get_root_CoordinateSystem(), 0_eV, 0_eV, 0_eV};
}
HEPEnergyType getEnergy() const { return 10_GeV; }
- unsigned int getNuclearA() const { return 1; }
};
// there is no real trajectory/track
@@ -454,27 +453,6 @@ TEST_CASE("ProcessSequence General", "ProcessSequence") {
sequence2_rv.getProcess2().getProcess2())>); // Process3
}
- SECTION("interaction length") {
- globalCount = 0;
- ContinuousProcess1 cp1(0, 1_m);
- Process2 m2(1);
- Process3 m3(2);
-
- DummyData particle;
-
- auto sequence2 = make_sequence(cp1, m2, m3);
- /*GrammageType const tot = sequence2.getInteractionLength(particle);
- InverseGrammageType const tot_inv = sequence2.getInverseInteractionLength(particle);
- CORSIKA_LOG_DEBUG(
- "lambda_tot={}"
- "; lambda_tot_inv={}",
- tot, tot_inv);
-
- CHECK(tot / 1_g * square(1_cm) == 12);
- CHECK(tot_inv * 1_g / square(1_cm) == 1. / 12);*/
- globalCount = 0;
- }
-
SECTION("lifetime") {
globalCount = 0;
ContinuousProcess1 cp1(0, 1_m);
@@ -640,12 +618,20 @@ TEST_CASE("SwitchProcessSequence", "ProcessSequence") {
auto sec1 = Secondaries1();
auto sec2 = Secondaries2();
- auto sequence1 = make_sequence(Process1(0), cp2, Decay1(0), sec1, Boundary1(1.0));
- auto sequence2 = make_sequence(cp3, Process2(0), Boundary1(-1.0), Decay2(0), sec2);
+ auto sequence1 =
+ make_sequence(Process1(0), cp2, Decay1(0), sec1, Boundary1(1.0)); // 10 mb
+ auto sequence2 =
+ make_sequence(cp3, Process2(0), Boundary1(-1.0), Decay2(0), sec2); // 20 mb
- auto sequence3 = make_sequence(cp1, Process3(0),
+ auto sequence3 = make_sequence(cp1, Process3(0), // 30 mb
SwitchProcessSequence(select1, sequence1, sequence2));
+ // it is even more typical to have just one sub-process inside the branches of
+ // SwitchProcessSequence
+ auto sequence3_short =
+ make_sequence(cp1, Process3(0), // 30 mb
+ SwitchProcessSequence(select1, Process1(0), Process2(0)));
+
auto sequence4 =
make_sequence(cp1, Boundary1(2.0), Process3(0),
SwitchProcessSequence(select1, sequence1, Boundary1(-1.0)));
@@ -778,15 +764,96 @@ TEST_CASE("SwitchProcessSequence", "ProcessSequence") {
CHECK(checkCont == 0);
CHECK(checkSec == 0);
- // check that large "select" value will correctly ignore the call
- cx_select = 1e5_mb;
- time_select = 1e5 / second;
- checkDecay = 0;
- checkInteract = 0;
- sequence3.selectInteraction(view, projectileP4, noComposition, rng, cx_select);
- sequence3.selectDecay(view, time_select);
- CHECK(checkInteract == 0);
- CHECK(checkDecay == 0);
+ // now check sequence3, which contains a SwitchProcessSequence that contains two
+ // longer sequences in each branch.
+ {
+ // check that large "select" value will correctly ignore the call
+ cx_select = 1e5_mb;
+ time_select = 1e5 / second;
+ checkDecay = 0;
+ checkInteract = 0;
+ sequence3.selectInteraction(view, projectileP4, noComposition, rng, cx_select);
+ sequence3.selectDecay(view, time_select);
+ CHECK(checkInteract == 0);
+ CHECK(checkDecay == 0);
+
+ // for a small cx_select selection must be sucessful
+ cx_select = 28_mb; // -> Process3
+ checkInteract = 0;
+ particle.data_[0] = -100; // data negative --> sequence2
+ CHECK(sequence3.getCrossSection(particle, Code::Oxygen,
+ {Oxygen::mass, {rootCS, {0_eV, 0_eV, 0_eV}}}) /
+ 1_mb ==
+ Approx(50.));
+ sequence3.selectInteraction(view, projectileP4, noComposition, rng, cx_select);
+ CHECK(checkInteract == 4); // 2^3
+
+ particle.data_[0] = 100; // data positive --> sequence1
+ checkInteract = 0;
+ CHECK(sequence3.getCrossSection(particle, Code::Oxygen,
+ {Oxygen::mass, {rootCS, {0_eV, 0_eV, 0_eV}}}) /
+ 1_mb ==
+ Approx(40.));
+ sequence3.selectInteraction(view, projectileP4, noComposition, rng, cx_select);
+ CHECK(checkInteract == 4); // 2^3
+
+ cx_select = 32_mb; // -> Process2 or Process1
+ checkInteract = 0;
+ particle.data_[0] = -100; // data negative --> Process2
+ sequence3.selectInteraction(view, projectileP4, noComposition, rng, cx_select);
+ CHECK(checkInteract == 2); // 2^2
+
+ particle.data_[0] = 100; // data positive --> Process1
+ checkInteract = 0;
+ sequence3.selectInteraction(view, projectileP4, noComposition, rng, cx_select);
+ CHECK(checkInteract == 1); // 2^1
+ }
+
+ // now check sequence3, which contains a SwitchProcessSequence that contains just two
+ // bare InteractionProcess-es in each branch.
+ {
+ // check that large "select" value will correctly ignore the call
+ cx_select = 1e5_mb;
+ checkInteract = 0;
+ sequence3_short.selectInteraction(view, projectileP4, noComposition, rng,
+ cx_select);
+ CHECK(checkInteract == 0);
+
+ // for a small cx_select selection must be sucessful
+ cx_select = 28_mb; // -> Process3
+ checkInteract = 0;
+ particle.data_[0] = -100; // data negative --> sequence2
+ CHECK(sequence3_short.getCrossSection(
+ particle, Code::Oxygen, {Oxygen::mass, {rootCS, {0_eV, 0_eV, 0_eV}}}) /
+ 1_mb ==
+ Approx(50.));
+ sequence3_short.selectInteraction(view, projectileP4, noComposition, rng,
+ cx_select);
+ CHECK(checkInteract == 4); // 2^3
+
+ particle.data_[0] = 100; // data positive --> sequence1
+ checkInteract = 0;
+ CHECK(sequence3_short.getCrossSection(
+ particle, Code::Oxygen, {Oxygen::mass, {rootCS, {0_eV, 0_eV, 0_eV}}}) /
+ 1_mb ==
+ Approx(40.));
+ sequence3_short.selectInteraction(view, projectileP4, noComposition, rng,
+ cx_select);
+ CHECK(checkInteract == 4); // 2^3
+
+ cx_select = 32_mb; // -> Process2 or Process1
+ checkInteract = 0;
+ particle.data_[0] = -100; // data negative --> Process2
+ sequence3_short.selectInteraction(view, projectileP4, noComposition, rng,
+ cx_select);
+ CHECK(checkInteract == 2); // 2^2
+
+ particle.data_[0] = 100; // data positive --> Process1
+ checkInteract = 0;
+ sequence3_short.selectInteraction(view, projectileP4, noComposition, rng,
+ cx_select);
+ CHECK(checkInteract == 1); // 2^1
+ }
}
SECTION("Check SecondariesProcesses in SwitchProcessSequence") {
--
GitLab