/*
* (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
*
* This software is distributed under the terms of the 3-clause BSD license.
* See file LICENSE for a full version of the license.
*/
#include <corsika/modules/QGSJetII.hpp>
#include <corsika/framework/core/ParticleProperties.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/geometry/Point.hpp>
#include <corsika/framework/random/RNGManager.hpp>
#include <SetupTestEnvironment.hpp>
#include <catch2/catch_all.hpp>
#include <string>
#include <cstdlib>
using namespace corsika;
using Catch::Approx;
using DummyEnvironmentInterface = IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>;
using DummyEnvironment = Environment<DummyEnvironmentInterface>;
template <typename TStackView>
auto sumCharge(TStackView const& view) {
int totalCharge = 0;
for (auto const& p : view) { totalCharge += get_charge_number(p.getPID()); }
return totalCharge;
}
template <typename TStackView>
auto sumMomentum(TStackView const& view, CoordinateSystemPtr const& vCS) {
Vector<hepenergy_d> sum{vCS, 0_eV, 0_eV, 0_eV};
for (auto const& p : view) { sum += p.getMomentum(); }
return sum;
}
TEST_CASE("QgsjetII", "[processes]") {
logging::set_level(logging::level::info);
RNGManager<>::getInstance().registerRandomStream("qgsjet");
SECTION("Corsika -> QgsjetII") {
CHECK(corsika::qgsjetII::convertToQgsjetII(PiMinus::code) ==
corsika::qgsjetII::QgsjetIICode::PiMinus);
CHECK(corsika::qgsjetII::convertToQgsjetIIRaw(Proton::code) == 2);
}
SECTION("QgsjetII -> Corsika") {
CHECK(Code::PiPlus == corsika::qgsjetII::convertFromQgsjetII(
corsika::qgsjetII::QgsjetIICode::PiPlus));
CHECK_THROWS(
corsika::qgsjetII::convertFromQgsjetII(corsika::qgsjetII::QgsjetIICode::Unknown));
}
SECTION("Corsika -> QgsjetII") {
CHECK(corsika::qgsjetII::convertToQgsjetII(Code::PiMinus) ==
corsika::qgsjetII::QgsjetIICode::PiMinus);
CHECK(corsika::qgsjetII::convertToQgsjetIIRaw(Code::Proton) == 2);
}
SECTION("canInteractInQgsjetII") {
CHECK(corsika::qgsjetII::canInteract(Code::Proton));
CHECK(corsika::qgsjetII::canInteract(Code::KPlus));
CHECK(corsika::qgsjetII::canInteract(Code::Nucleus));
CHECK(corsika::qgsjetII::canInteract(Code::Rho0));
// CHECK(corsika::qgsjetII::canInteract(Helium::getCode()));
CHECK_FALSE(corsika::qgsjetII::canInteract(Code::EtaC));
CHECK_FALSE(corsika::qgsjetII::canInteract(Code::SigmaC_2455_0));
}
SECTION("cross-section type") {
CHECK(corsika::qgsjetII::getQgsjetIIXSCode(Code::Neutron) ==
corsika::qgsjetII::QgsjetIIXSClass::Baryons);
CHECK(corsika::qgsjetII::getQgsjetIIXSCode(Code::K0Long) ==
corsika::qgsjetII::QgsjetIIXSClass::Kaons);
CHECK(corsika::qgsjetII::getQgsjetIIXSCode(Code::Proton) ==
corsika::qgsjetII::QgsjetIIXSClass::Baryons);
CHECK(corsika::qgsjetII::getQgsjetIIXSCode(Code::PiMinus) ==
corsika::qgsjetII::QgsjetIIXSClass::LightMesons);
CHECK(corsika::qgsjetII::getQgsjetIIXSCode(Code::Helium) ==
corsika::qgsjetII::QgsjetIIXSClass::Baryons);
}
SECTION("valid") {
corsika::qgsjetII::InteractionModel model;
CHECK_FALSE(model.isValid(Code::Electron, Code::Proton, 1_TeV));
CHECK_FALSE(model.isValid(Code::Proton, Code::Electron, 1_TeV));
CHECK_FALSE(model.isValid(Code::Proton, Code::Proton, 1_GeV));
CHECK(model.isValid(Code::Proton, Code::Helium, 1_TeV));
CHECK_FALSE(model.isValid(Code::Proton, Code::Helium, 1_GeV));
}
}
#include <corsika/framework/geometry/Point.hpp>
#include <corsika/framework/geometry/RootCoordinateSystem.hpp>
#include <corsika/framework/geometry/Vector.hpp>
#include <corsika/framework/core/ParticleProperties.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/setup/SetupStack.hpp>
#include <corsika/setup/SetupTrajectory.hpp>
#include <corsika/media/Environment.hpp>
#include <corsika/media/HomogeneousMedium.hpp>
#include <corsika/media/NuclearComposition.hpp>
#include <SetupTestEnvironment.hpp>
#include <SetupTestStack.hpp>
TEST_CASE("QgsjetIIInterface", "interaction,processes") {
logging::set_level(logging::level::info);
RNGManager<>::getInstance().registerRandomStream("qgsjet");
auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Oxygen);
auto const& cs = *csPtr;
[[maybe_unused]] auto const& env_dummy = env;
[[maybe_unused]] auto const& node_dummy = nodePtr;
corsika::qgsjetII::InteractionModel model;
SECTION("cross-sections") {
auto projCode = GENERATE(Code::PiPlus, Code::Proton, Code::K0Long, Code::Iron,
Code::Nitrogen, Code::Helium);
auto targetCode = GENERATE(Code::Oxygen, Code::Nitrogen);
auto projEnergy = GENERATE(1_PeV, 1e18_eV);
auto momMagnitude = calculate_momentum(projEnergy, get_mass(projCode));
MomentumVector const projMomentum{*csPtr, 0_eV, momMagnitude, 0_eV};
REQUIRE(model.getCrossSection(
projCode, targetCode, FourMomentum{projEnergy, projMomentum},
FourMomentum{get_mass(targetCode), {*csPtr, 0_eV, 0_eV, 0_eV}}) /
1_mb >
0);
}
SECTION("InteractionInterface") {
auto projCode = GENERATE(Code::PiPlus, Code::Proton, Code::K0Long, Code::Iron,
Code::Nitrogen, Code::Helium);
auto targetCode = GENERATE(Code::Oxygen, Code::Nitrogen);
auto projMomentum = is_nucleus(projCode)
? GENERATE(1e20_eV, 1_PeV, 100_TeV, 10_TeV)
: GENERATE(1e20_eV, 1_PeV, 100_TeV, 10_TeV, 1_TeV, 100_GeV);
auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
Code::Proton, projMomentum, (DummyEnvironment::BaseNodeType* const)nodePtr,
*csPtr);
test::StackView& view = *(secViewPtr.get());
auto projectile = secViewPtr->getProjectile();
auto const projectileMomentum = projectile.getMomentum();
model.doInteraction(
view, projCode, targetCode,
FourMomentum{calculate_total_energy(projMomentum, get_mass(projCode)),
projectileMomentum},
FourMomentum{get_mass(targetCode), MomentumVector{cs, {0_eV, 0_eV, 0_eV}}});
/* **********************************
As it turned out already twice (#291 and #307), the detailed output of
QGSJetII event generation depends on the gfortran version used. This is not reliable
and cannot be tested in a unit test here. One related problem was already found
(#291) and is realted to undefined behaviour in the evaluation of functions in logical
expressions. It is not clear if #307 is the same issue.
CHECK(view.getSize() == 14);
CHECK(sumCharge(view) == 2);
bit us again, lets set the size to max. size of stack in qgsjet
*********************************** */
auto const secMomSum = sumMomentum(view, projectileMomentum.getCoordinateSystem());
CHECK((secMomSum - projectileMomentum).getNorm() / projectileMomentum.getNorm() ==
Approx(0).margin(1e-2));
CHECK(view.getSize() == Approx(95000).margin(94998));
}
SECTION("InteractionInterface Nuclei") {
HEPEnergyType const P0 = 20100_GeV;
MomentumVector const plab = MomentumVector(cs, {P0, 0_eV, 0_eV});
Code const pid = get_nucleus_code(60, 30);
auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
pid, P0, (DummyEnvironment::BaseNodeType* const)nodePtr, *csPtr);
test::StackView& view = *(secViewPtr.get());
HEPEnergyType const Elab = sqrt(static_pow<2>(P0) + static_pow<2>(get_mass(pid)));
FourMomentum const projectileP4(Elab, plab);
FourMomentum const targetP4(Oxygen::mass, MomentumVector(cs, {0_eV, 0_eV, 0_eV}));
view.clear();
model.doInteraction(view, pid, Code::Oxygen, projectileP4,
targetP4); // this also should produce some fragments
CHECK(view.getSize() ==
Approx(95000).margin(94998)); // this is not physics validation
int countFragments = 0;
for (auto const& sec : view) { countFragments += (is_nucleus(sec.getPID())); }
CHECK(countFragments == Approx(4).margin(3)); // this is not physics validation
}
SECTION("Heavy nuclei") {
auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
get_nucleus_code(1000, 1000), 1100_GeV,
(DummyEnvironment::BaseNodeType* const)nodePtr, *csPtr);
test::StackView& view = *(secViewPtr.get());
auto projectile = secViewPtr->getProjectile();
auto const projectileMomentum = projectile.getMomentum();
FourMomentum const aP4(100_GeV, {cs, 99_GeV, 0_GeV, 0_GeV});
FourMomentum const bP4(1_TeV, {cs, 0.9_TeV, 0_GeV, 0_GeV});
CHECK(model.getCrossSection(get_nucleus_code(10, 5), get_nucleus_code(1000, 500), aP4,
bP4) == 0_mb);
CHECK(model.getCrossSection(Code::Nucleus, Code::Nucleus, aP4, bP4) == 0_mb);
CHECK_THROWS(
model.doInteraction(view, get_nucleus_code(1000, 500), Code::Oxygen, aP4, bP4));
}
SECTION("Allowed Particles") {
HEPEnergyType const projMomentum = 500_GeV;
// pi0 is internally converted into pi+/pi-
// rho0 is internally converted into pi-/pi+
// (anti-)Lambda is internally converted into (anti-)neutron
auto pid = GENERATE(Code::Pi0, Code::Rho0, Code::Lambda, Code::LambdaBar);
auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
Code::Proton, projMomentum, (DummyEnvironment::BaseNodeType* const)nodePtr,
*csPtr);
test::StackView& view = *(secViewPtr.get());
auto projectile = secViewPtr->getProjectile();
auto const projectileMomentum = projectile.getMomentum();
CHECK_NOTHROW(model.doInteraction(
view, pid, Code::Oxygen,
FourMomentum{calculate_total_energy(projMomentum, get_mass(pid)),
projectileMomentum},
FourMomentum{get_mass(Code::Oxygen), MomentumVector{cs, {0_eV, 0_eV, 0_eV}}}));
CHECK(view.getSize() == Approx(70).margin(69));
}
}