From 9cc0c321e6e54d3e06c1f1bd45388c092fb3aada Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Wed, 7 Oct 2020 15:09:24 +0200
Subject: [PATCH] cleanup of process tests
---
Processes/Pythia/testPythia8.cc | 55 +++--------------
Processes/QGSJetII/testQGSJetII.cc | 40 ++-----------
Processes/Sibyll/testSibyll.cc | 96 ++++--------------------------
Processes/UrQMD/testUrQMD.cc | 14 ++---
Setup/SetupStack.h | 23 +++----
5 files changed, 47 insertions(+), 181 deletions(-)
diff --git a/Processes/Pythia/testPythia8.cc b/Processes/Pythia/testPythia8.cc
index 4fd221bd6..a0af2a46a 100644
--- a/Processes/Pythia/testPythia8.cc
+++ b/Processes/Pythia/testPythia8.cc
@@ -101,46 +101,18 @@ auto sumMomentum(TStackView const& view, geometry::CoordinateSystem const& vCS)
TEST_CASE("pythia process") {
- // setup environment, geometry
- setup::Environment env;
- auto& universe = *(env.GetUniverse());
- using EnvironmentModel = environment::UniformMediumType<environment::UniformMagneticField<environment::HomogeneousMedium<setup::IEnvironment>>>;
-
- auto theMedium =
- setup::Environment::CreateNode<geometry::Sphere>(
- geometry::Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
- 1_km * std::numeric_limits<double>::infinity());
-
- theMedium->SetModelProperties<EnvironmentModel>(
- environment::EMediumType::eAir,
- geometry::Vector(env.GetCoordinateSystem(), 0_T, 0_T, 0_T),
- 1_kg / (1_m * 1_m * 1_m),
- environment::NuclearComposition(
- std::vector<particles::Code>{particles::Code::Oxygen}, std::vector<float>{1.}));
-
- auto const* nodePtr = theMedium.get();
- universe.AddChild(std::move(theMedium));
-
- const geometry::CoordinateSystem& cs = env.GetCoordinateSystem();
-
+ auto [env, csPtr, nodePtr] = testing::setupEnvironment(particles::Code::Oxygen);
+ auto const& cs = *csPtr;
+ [[maybe_unused]] auto const& env_dummy = env;
+ [[maybe_unused]] auto const& node_dummy = nodePtr;
SECTION("pythia decay") {
feenableexcept(FE_INVALID);
- setup::Stack stack;
- const HEPEnergyType E0 = 10_GeV;
- HEPMomentumType P0 =
- sqrt(E0 * E0 - particles::PiPlus::GetMass() * particles::PiPlus::GetMass());
- auto plab = corsika::stack::MomentumVector(cs, {0_GeV, 0_GeV, P0});
- geometry::Point pos(cs, 0_m, 0_m, 0_m);
- auto particle = stack.AddParticle(
- std::tuple<particles::Code, units::si::HEPEnergyType,
- corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
- particles::Code::PiPlus, E0, plab, pos, 0_ns});
+ auto [stackPtr, secViewPtr] = testing::setupStack(code::PiPlus, 0, 0, 10_GeV, nodePtr, *csPtr);
+ auto projectile = secViewPtr->GetProjectile();
random::RNGManager::GetInstance().RegisterRandomStream("pythia");
- setup::StackView view(particle);
-
process::pythia::Decay model;
[[maybe_unused]] const TimeType time = model.GetLifetime(particle);
@@ -177,18 +149,9 @@ TEST_CASE("pythia process") {
SECTION("pythia interaction") {
- setup::Stack stack;
- const HEPEnergyType E0 = 100_GeV;
- HEPMomentumType P0 =
- sqrt(E0 * E0 - particles::PiPlus::GetMass() * particles::PiPlus::GetMass());
- auto plab = corsika::stack::MomentumVector(cs, {0_GeV, 0_GeV, -P0});
- geometry::Point pos(cs, 0_m, 0_m, 0_m);
- auto particle = stack.AddParticle(
- std::tuple<particles::Code, units::si::HEPEnergyType,
- corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
- particles::Code::PiPlus, E0, plab, pos, 0_ns});
- particle.SetNode(nodePtr);
- setup::StackView view(particle);
+ feenableexcept(FE_INVALID);
+ auto [stackPtr, secViewPtr] = testing::setupStack(code::PiPlus, 0, 0, 100_GeV, nodePtr, *csPtr);
+ auto projectile = secViewPtr->GetProjectile();
process::pythia::Interaction model;
diff --git a/Processes/QGSJetII/testQGSJetII.cc b/Processes/QGSJetII/testQGSJetII.cc
index d1e9c7783..7f9af0675 100644
--- a/Processes/QGSJetII/testQGSJetII.cc
+++ b/Processes/QGSJetII/testQGSJetII.cc
@@ -125,45 +125,17 @@ using namespace corsika::units;
TEST_CASE("QgsjetIIInterface", "[processes]") {
- // setup environment, geometry
- setup::Environment env;
- auto& universe = *(env.GetUniverse());
- using EnvironmentModel = environment::UniformMediumType<environment::UniformMagneticField<environment::HomogeneousMedium<setup::IEnvironment>>>;
-
- auto theMedium =
- setup::Environment::CreateNode<geometry::Sphere>(
- geometry::Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
- 1_km * std::numeric_limits<double>::infinity());
-
- theMedium->SetModelProperties<EnvironmentModel>(
- environment::EMediumType::eAir,
- geometry::Vector(env.GetCoordinateSystem(), 0_T, 0_T, 0_T),
- 1_kg / (1_m * 1_m * 1_m),
- environment::NuclearComposition(
- std::vector<particles::Code>{particles::Code::Oxygen}, std::vector<float>{1.}));
-
- auto const* nodePtr = theMedium.get();
- universe.AddChild(std::move(theMedium));
-
- const geometry::CoordinateSystem& cs = env.GetCoordinateSystem();
+ auto [env, csPtr, nodePtr] = testing::setupEnvironment(particles::Code::Oxygen);
+ auto const& cs = *csPtr;
+ [[maybe_unused]] auto const& env_dummy = env;
+ [[maybe_unused]] auto const& node_dummy = nodePtr;
corsika::random::RNGManager::GetInstance().RegisterRandomStream("qgsjet");
SECTION("InteractionInterface") {
- setup::Stack stack;
- const HEPEnergyType E0 = 100_GeV;
- HEPMomentumType P0 =
- sqrt(E0 * E0 - particles::Proton::GetMass() * particles::Proton::GetMass());
- auto plab = corsika::stack::MomentumVector(cs, {0_GeV, 0_GeV, -P0});
- geometry::Point pos(cs, 0_m, 0_m, 0_m);
- auto particle = stack.AddParticle(
- std::tuple<particles::Code, units::si::HEPEnergyType,
- corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
- particles::Code::Proton, E0, plab, pos, 0_ns});
-
- particle.SetNode(nodePtr);
- setup::StackView view(particle);
+ auto [stackPtr, secViewPtr] =
+ testing::setupStack(particles::Code::Proton, 0,0, 110_GeV, nodePtr, *csPtr);
auto projectile = view.GetProjectile();
auto const projectileMomentum = projectile.GetMomentum();
diff --git a/Processes/Sibyll/testSibyll.cc b/Processes/Sibyll/testSibyll.cc
index 755237efa..ccdb4b651 100644
--- a/Processes/Sibyll/testSibyll.cc
+++ b/Processes/Sibyll/testSibyll.cc
@@ -74,15 +74,15 @@ TEST_CASE("Sibyll", "[processes]") {
#include <corsika/units/PhysicalUnits.h>
#include <corsika/particles/ParticleProperties.h>
-#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupEnvironment.h>
+#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/environment/Environment.h>
#include <corsika/environment/HomogeneousMedium.h>
#include <corsika/environment/NuclearComposition.h>
-#include <corsika/environment/UniformMediumType.h>
#include <corsika/environment/UniformMagneticField.h>
+#include <corsika/environment/UniformMediumType.h>
using namespace corsika::units::si;
using namespace corsika::units;
@@ -96,42 +96,17 @@ auto sumMomentum(TStackView const& view, geometry::CoordinateSystem const& vCS)
TEST_CASE("SibyllInterface", "[processes]") {
- // setup environment, geometry
- setup::Environment env;
- auto& universe = *(env.GetUniverse());
- using EnvironmentModel = environment::UniformMediumType<environment::UniformMagneticField<environment::HomogeneousMedium<setup::IEnvironment>>>;
-
- auto theMedium =
- setup::Environment::CreateNode<geometry::Sphere>(
- geometry::Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
- 1_km * std::numeric_limits<double>::infinity());
-
- theMedium->SetModelProperties<EnvironmentModel>(
- environment::EMediumType::eAir,
- geometry::Vector(env.GetCoordinateSystem(), 0_T, 0_T, 0_T),
- 1_kg / (1_m * 1_m * 1_m),
- environment::NuclearComposition(
- std::vector<particles::Code>{particles::Code::Oxygen}, std::vector<float>{1.}));
-
- auto const* nodePtr = theMedium.get();
- universe.AddChild(std::move(theMedium));
-
- const geometry::CoordinateSystem& cs = env.GetCoordinateSystem();
+ auto [env, csPtr, nodePtr] = testing::setupEnvironment(particles::Code::Oxygen);
+ auto const& cs = *csPtr;
+ [[maybe_unused]] auto const& env_dummy = env;
+ [[maybe_unused]] auto const& node_dummy = nodePtr;
random::RNGManager::GetInstance().RegisterRandomStream("sibyll");
SECTION("InteractionInterface - low energy") {
- setup::Stack stack;
- const HEPEnergyType E0 = 60_GeV;
- HEPMomentumType P0 =
- sqrt(E0 * E0 - particles::Proton::GetMass() * particles::Proton::GetMass());
- auto plab = corsika::stack::MomentumVector(cs, {P0, 0_eV, 0_eV});
- geometry::Point pos(cs, 0_m, 0_m, 0_m);
- auto particle =
- stack.AddParticle(std::make_tuple(particles::Code::Proton, E0, plab, pos, 0_ns));
- particle.SetNode(nodePtr);
- corsika::setup::StackView view(particle);
+ auto [stack, view] = testing::setupStack(particles::Code::Proton, 60_GeV, nodePtr, cs);
+ auto projectile = view.GetProjectile();
Interaction model;
@@ -203,45 +178,10 @@ TEST_CASE("SibyllInterface", "[processes]") {
[[maybe_unused]] const GrammageType length = model.GetInteractionLength(particle);
}
- SECTION("InteractionInterface - high energy") {
-
- setup::Stack stack;
- const HEPEnergyType E0 = 60_EeV;
- HEPMomentumType P0 =
- sqrt(E0 * E0 - particles::Proton::GetMass() * particles::Proton::GetMass());
- auto plab = corsika::stack::MomentumVector(cs, {P0, 0_eV, 0_eV});
- geometry::Point pos(cs, 0_m, 0_m, 0_m);
- auto particle =
- stack.AddParticle(std::make_tuple(particles::Code::Proton, E0, plab, pos, 0_ns));
- particle.SetNode(nodePtr);
- corsika::setup::StackView view(particle);
-
- Interaction model;
-
- [[maybe_unused]] const process::EProcessReturn ret = model.DoInteraction(view);
- auto const pSum = sumMomentum(view, cs);
- CHECK(pSum.GetComponents(cs).GetX() / P0 == Approx(1).margin(0.001));
- CHECK(pSum.GetComponents(cs).GetY() / 1_GeV == Approx(0).margin(1e-4));
- CHECK(pSum.GetComponents(cs).GetZ() / 1_GeV == Approx(0).margin(1e-4));
-
- CHECK((pSum - plab).norm() / 1_GeV == Approx(0).margin(plab.norm() * 0.001 / 1_GeV));
- CHECK(pSum.norm() / P0 == Approx(1).margin(0.05));
- [[maybe_unused]] const GrammageType length = model.GetInteractionLength(particle);
- }
-
SECTION("NuclearInteractionInterface") {
- setup::Stack stack;
- const HEPEnergyType E0 = 400_GeV;
- HEPMomentumType P0 =
- sqrt(E0 * E0 - particles::Proton::GetMass() * particles::Proton::GetMass());
- auto plab = corsika::stack::MomentumVector(cs, {0_GeV, 0_GeV, -P0});
- geometry::Point pos(cs, 0_m, 0_m, 0_m);
-
- auto particle = stack.AddParticle(
- std::make_tuple(particles::Code::Nucleus, E0, plab, pos, 0_ns, 4, 2));
- particle.SetNode(nodePtr);
- corsika::setup::StackView view(particle);
+ auto [stack, view] = testing::setupStack(particles::Code::Nucleus, 4, 2, 100_GeV, nodePtr, cs);
+ auto projectile = view.GetProjectile();
Interaction hmodel;
NuclearInteraction model(hmodel, env);
@@ -252,24 +192,14 @@ TEST_CASE("SibyllInterface", "[processes]") {
SECTION("DecayInterface") {
- setup::Stack stack;
- const HEPEnergyType E0 = 10_GeV;
- HEPMomentumType P0 =
- sqrt(E0 * E0 - particles::Proton::GetMass() * particles::Proton::GetMass());
- auto plab = corsika::stack::MomentumVector(cs, {0_GeV, 0_GeV, -P0});
- geometry::Point pos(cs, 0_m, 0_m, 0_m);
- auto particle =
- stack.AddParticle(std::make_tuple(particles::Code::Lambda0, E0, plab, pos, 0_ns));
- corsika::setup::StackView view(particle);
+ auto [stack, view] = testing::setupStack(particles::Code::Lambda0, 0,0, 10_GeV, nodePtr, cs);
+ auto projectile = view.GetProjectile();
Decay model;
-
model.PrintDecayConfig();
-
[[maybe_unused]] const TimeType time = model.GetLifetime(particle);
- /*[[maybe_unused]] const process::EProcessReturn ret =*/model.DoDecay(view);
-
+ /*[[maybe_unused]] const process::EProcessReturn ret =*/model.DoDecay(projectile);
// run checks
// lambda decays into proton and pi- or neutron and pi+
CHECK(stack.getEntries() == 3);
diff --git a/Processes/UrQMD/testUrQMD.cc b/Processes/UrQMD/testUrQMD.cc
index 537b656c6..7bd929e3e 100644
--- a/Processes/UrQMD/testUrQMD.cc
+++ b/Processes/UrQMD/testUrQMD.cc
@@ -67,7 +67,7 @@ TEST_CASE("UrQMD") {
UrQMD urqmd;
SECTION("interaction length") {
- auto [env, csPtr, nodePtr] = setupEnvironment(particles::Code::Nitrogen);
+ auto [env, csPtr, nodePtr] = testing::setupEnvironment(particles::Code::Nitrogen);
auto const& cs = *csPtr;
[[maybe_unused]] auto const& env_dummy = env;
[[maybe_unused]] auto const& node_dummy = nodePtr;
@@ -89,12 +89,12 @@ TEST_CASE("UrQMD") {
}
SECTION("nucleus projectile") {
- auto [env, csPtr, nodePtr] = setupEnvironment(particles::Code::Oxygen);
+ auto [env, csPtr, nodePtr] = testing::setupEnvironment(particles::Code::Oxygen);
[[maybe_unused]] auto const& env_dummy = env; // against warnings
[[maybe_unused]] auto const& node_dummy = nodePtr; // against warnings
unsigned short constexpr A = 14, Z = 7;
- auto [stackPtr, secViewPtr] = setupStack(A, Z, 400_GeV, nodePtr, *csPtr);
+ auto [stackPtr, secViewPtr] = setupStack(code::Nucleus, A, Z, 400_GeV, nodePtr, *csPtr);
REQUIRE(stackPtr->getEntries() == 1);
REQUIRE(secViewPtr->getEntries() == 0);
@@ -113,12 +113,12 @@ TEST_CASE("UrQMD") {
}
SECTION("\"special\" projectile") {
- auto [env, csPtr, nodePtr] = setupEnvironment(particles::Code::Oxygen);
+ auto [env, csPtr, nodePtr] = testing::setupEnvironment(particles::Code::Oxygen);
[[maybe_unused]] auto const& env_dummy = env; // against warnings
[[maybe_unused]] auto const& node_dummy = nodePtr; // against warnings
auto [stackPtr, secViewPtr] =
- setupStack(particles::Code::PiPlus, 400_GeV, nodePtr, *csPtr);
+ setupStack(particles::Code::PiPlus, 0,0, 400_GeV, nodePtr, *csPtr);
REQUIRE(stackPtr->getEntries() == 1);
REQUIRE(secViewPtr->getEntries() == 0);
@@ -139,12 +139,12 @@ TEST_CASE("UrQMD") {
}
SECTION("K0Long projectile") {
- auto [env, csPtr, nodePtr] = setupEnvironment(particles::Code::Oxygen);
+ auto [env, csPtr, nodePtr] = testing::setupEnvironment(particles::Code::Oxygen);
[[maybe_unused]] auto const& env_dummy = env; // against warnings
[[maybe_unused]] auto const& node_dummy = nodePtr; // against warnings
auto [stackPtr, secViewPtr] =
- setupStack(particles::Code::K0Long, 400_GeV, nodePtr, *csPtr);
+ setupStack(particles::Code::K0Long,0,0, 400_GeV, nodePtr, *csPtr);
REQUIRE(stackPtr->getEntries() == 1);
REQUIRE(secViewPtr->getEntries() == 0);
diff --git a/Setup/SetupStack.h b/Setup/SetupStack.h
index b54911bb0..91aa03417 100644
--- a/Setup/SetupStack.h
+++ b/Setup/SetupStack.h
@@ -147,28 +147,29 @@ namespace corsika::setup::testing {
using namespace corsika::units::si;
auto stack = std::make_unique<setup::Stack>();
- auto constexpr mN = corsika::units::constants::nucleonMass;
geometry::Point const origin(cs, {0_m, 0_m, 0_m});
corsika::stack::MomentumVector const pLab(cs, {vMomentum, 0_GeV, 0_GeV});
-
+
if (vProjectileType == particles::Code::Nucleus) {
+ auto constexpr mN = corsika::units::constants::nucleonMass;
HEPEnergyType const E0 = sqrt(units::static_pow<2>(mN * vA) + pLab.squaredNorm());
auto particle = stack->AddParticle(
std::make_tuple(particles::Code::Nucleus, E0, pLab, origin, 0_ns, vA, vZ));
+ particle.SetNode(vNodePtr);
+ return std::make_tuple(
+ std::move(stack),
+ std::make_unique<decltype(corsika::stack::SecondaryView(particle))>(particle));
} else { // not a nucleus
HEPEnergyType const E0 = sqrt(
units::static_pow<2>(particles::GetMass(vProjectileType)) + pLab.squaredNorm());
- auto particle = stack->AddParticle(
- std::tuple<particles::Code, units::si::HEPEnergyType,
- corsika::stack::MomentumVector, geometry::Point,
- units::si::TimeType>{vProjectileType, E0, pLab, origin, 0_ns});
+ auto particle =
+ stack->AddParticle(std::make_tuple(vProjectileType, E0, pLab, origin, 0_ns));
+ particle.SetNode(vNodePtr);
+ return std::make_tuple(
+ std::move(stack),
+ std::make_unique<decltype(corsika::stack::SecondaryView(particle))>(particle));
}
-
- particle.SetNode(vNodePtr);
- return std::make_tuple(
- std::move(stack),
- std::make_unique<decltype(corsika::stack::SecondaryView(particle))>(particle));
}
} // namespace corsika::setup::testing
--
GitLab