diff --git a/corsika/detail/modules/pythia8/Interaction.inl b/corsika/detail/modules/pythia8/Interaction.inl index 6a62a468fbb4e5d02bf2c8946457f31d57fa99c2..98081cc0fbea41b7bcec7f0966d7f07e31737656 100644 --- a/corsika/detail/modules/pythia8/Interaction.inl +++ b/corsika/detail/modules/pythia8/Interaction.inl @@ -146,16 +146,22 @@ namespace corsika::pythia8 { } inline bool Interaction::isValid(Code const projectileId, Code const targetId, - HEPEnergyType const sqrtS) const { - CORSIKA_LOG_DEBUG("pythia isValid: {}+{} at sqrtS = {} GeV", projectileId, targetId, - sqrtS / 1_GeV); + HEPEnergyType const sqrtSNN) const { + + CORSIKA_LOG_DEBUG("pythia isValid: {} + {} at sqrtSNN = {} GeV", projectileId, + targetId, sqrtSNN / 1_GeV); + if (is_nucleus(targetId)) + CORSIKA_LOG_DEBUG("nucleus = {}-{}", get_nucleus_A(targetId), + get_nucleus_Z(targetId)); if (is_nucleus(projectileId)) // not yet possible with Pythia return false; if (!canInteract(projectileId)) return false; - HEPEnergyType const labE = calculate_lab_energy( - static_pow<2>(sqrtS), get_mass(projectileId), get_mass(targetId)); + auto const mass_target = + get_mass(targetId) / (is_nucleus(targetId) ? get_nucleus_A(targetId) : 1.); + HEPEnergyType const labE = + calculate_lab_energy(static_pow<2>(sqrtSNN), get_mass(projectileId), mass_target); if (labE < eKinMinLab_) return false; bool const validProjectile = @@ -220,7 +226,7 @@ namespace corsika::pythia8 { CORSIKA_LOG_DEBUG( "Pythia::Interaction: " - "doInteraction: {} interaction? ", + "doInteraction: {} interaction? {} ", projectileId, corsika::pythia8::Interaction::canInteract(projectileId)); // define system (this is Nucleus-Nucleus!!) @@ -230,10 +236,6 @@ namespace corsika::pythia8 { HEPEnergyType const eProjectileLab = calculate_lab_energy( sqrtS2NN, get_mass(projectileId), get_mass(targetId) / get_nucleus_A(targetId)); - if (!isValid(projectileId, targetId, sqrtSNN)) { - throw std::runtime_error("invalid target,projectile,energy combination."); - } - CORSIKA_LOG_DEBUG("Interaction: ebeam lab: {} GeV", eProjectileLab / 1_GeV); CORSIKA_LOG_DEBUG( "Interaction: " @@ -241,6 +243,10 @@ namespace corsika::pythia8 { " Ecm_NN(GeV): {}", eProjectileLab / 1_GeV, sqrtSNN / 1_GeV); + if (!isValid(projectileId, targetId, sqrtSNN)) { + throw std::runtime_error("invalid target,projectile,energy combination."); + } + COMBoost const labFrameBoost{targetP4.getSpaceLikeComponents(), get_mass(targetId)}; // auto const proj4MomLab = labFrameBoost.toCoM(projectileP4); auto const& rotCS = labFrameBoost.getRotatedCS(); @@ -250,8 +256,8 @@ namespace corsika::pythia8 { double const eCM_pythia = sqrtSNN / 1_GeV; // CoM energy hadron-Nucleon double const idA_pythia = static_cast<double>(get_PDG(projectileId)); double const idB_pythia = static_cast<double>(get_PDG(targetId)); - CORSIKA_LOG_INFO("re-configuring pythia idA={}, idB={}, ecm={} GeV", idA_pythia, - idB_pythia, eCM_pythia); + CORSIKA_LOG_DEBUG("re-configuring pythia idA={}, idB={}, ecm={} GeV", idA_pythia, + idB_pythia, eCM_pythia); // configure event pythia_.setBeamIDs(idA_pythia, idB_pythia); pythia_.setKinematics(eCM_pythia); diff --git a/corsika/modules/pythia8/Interaction.hpp b/corsika/modules/pythia8/Interaction.hpp index 15d07e97cdede2ab5322ccde55584cc77fefea75..df5ece90c57f9f1672190d38388e78e05242ab4a 100644 --- a/corsika/modules/pythia8/Interaction.hpp +++ b/corsika/modules/pythia8/Interaction.hpp @@ -122,7 +122,7 @@ namespace corsika::pythia8 { bool const print_listing_ = false; HEPEnergyType const eMaxLab_ = 1e21_eV; // Cross-section tables tabulated up to 10^21 eV - HEPEnergyType const eKinMinLab_ = 0.2_GeV; + HEPEnergyType const eKinMinLab_ = 100_GeV; }; } // namespace corsika::pythia8 diff --git a/tests/modules/testPythia8Interaction.inl b/tests/modules/testPythia8Interaction.inl index 1ccc00fcc7a1a26a03f331b8705614e8ee14415f..77c3643a3e63e874e2cfed0b5a76dd22874c33f8 100644 --- a/tests/modules/testPythia8Interaction.inl +++ b/tests/modules/testPythia8Interaction.inl @@ -81,7 +81,7 @@ SECTION("pythia cross section bug") { Code const target = Code::Oxygen; Code const projectile = Code::Proton; - corsika::units::si::HEPMomentumType P1 = 80_GeV; + corsika::units::si::HEPMomentumType P1 = 100_GeV; CORSIKA_LOG_INFO("testing: {} - {}", projectile, target); REQUIRE( @@ -101,11 +101,10 @@ SECTION("pythia too low energy") { corsika::pythia8::Interaction collision; // 5 MeV lab is too low, 0 mb expected - REQUIRE_THROWS( - collision.getCrossSectionInelEla( - Code::Proton, Code::Proton, - {calculate_total_energy(Proton::mass, 5_MeV), {rootCS, 0_eV, 0_eV, 5_MeV}}, - {Proton::mass, {rootCS, 0_eV, 0_eV, 0_eV}}) == std::tuple{0_mb, 0_mb}); + REQUIRE(collision.getCrossSectionInelEla( + Code::Proton, Code::Proton, + {calculate_total_energy(Proton::mass, 5_MeV), {rootCS, 0_eV, 0_eV, 5_MeV}}, + {Proton::mass, {rootCS, 0_eV, 0_eV, 0_eV}}) == std::tuple{0_mb, 0_mb}); REQUIRE_THROWS(collision.doInteraction( view, Code::Neutron, Code::Proton,