IAP GITLAB

Skip to content
Snippets Groups Projects
Commit 9cc0c321 authored by ralfulrich's avatar ralfulrich Committed by Ralf Ulrich
Browse files

cleanup of process tests

parent 941f15a8
No related branches found
No related tags found
1 merge request!265Add medium type, as new property (air, water, rock...)
Hide whitespace changes
Inline Side-by-side
Processes/Pythia/testPythia8.cc
+ 9
− 46
View file @ 9cc0c321
......@@ -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;
......
Processes/QGSJetII/testQGSJetII.cc
+ 6
− 34
View file @ 9cc0c321
......@@ -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();
......
Processes/Sibyll/testSibyll.cc
+ 13
− 83
View file @ 9cc0c321
......@@ -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);
......
Processes/UrQMD/testUrQMD.cc
+ 7
− 7
View file @ 9cc0c321
......@@ -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);
......
Setup/SetupStack.h
+ 12
− 11
View file @ 9cc0c321
......@@ -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
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

data privacy - impressum