diff --git a/corsika/detail/modules/pythia8/Interaction.inl b/corsika/detail/modules/pythia8/Interaction.inl
index b2580ba752790441019fbaa2eeb669de7ac721e2..01e88c6d2700c2d3a9e68fff8acf4d935b79c877 100644
--- a/corsika/detail/modules/pythia8/Interaction.inl
+++ b/corsika/detail/modules/pythia8/Interaction.inl
@@ -43,8 +43,10 @@ namespace corsika::pythia8 {
// projectile and target init to p Nitrogen to initialize pythia with angantyr
pythia_.readString("Beams:allowIDASwitch = on");
- pythia_.settings.mode("Beams:idA", 2212);
- pythia_.settings.mode("Beams:idB", 1000070140);
+ pythia_.settings.mode("Beams:idA",
+ static_cast<PDGCodeIntType>(get_PDG(Code::Proton)));
+ pythia_.settings.mode("Beams:idB",
+ static_cast<PDGCodeIntType>(get_PDG(Code::Nitrogen)));
pythia_.readString("Beams:allowVariableEnergy = on");
pythia_.readString("Beams:frameType = 1"); // CoM
@@ -72,7 +74,7 @@ namespace corsika::pythia8 {
pythia_.readString("HadronLevel:Decay = off");
// Reduce printout and relax energy-momentum conservation.
- // pythia_.readString("Print:quiet = on");
+ pythia_.readString("Print:quiet = on");
pythia_.readString("Check:epTolErr = 0.1");
pythia_.readString("Check:epTolWarn = 0.0001");
pythia_.readString("Check:mTolErr = 0.01");
@@ -152,10 +154,11 @@ namespace corsika::pythia8 {
static_pow<2>(sqrtS), get_mass(projectileId), get_mass(targetId));
if (labE < eKinMinLab_) return false;
- bool const validProjectile = std::find(validProjectiles_.begin(), validProjectiles_.end(), projectileId) !=
- validProjectiles_.end();
- bool const validTarget = std::find(validTargets_.begin(), validTargets_.end(), targetId) !=
- validTargets_.end();
+ bool const validProjectile =
+ std::find(validProjectiles_.begin(), validProjectiles_.end(), projectileId) !=
+ validProjectiles_.end();
+ bool const validTarget = std::find(validTargets_.begin(), validTargets_.end(),
+ targetId) != validTargets_.end();
return validProjectile && validTarget;
}
@@ -167,9 +170,14 @@ namespace corsika::pythia8 {
Interaction::getCrossSectionInelEla(Code const projectileId, Code const targetId,
FourMomentum const& projectileP4,
FourMomentum const& targetP4) const {
+ // elastic cross sections only available for proton
if (projectileId != Code::Proton || targetId != Code::Proton) {
return std::tuple{CrossSectionType::zero(), CrossSectionType::zero()};
}
+ auto const sqrtS2 = (projectileP4 + targetP4).getNormSqr();
+ HEPEnergyType const sqrtS = sqrt(sqrtS2);
+ if (!isValid(projectileId, targetId, sqrtS))
+ return std::tuple{CrossSectionType::zero(), CrossSectionType::zero()};
HEPEnergyType const Elab =
calculate_lab_energy((projectileP4 + targetP4).getNormSqr(),
@@ -184,17 +192,17 @@ namespace corsika::pythia8 {
Code const targetId,
FourMomentum const& projectileP4,
FourMomentum const& targetP4) const {
-
+
// define system (this is Nucleus-Nucleus!!)
auto const sqrtS2 = (projectileP4 + targetP4).getNormSqr();
- HEPEnergyType const sqrtS = sqrt(sqrtS2);
- if(!isValid(projectileId, targetId, sqrtS)) return CrossSectionType::zero();
-
+ HEPEnergyType const sqrtS = sqrt(sqrtS2);
+ if (!isValid(projectileId, targetId, sqrtS)) return CrossSectionType::zero();
+
auto const it = xs_map_.find(projectileId);
auto const mapped_projectile = (it == xs_map_.end()) ? projectileId : it->second;
auto const& table = crossSectionTables_.at(std::pair{mapped_projectile, targetId});
HEPEnergyType const Elab =
- calculate_lab_energy(sqrtS2,get_mass(projectileId), get_mass(targetId));
+ calculate_lab_energy(sqrtS2, get_mass(projectileId), get_mass(targetId));
return table.interpolate(Elab);
}
diff --git a/corsika/modules/pythia8/Interaction.hpp b/corsika/modules/pythia8/Interaction.hpp
index 2fc27cc1c98fed35813f1195823cd7def8335e61..15d07e97cdede2ab5322ccde55584cc77fefea75 100644
--- a/corsika/modules/pythia8/Interaction.hpp
+++ b/corsika/modules/pythia8/Interaction.hpp
@@ -120,7 +120,8 @@ namespace corsika::pythia8 {
crossSectionPPInelastic_;
bool const print_listing_ = false;
- HEPEnergyType const eMaxLab_ = 1e21_eV; // Cross-section tables tabulated up to 10^21 eV
+ HEPEnergyType const eMaxLab_ =
+ 1e21_eV; // Cross-section tables tabulated up to 10^21 eV
HEPEnergyType const eKinMinLab_ = 0.2_GeV;
};
diff --git a/tests/modules/testPythia8Interaction.inl b/tests/modules/testPythia8Interaction.inl
index f38917ead720deef0c2e33d0a9505704afbb1376..991cbfd970657bf27249f94ed29003d4ea3efd83 100644
--- a/tests/modules/testPythia8Interaction.inl
+++ b/tests/modules/testPythia8Interaction.inl
@@ -17,18 +17,17 @@
#include "corsika/framework/core/Logging.hpp"
#include "corsika/framework/core/PhysicalUnits.hpp"
- logging::set_level(logging::level::debug);
+logging::set_level(logging::level::debug);
- // this will be a p-p collision at sqrts=3.5TeV -> no problem for pythia
- auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
- Code::Proton, 7_TeV, (DummyEnvironment::BaseNodeType* const)nodePtr, *csPtr);
- auto& view = *secViewPtr;
+// this will be a p-p collision at sqrts=3.5TeV -> no problem for pythia
+auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
+ Code::Proton, 7_TeV, (DummyEnvironment::BaseNodeType* const)nodePtr, *csPtr);
+auto& view = *secViewPtr;
- corsika::pythia8::Interaction collision;
+corsika::pythia8::Interaction collision;
SECTION("pythia interaction configurations") {
-
REQUIRE(collision.canInteract(Code::Proton));
REQUIRE(collision.canInteract(Code::AntiProton));
REQUIRE(collision.canInteract(Code::Neutron));
@@ -44,8 +43,8 @@ SECTION("pythia interaction: proton") {
// test some combinations of valid target and projectile particles
// so far only hadron-hadron and hadron-Nucleus is allowed
Code const target = GENERATE(Code::Proton, Code::Nitrogen, Code::Oxygen, Code::Argon);
- Code const projectile = Code::Proton;
-
+ Code const projectile = Code::Proton;
+
CORSIKA_LOG_INFO("testing: {}-{}", projectile, target);
REQUIRE(
collision.getCrossSection(
@@ -64,8 +63,8 @@ SECTION("pythia interaction: PiPlus") {
// test some combinations of valid target and projectile particles
// so far only hadron-hadron and hadron-Nucleus is allowed
Code const target = GENERATE(Code::Proton, Code::Nitrogen, Code::Oxygen, Code::Argon);
- Code const projectile = Code::PiPlus;
-
+ Code const projectile = Code::PiPlus;
+
CORSIKA_LOG_INFO("testing: {}-{}", projectile, target);
REQUIRE(
collision.getCrossSection(
@@ -84,8 +83,8 @@ SECTION("pythia interaction: KPlus") {
// test some combinations of valid target and projectile particles
// so far only hadron-hadron and hadron-Nucleus is allowed
Code const target = GENERATE(Code::Proton, Code::Nitrogen, Code::Oxygen, Code::Argon);
- Code const projectile = Code::KPlus;
-
+ Code const projectile = Code::KPlus;
+
CORSIKA_LOG_INFO("testing: {}-{}", projectile, target);
REQUIRE(
collision.getCrossSection(
@@ -180,9 +179,11 @@ SECTION("pythia wrong projectile") {
{Nitrogen::mass, {rootCS, {0_eV, 0_eV, 0_eV}}}));
// gamma+p not possible
- REQUIRE(collision.getCrossSection(
- Code::H0, Code::Proton, {calculate_total_energy(100_GeV, H0::mass), {rootCS, {0_eV, 0_eV, 100_GeV}}},
- {Proton::mass, {rootCS, {0_eV, 0_eV, 0_eV}}}) == CrossSectionType::zero());
+ REQUIRE(
+ collision.getCrossSection(
+ Code::H0, Code::Proton,
+ {calculate_total_energy(100_GeV, H0::mass), {rootCS, {0_eV, 0_eV, 100_GeV}}},
+ {Proton::mass, {rootCS, {0_eV, 0_eV, 0_eV}}}) == CrossSectionType::zero());
REQUIRE(collision.getCrossSectionInelEla(
Code::Photon, Code::Proton, {100_GeV, {rootCS, {0_eV, 0_eV, 100_GeV}}},