IAP GITLAB

Skip to content
Snippets Groups Projects
Commit 1143a7bd authored by Maximilian Reininghaus's avatar Maximilian Reininghaus :vulcan: Committed by Maximilian Reininghaus
Browse files

re-enabled test-cases that should work

parent 0cfb67a0
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
corsika/detail/modules/qgsjetII/InteractionModel.inl
+ 20
− 18
View file @ 1143a7bd
......@@ -28,7 +28,7 @@ namespace corsika::qgsjetII {
// initialize QgsjetII
static bool initialized = false;
if (!initialized) {
CORSIKA_LOG_DEBUG("Reading QGSJetII data tables from {}", dataPath);
CORSIKA_LOG_DEBUG("Reading QGSJetII data tables from {}", dataPath);
qgset_();
datadir DIR(dataPath.string() + "/");
qgaini_(DIR.data);
......@@ -87,7 +87,7 @@ namespace corsika::qgsjetII {
int const iBeam = static_cast<QgsjetIIXSClassIntType>(
corsika::qgsjetII::getQgsjetIIXSCode(projectileId));
CORSIKA_LOG_DEBUG(
"QgsjetII::getCrossSection Elab= {} GeV iBeam= {}"
" iProjectile= {} iTarget= {}",
......@@ -121,15 +121,17 @@ namespace corsika::qgsjetII {
if (!corsika::qgsjetII::canInteract(projectileId) ||
!isValid(projectileId, targetId, sqrtSNN)) {
throw std::runtime_error(fmt::format("invalid target [{}]/ projectile [{}] /energy [{} GeV] combination.", get_name(targetId, full_name{}), get_name(projectileId, full_name{}), sqrtSNN/1_GeV));
throw std::runtime_error(fmt::format(
"invalid target [{}]/ projectile [{}] /energy [{} GeV] combination.",
get_name(targetId, full_name{}), get_name(projectileId, full_name{}),
sqrtSNN / 1_GeV));
}
auto const projMass = get_mass(projectileId);
auto const targetMass = get_mass(targetId);
// lab-frame energy per projectile nucleon
HEPEnergyType const Elab =
calculate_lab_energy(S, projMass, targetMass);
HEPEnergyType const Elab = calculate_lab_energy(S, projMass, targetMass);
auto const ElabN = Elab / AfactorProjectile;
CORSIKA_LOG_DEBUG("ebeam lab: {} GeV per projectile nucleon", ElabN / 1_GeV);
......@@ -142,7 +144,8 @@ namespace corsika::qgsjetII {
QgsjetIIHadronType qgsjet_hadron_type = qgsjetII::getQgsjetIIHadronType(projectileId);
if (qgsjet_hadron_type == QgsjetIIHadronType::NucleusType) {
projectileMassNumber = get_nucleus_A(projectileId);
qgsjet_hadron_type = bernoulli_(rng_) ? QgsjetIIHadronType::ProtonType : QgsjetIIHadronType::NeutronType;
qgsjet_hadron_type = bernoulli_(rng_) ? QgsjetIIHadronType::ProtonType
: QgsjetIIHadronType::NeutronType;
} else if (qgsjet_hadron_type == QgsjetIIHadronType::NeutralLightMesonType) {
// from conex: replace pi0 or rho0 with pi+/pi- in alternating sequence
qgsjet_hadron_type = alternate_;
......@@ -152,7 +155,8 @@ namespace corsika::qgsjetII {
}
count_++;
int const qgsjet_hadron_type_int = static_cast<QgsjetIICodeIntType>(qgsjet_hadron_type);
int const qgsjet_hadron_type_int =
static_cast<QgsjetIICodeIntType>(qgsjet_hadron_type);
CORSIKA_LOG_DEBUG(
"qgsjet_hadron_type_int={} projectileMassNumber={} targetMassNumber={}",
qgsjet_hadron_type_int, projectileMassNumber, targetMassNumber);
......@@ -195,17 +199,18 @@ namespace corsika::qgsjetII {
HEPMassType const nucleonMass = get_mass(idFragm);
// no pT, fragments just go forward
MomentumVector const momentum{csPrime, {0_eV, 0_eV, calculate_momentum(ElabN, nucleonMass)}};
MomentumVector const momentum{
csPrime, {0_eV, 0_eV, calculate_momentum(ElabN, nucleonMass)}};
// this is not "CoM" here, but rather the system defined by projectile+target,
// which in Cascade-mode is already lab
auto const P4com =
boostInternal.toCoM(FourVector{ElabN, momentum});
auto const P4com = boostInternal.toCoM(FourVector{ElabN, momentum});
auto const P4output = boost.fromCoM(P4com);
auto p3output = P4output.getSpaceLikeComponents();
p3output.rebase(originalCS); // transform back into standard lab frame
HEPEnergyType const Ekin = calculate_kinetic_energy(p3output.getNorm(), nucleonMass);
HEPEnergyType const Ekin =
calculate_kinetic_energy(p3output.getNorm(), nucleonMass);
CORSIKA_LOG_DEBUG(
"secondary fragment> id= {}"
......@@ -239,15 +244,12 @@ namespace corsika::qgsjetII {
Code const idFragm = get_nucleus_code(A, Z);
HEPEnergyType const mass = get_mass(idFragm);
// no pT, frgments just go forward
MomentumVector momentum{
csPrime,
{0.0_GeV, 0.0_GeV,
calculate_momentum(ElabN * A, mass)}};
MomentumVector momentum{csPrime,
{0.0_GeV, 0.0_GeV, calculate_momentum(ElabN * A, mass)}};
// this is not "CoM" here, but rather the system defined by projectile+target,
// which in Cascade-mode is already lab
auto const P4com =
boostInternal.toCoM(FourVector{ElabN * A, momentum});
auto const P4com = boostInternal.toCoM(FourVector{ElabN * A, momentum});
auto const P4output = boost.fromCoM(P4com);
auto p3output = P4output.getSpaceLikeComponents();
p3output.rebase(originalCS); // transform back into standard lab frame
......
corsika/framework/utility/COMBoost.hpp
+ 2
− 1
View file @ 1143a7bd
......@@ -56,7 +56,8 @@ namespace corsika {
COMBoost(FourMomentum const& P4projectile, HEPEnergyType const massTarget);
/**
* Construct a COMBoost to boost into the rest frame of a particle given its 3-momentum and mass.
* Construct a COMBoost to boost into the rest frame of a particle given its
* 3-momentum and mass.
*/
COMBoost(MomentumVector const& momentum, HEPEnergyType const mass);
......
tests/modules/testQGSJetII.cpp
+ 110
− 105
View file @ 1143a7bd
......@@ -142,7 +142,8 @@ TEST_CASE("QgsjetIIInterface", "interaction,processes") {
corsika::qgsjetII::InteractionModel model;
SECTION("cross-sections") {
auto projCode = GENERATE(Code::PiPlus, Code::Proton, Code::K0Long, Code::Nitrogen, Code::Helium);
auto projCode =
GENERATE(Code::PiPlus, Code::Proton, Code::K0Long, Code::Nitrogen, Code::Helium);
auto targetCode = GENERATE(Code::Oxygen, Code::Nitrogen);
auto projEnergy = GENERATE(1_PeV, 1e18_eV);
......@@ -152,24 +153,28 @@ TEST_CASE("QgsjetIIInterface", "interaction,processes") {
REQUIRE(model.getCrossSection(
projCode, targetCode, FourMomentum{projEnergy, projMomentum},
FourMomentum{get_mass(targetCode), {*csPtr, 0_eV, 0_eV, 0_eV}}) /
1_mb > 0);
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 = GENERATE(1_PeV); //, 1e20_eV);
auto projCode =
GENERATE(Code::PiPlus, Code::Proton, Code::K0Long, Code::Nitrogen, Code::Helium);
auto targetCode = GENERATE(Code::Oxygen, Code::Nitrogen);
auto projMomentum = GENERATE(100_GeV, 1_PeV, 1e20_eV);
auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
Code::Proton, projMomentum, (DummyEnvironment::BaseNodeType* const)nodePtr, *csPtr);
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}}});
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
......@@ -186,99 +191,99 @@ TEST_CASE("QgsjetIIInterface", "interaction,processes") {
Approx(0).margin(1e-2));
}
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(150).margin(150)); // this is not physics validation
int countFragments = 0;
for (auto const& sec : view) { countFragments += (is_nucleus(sec.getPID())); }
CHECK(countFragments == Approx(4).margin(2)); // 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) /
1_mb ==
Approx(0));
CHECK(model.getCrossSection(Code::Nucleus, Code::Nucleus, aP4, bP4) / 1_mb ==
Approx(0));
CHECK_THROWS(
model.doInteraction(view, get_nucleus_code(1000, 500), Code::Oxygen, aP4, bP4));
}
SECTION("Allowed Particles") {
{ // pi0 is internally converted into pi+/pi-
auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
Code::Pi0, 1000_GeV, (DummyEnvironment::BaseNodeType* const)nodePtr, *csPtr);
[[maybe_unused]] test::StackView& view = *(secViewPtr.get());
[[maybe_unused]] auto particle = stackPtr->first();
model.doInteraction(view, Code::Pi0, Code::Oxygen,
{sqrt(static_pow<2>(1_TeV) + static_pow<2>(Pi0::mass)),
MomentumVector{cs, 1_TeV, 0_GeV, 0_GeV}},
{Oxygen::mass, MomentumVector{cs, 0_eV, 0_eV, 0_eV}});
CHECK(view.getSize() == Approx(20).margin(20)); // this is not physics validation
}
{ // rho0 is internally converted into pi-/pi+
auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
Code::Rho0, 1000_GeV, (DummyEnvironment::BaseNodeType* const)nodePtr, *csPtr);
[[maybe_unused]] test::StackView& view = *(secViewPtr.get());
[[maybe_unused]] auto particle = stackPtr->first();
model.doInteraction(view, Code::Rho0, Code::Oxygen,
{sqrt(static_pow<2>(1_TeV) + static_pow<2>(Rho0::mass)),
MomentumVector{cs, 1_TeV, 0_GeV, 0_GeV}},
{Oxygen::mass, MomentumVector{cs, 0_eV, 0_eV, 0_eV}});
CHECK(view.getSize() == Approx(50).margin(50)); // this is not physics validation
}
{ // Lambda is internally converted into neutron
auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
Code::Lambda0, 100_GeV, (DummyEnvironment::BaseNodeType* const)nodePtr, *csPtr);
[[maybe_unused]] test::StackView& view = *(secViewPtr.get());
[[maybe_unused]] auto particle = stackPtr->first();
model.doInteraction(view, Code::Lambda0, Code::Oxygen,
{sqrt(static_pow<2>(100_GeV) + static_pow<2>(Lambda0::mass)),
MomentumVector{cs, 100_GeV, 0_GeV, 0_GeV}},
{Oxygen::mass, MomentumVector{cs, 0_eV, 0_eV, 0_eV}});
CHECK(view.getSize() == Approx(50).margin(50)); // this is not physics validation
}
{ // AntiLambda is internally converted into anti neutron
auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
Code::Lambda0Bar, 1000_GeV, (DummyEnvironment::BaseNodeType* const)nodePtr,
*csPtr);
[[maybe_unused]] test::StackView& view = *(secViewPtr.get());
[[maybe_unused]] auto particle = stackPtr->first();
model.doInteraction(view, Code::Lambda0Bar, Code::Oxygen,
{sqrt(static_pow<2>(1_TeV) + static_pow<2>(Lambda0Bar::mass)),
MomentumVector{cs, 1_TeV, 0_GeV, 0_GeV}},
{Oxygen::mass, MomentumVector{cs, 0_eV, 0_eV, 0_eV}});
CHECK(view.getSize() == Approx(70).margin(67)); // this is not physics validation
}
}
//~ 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(150).margin(150)); // this is not physics validation
//~ int countFragments = 0;
//~ for (auto const& sec : view) { countFragments += (is_nucleus(sec.getPID())); }
//~ CHECK(countFragments == Approx(4).margin(2)); // 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) /
//~ 1_mb ==
//~ Approx(0));
//~ CHECK(model.getCrossSection(Code::Nucleus, Code::Nucleus, aP4, bP4) / 1_mb ==
//~ Approx(0));
//~ CHECK_THROWS(
//~ model.doInteraction(view, get_nucleus_code(1000, 500), Code::Oxygen, aP4, bP4));
//~ }
//~ SECTION("Allowed Particles") {
//~ { // pi0 is internally converted into pi+/pi-
//~ auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
//~ Code::Pi0, 1000_GeV, (DummyEnvironment::BaseNodeType* const)nodePtr, *csPtr);
//~ [[maybe_unused]] test::StackView& view = *(secViewPtr.get());
//~ [[maybe_unused]] auto particle = stackPtr->first();
//~ model.doInteraction(view, Code::Pi0, Code::Oxygen,
//~ {sqrt(static_pow<2>(1_TeV) + static_pow<2>(Pi0::mass)),
//~ MomentumVector{cs, 1_TeV, 0_GeV, 0_GeV}},
//~ {Oxygen::mass, MomentumVector{cs, 0_eV, 0_eV, 0_eV}});
//~ CHECK(view.getSize() == Approx(20).margin(20)); // this is not physics validation
//~ }
//~ { // rho0 is internally converted into pi-/pi+
//~ auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
//~ Code::Rho0, 1000_GeV, (DummyEnvironment::BaseNodeType* const)nodePtr, *csPtr);
//~ [[maybe_unused]] test::StackView& view = *(secViewPtr.get());
//~ [[maybe_unused]] auto particle = stackPtr->first();
//~ model.doInteraction(view, Code::Rho0, Code::Oxygen,
//~ {sqrt(static_pow<2>(1_TeV) + static_pow<2>(Rho0::mass)),
//~ MomentumVector{cs, 1_TeV, 0_GeV, 0_GeV}},
//~ {Oxygen::mass, MomentumVector{cs, 0_eV, 0_eV, 0_eV}});
//~ CHECK(view.getSize() == Approx(50).margin(50)); // this is not physics validation
//~ }
//~ { // Lambda is internally converted into neutron
//~ auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
//~ Code::Lambda0, 100_GeV, (DummyEnvironment::BaseNodeType* const)nodePtr, *csPtr);
//~ [[maybe_unused]] test::StackView& view = *(secViewPtr.get());
//~ [[maybe_unused]] auto particle = stackPtr->first();
//~ model.doInteraction(view, Code::Lambda0, Code::Oxygen,
//~ {sqrt(static_pow<2>(100_GeV) + static_pow<2>(Lambda0::mass)),
//~ MomentumVector{cs, 100_GeV, 0_GeV, 0_GeV}},
//~ {Oxygen::mass, MomentumVector{cs, 0_eV, 0_eV, 0_eV}});
//~ CHECK(view.getSize() == Approx(50).margin(50)); // this is not physics validation
//~ }
//~ { // AntiLambda is internally converted into anti neutron
//~ auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
//~ Code::Lambda0Bar, 1000_GeV, (DummyEnvironment::BaseNodeType* const)nodePtr,
//~ *csPtr);
//~ [[maybe_unused]] test::StackView& view = *(secViewPtr.get());
//~ [[maybe_unused]] auto particle = stackPtr->first();
//~ model.doInteraction(view, Code::Lambda0Bar, Code::Oxygen,
//~ {sqrt(static_pow<2>(1_TeV) + static_pow<2>(Lambda0Bar::mass)),
//~ MomentumVector{cs, 1_TeV, 0_GeV, 0_GeV}},
//~ {Oxygen::mass, MomentumVector{cs, 0_eV, 0_eV, 0_eV}});
//~ CHECK(view.getSize() == Approx(70).margin(67)); // this is not physics validation
//~ }
//~ }
}
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