IAP GITLAB

Skip to content
Snippets Groups Projects
Commit a54a72fc authored by Felix Riehn's avatar Felix Riehn
Browse files

set threshold in pythia to 100GeV, reduce verbosity, make debug output more sensible

parent 4935690a
No related branches found
No related tags found
1 merge request!549Draft: "Use Pythia for interactions"
Pipeline #13645 passed
...@@ -146,16 +146,22 @@ namespace corsika::pythia8 { ...@@ -146,16 +146,22 @@ namespace corsika::pythia8 {
} }
inline bool Interaction::isValid(Code const projectileId, Code const targetId, inline bool Interaction::isValid(Code const projectileId, Code const targetId,
HEPEnergyType const sqrtS) const { HEPEnergyType const sqrtSNN) const {
CORSIKA_LOG_DEBUG("pythia isValid: {}+{} at sqrtS = {} GeV", projectileId, targetId,
sqrtS / 1_GeV); 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 if (is_nucleus(projectileId)) // not yet possible with Pythia
return false; return false;
if (!canInteract(projectileId)) return false; if (!canInteract(projectileId)) return false;
HEPEnergyType const labE = calculate_lab_energy( auto const mass_target =
static_pow<2>(sqrtS), get_mass(projectileId), get_mass(targetId)); 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; if (labE < eKinMinLab_) return false;
bool const validProjectile = bool const validProjectile =
...@@ -220,7 +226,7 @@ namespace corsika::pythia8 { ...@@ -220,7 +226,7 @@ namespace corsika::pythia8 {
CORSIKA_LOG_DEBUG( CORSIKA_LOG_DEBUG(
"Pythia::Interaction: " "Pythia::Interaction: "
"doInteraction: {} interaction? ", "doInteraction: {} interaction? {} ",
projectileId, corsika::pythia8::Interaction::canInteract(projectileId)); projectileId, corsika::pythia8::Interaction::canInteract(projectileId));
// define system (this is Nucleus-Nucleus!!) // define system (this is Nucleus-Nucleus!!)
...@@ -230,10 +236,6 @@ namespace corsika::pythia8 { ...@@ -230,10 +236,6 @@ namespace corsika::pythia8 {
HEPEnergyType const eProjectileLab = calculate_lab_energy( HEPEnergyType const eProjectileLab = calculate_lab_energy(
sqrtS2NN, get_mass(projectileId), get_mass(targetId) / get_nucleus_A(targetId)); 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: ebeam lab: {} GeV", eProjectileLab / 1_GeV);
CORSIKA_LOG_DEBUG( CORSIKA_LOG_DEBUG(
"Interaction: " "Interaction: "
...@@ -241,6 +243,10 @@ namespace corsika::pythia8 { ...@@ -241,6 +243,10 @@ namespace corsika::pythia8 {
" Ecm_NN(GeV): {}", " Ecm_NN(GeV): {}",
eProjectileLab / 1_GeV, sqrtSNN / 1_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)}; COMBoost const labFrameBoost{targetP4.getSpaceLikeComponents(), get_mass(targetId)};
// auto const proj4MomLab = labFrameBoost.toCoM(projectileP4); // auto const proj4MomLab = labFrameBoost.toCoM(projectileP4);
auto const& rotCS = labFrameBoost.getRotatedCS(); auto const& rotCS = labFrameBoost.getRotatedCS();
...@@ -250,8 +256,8 @@ namespace corsika::pythia8 { ...@@ -250,8 +256,8 @@ namespace corsika::pythia8 {
double const eCM_pythia = sqrtSNN / 1_GeV; // CoM energy hadron-Nucleon double const eCM_pythia = sqrtSNN / 1_GeV; // CoM energy hadron-Nucleon
double const idA_pythia = static_cast<double>(get_PDG(projectileId)); double const idA_pythia = static_cast<double>(get_PDG(projectileId));
double const idB_pythia = static_cast<double>(get_PDG(targetId)); double const idB_pythia = static_cast<double>(get_PDG(targetId));
CORSIKA_LOG_INFO("re-configuring pythia idA={}, idB={}, ecm={} GeV", idA_pythia, CORSIKA_LOG_DEBUG("re-configuring pythia idA={}, idB={}, ecm={} GeV", idA_pythia,
idB_pythia, eCM_pythia); idB_pythia, eCM_pythia);
// configure event // configure event
pythia_.setBeamIDs(idA_pythia, idB_pythia); pythia_.setBeamIDs(idA_pythia, idB_pythia);
pythia_.setKinematics(eCM_pythia); pythia_.setKinematics(eCM_pythia);
......
...@@ -122,7 +122,7 @@ namespace corsika::pythia8 { ...@@ -122,7 +122,7 @@ namespace corsika::pythia8 {
bool const print_listing_ = false; bool const print_listing_ = false;
HEPEnergyType const eMaxLab_ = HEPEnergyType const eMaxLab_ =
1e21_eV; // Cross-section tables tabulated up to 10^21 eV 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 } // namespace corsika::pythia8
......
...@@ -81,7 +81,7 @@ SECTION("pythia cross section bug") { ...@@ -81,7 +81,7 @@ SECTION("pythia cross section bug") {
Code const target = Code::Oxygen; Code const target = Code::Oxygen;
Code const projectile = Code::Proton; 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); CORSIKA_LOG_INFO("testing: {} - {}", projectile, target);
REQUIRE( REQUIRE(
...@@ -101,11 +101,10 @@ SECTION("pythia too low energy") { ...@@ -101,11 +101,10 @@ SECTION("pythia too low energy") {
corsika::pythia8::Interaction collision; corsika::pythia8::Interaction collision;
// 5 MeV lab is too low, 0 mb expected // 5 MeV lab is too low, 0 mb expected
REQUIRE_THROWS( REQUIRE(collision.getCrossSectionInelEla(
collision.getCrossSectionInelEla( Code::Proton, Code::Proton,
Code::Proton, Code::Proton, {calculate_total_energy(Proton::mass, 5_MeV), {rootCS, 0_eV, 0_eV, 5_MeV}},
{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});
{Proton::mass, {rootCS, 0_eV, 0_eV, 0_eV}}) == std::tuple{0_mb, 0_mb});
REQUIRE_THROWS(collision.doInteraction( REQUIRE_THROWS(collision.doInteraction(
view, Code::Neutron, Code::Proton, view, Code::Neutron, Code::Proton,
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment