diff --git a/Processes/Pythia/testPythia8.cc b/Processes/Pythia/testPythia8.cc index 4fd221bd66cdb45378d0f21bd162336d97c82084..a0af2a46a138e92c7d14032acb17400f5ebf8ca6 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 d1e9c7783fd1a57105e01d6d3837b907a36b59ad..7f9af0675f759dce9b3b2a5293c062e45ca4503f 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 755237efa5d346f133bd0b9f999322e35ba4c5fd..ccdb4b651482f85fb1ef5b63f90f0c73c5eb731d 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 537b656c664c3b4f98d012f8b123bc41b34b4f72..7bd929e3ec8c9177cb36af276fc9f4aaac5509f2 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 b54911bb03e146ab352b08e401cd6e05bfab07c6..91aa034179d2402e91bbdc429b6a711f347f9b1c 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