From 4935690a83eeb19261085f4148972dd5ccfa5681 Mon Sep 17 00:00:00 2001 From: Felix Riehn <friehn@lip.pt> Date: Thu, 1 Aug 2024 13:29:03 +0200 Subject: [PATCH] check that threshold is high enough for pythia in EAS simulation --- applications/c8_air_shower.cpp | 17 ++++++++++++----- 1 file changed, 12 insertions(+), 5 deletions(-) diff --git a/applications/c8_air_shower.cpp b/applications/c8_air_shower.cpp index ddb9d71ff..ba19ac20f 100644 --- a/applications/c8_air_shower.cpp +++ b/applications/c8_air_shower.cpp @@ -416,6 +416,13 @@ int main(int argc, char** argv) { // have SIBYLL always for PROPOSAL photo-hadronic interactions auto sibyll = std::make_shared<corsika::sibyll::Interaction>(all_elements); + // energy threshold for high energy hadronic model. Affects LE/HE switch for + // hadron interactions and the hadronic photon model in proposal + // valid range varies a bit between models. pythia only has cross sections tabulated + // down to 100GeV. + HEPEnergyType const heHadronModelThreshold = + 1_GeV * app["--hadronModelTransitionEnergy"]->as<double>(); + if (auto const modelStr = app["--hadronModel"]->as<std::string>(); modelStr == "SIBYLL-2.3d") { heModel = DynamicInteractionProcess<StackType>{sibyll}; @@ -428,6 +435,11 @@ int main(int argc, char** argv) { } else if (modelStr == "Pythia8") { heModel = DynamicInteractionProcess<StackType>{ std::make_shared<corsika::pythia8::Interaction>()}; + if (heHadronModelThreshold < 100_GeV) { + CORSIKA_LOG_CRITICAL( + "threshold too low for Pythia! use --hadronModelTransitionEnergy>=100_GeV!"); + return EXIT_FAILURE; + } } else { CORSIKA_LOG_CRITICAL("invalid choice \"{}\"; also check argument parser", modelStr); return EXIT_FAILURE; @@ -471,11 +483,6 @@ int main(int argc, char** argv) { set_energy_production_threshold(Code::TauMinus, std::min({emcut, hadcut, mucut})); set_energy_production_threshold(Code::TauPlus, std::min({emcut, hadcut, mucut})); - // energy threshold for high energy hadronic model. Affects LE/HE switch for - // hadron interactions and the hadronic photon model in proposal - HEPEnergyType const heHadronModelThreshold = - 1_GeV * app["--hadronModelTransitionEnergy"]->as<double>(); - corsika::proposal::Interaction emCascade( env, sophia, sibyll->getHadronInteractionModel(), heHadronModelThreshold); -- GitLab