diff --git a/applications/c8_air_shower.cpp b/applications/c8_air_shower.cpp
index 8dc0f00ecbb50eed6d749f04951d9cebdba3b987..e215f60228c9ce30b60ef0b9097deb3b37681101 100644
--- a/applications/c8_air_shower.cpp
+++ b/applications/c8_air_shower.cpp
@@ -406,7 +406,8 @@ int main(int argc, char** argv) {
auto const all_elements = corsika::get_all_elements_in_universe(env);
// have SIBYLL always for PROPOSAL photo-hadronic interactions
- auto sibyll = std::make_shared<corsika::sibyll::Interaction>(all_elements);
+ auto sibyll = std::make_shared<corsika::sibyll::Interaction>(
+ all_elements, corsika::setup::C7trackedParticles);
if (auto const modelStr = app["--hadronModel"]->as<std::string>();
modelStr == "SIBYLL-2.3d") {
diff --git a/corsika/detail/modules/sibyll/HadronInteractionModel.inl b/corsika/detail/modules/sibyll/HadronInteractionModel.inl
index 45a5fa957428d7da2c78b3d099a616ddd7144f91..b1b13c676ba12030317f25f5201c0c2f710b5782 100644
--- a/corsika/detail/modules/sibyll/HadronInteractionModel.inl
+++ b/corsika/detail/modules/sibyll/HadronInteractionModel.inl
@@ -26,7 +26,9 @@ namespace corsika::sibyll {
}
inline HadronInteractionModel::HadronInteractionModel()
- : sibyll_listing_(false) {
+ : sibyll_listing_(false)
+ , internal_decays_(false)
+ , stable_particles_{} {
// initialize Sibyll
corsika::connect_random_stream("sibyll", ::sibyll::set_rng_function);
static bool initialized = false;
@@ -36,6 +38,40 @@ namespace corsika::sibyll {
}
}
+ inline HadronInteractionModel::HadronInteractionModel(std::set<Code> vlist)
+ : sibyll_listing_(false)
+ , internal_decays_(true)
+ , stable_particles_(vlist) {
+ // initialize Sibyll
+ corsika::connect_random_stream("sibyll", ::sibyll::set_rng_function);
+ static bool initialized = false;
+ if (!initialized) {
+ sibyll_ini_();
+ initialized = true;
+ }
+ }
+
+ inline void HadronInteractionModel::setAllParticlesUnstable() {
+ // activate all decays in SIBYLL
+ CORSIKA_LOGGER_DEBUG(logger_, "setting all particles \"unstable\".");
+ for (int i = 0; i < 99; ++i) s_csydec_.idb[i] = abs(s_csydec_.idb[i]);
+ }
+
+ inline void HadronInteractionModel::setParticleListStable(std::set<Code> vList) {
+ // de-activate specific decays in SIBYLL
+ for (auto p : vList) {
+ CORSIKA_LOGGER_DEBUG(logger_, "setting {} as \"stable\". ", p);
+ auto const sib_code = sibyll::convertToSibyll(p);
+ if (sib_code != corsika::sibyll::SibyllCode::Unknown) {
+ int const s_id = abs(static_cast<int>(sib_code));
+ s_csydec_.idb[s_id - 1] = (-1) * abs(s_csydec_.idb[s_id - 1]);
+ } else {
+ CORSIKA_LOGGER_WARN(logger_,
+ "particle {} not known to SIBYLL! Cannot set stable!", p);
+ }
+ }
+ }
+
inline HadronInteractionModel::~HadronInteractionModel() {
CORSIKA_LOGGER_DEBUG(logger_, "Sibyll::Model n={}, Nnuc={}", count_, nucCount_);
}
@@ -99,7 +135,7 @@ namespace corsika::sibyll {
sib_sigma_hnuc_(iBeam, iTarget, dEcm, sigProd, dummy, sigEla);
}
return {sigProd * 1_mb, sigEla * 1_mb};
- } // namespace corsika::sibyll
+ }
/**
* In this function SIBYLL is called to produce one event. The
@@ -135,8 +171,13 @@ namespace corsika::sibyll {
count_++;
// Sibyll does not know about units..
double const sqs = sqrtSnn / 1_GeV;
+ if (internal_decays_) {
+ setAllParticlesUnstable();
+ setParticleListStable(stable_particles_);
+ }
// running sibyll, filling stack
sibyll_(projectileSibyllCode, targetSibCode, sqs);
+ if (internal_decays_) decsib_();
if (sibyll_listing_) {
// print final state
@@ -159,9 +200,13 @@ namespace corsika::sibyll {
for (auto& psib : ss) {
// abort on particles that have decayed in Sibyll. Should not happen!
if (psib.hasDecayed()) { // LCOV_EXCL_START
- throw std::runtime_error("found particle that decayed in SIBYLL!");
- } // LCOV_EXCL_STOP
-
+ if (internal_decays_) {
+ continue;
+ } else {
+ throw std::runtime_error(
+ "internal decayse switched off! found particle that decayed in SIBYLL!");
+ } // LCOV_EXCL_STOP
+ }
// transform 4-momentum to lab. frame
// note that the momentum needs to be rotated back
auto const tmp = psib.getMomentum().getComponents();
diff --git a/corsika/detail/modules/sibyll/InteractionModel.inl b/corsika/detail/modules/sibyll/InteractionModel.inl
index 5c2be73da315883e2423347519aac6beb4021938..52bd0c6e7ed5bf7fafa2179ceb6a76c3b5af4137 100644
--- a/corsika/detail/modules/sibyll/InteractionModel.inl
+++ b/corsika/detail/modules/sibyll/InteractionModel.inl
@@ -15,10 +15,11 @@
#include <corsika/modules/sibyll/NuclearInteractionModel.hpp>
namespace corsika::sibyll {
- template <typename TEnvironment>
- inline InteractionModel::InteractionModel(TEnvironment const& environment)
- : hadronSibyll_{}
- , nuclearSibyll_{hadronSibyll_, environment} {}
+
+ inline InteractionModel::InteractionModel(std::set<Code> const& nuccomp,
+ std::set<Code> const& list)
+ : hadronSibyll_{list}
+ , nuclearSibyll_{hadronSibyll_, nuccomp} {}
inline HadronInteractionModel& InteractionModel::getHadronInteractionModel() {
return hadronSibyll_;
diff --git a/corsika/modules/Sibyll.hpp b/corsika/modules/Sibyll.hpp
index b55dc7ec0c3acf43e6622c365dc9cb9074c9c95c..51d0d789db5b2a1de3b117d1406379bb0c999e6b 100644
--- a/corsika/modules/Sibyll.hpp
+++ b/corsika/modules/Sibyll.hpp
@@ -32,8 +32,8 @@ namespace corsika::sibyll {
* to provide all the functions for ProcessSequence.
*/
struct Interaction : public InteractionModel, public InteractionProcess<Interaction> {
- Interaction(std::set<Code> const& nuccomp)
- : InteractionModel{nuccomp} {}
+ Interaction(std::set<Code> const& nuccomp, std::set<Code> const& stablehad)
+ : InteractionModel{nuccomp, stablehad} {}
};
/**
diff --git a/corsika/modules/sibyll/HadronInteractionModel.hpp b/corsika/modules/sibyll/HadronInteractionModel.hpp
index 5eadc97078b6a67088493520ffedb3324045cb45..128d2e500132346a572e020216082dbef6e013ea 100644
--- a/corsika/modules/sibyll/HadronInteractionModel.hpp
+++ b/corsika/modules/sibyll/HadronInteractionModel.hpp
@@ -13,6 +13,7 @@
#include <corsika/framework/random/RNGManager.hpp>
#include <corsika/framework/geometry/FourVector.hpp>
+#include <set>
#include <tuple>
namespace corsika::sibyll {
@@ -27,6 +28,13 @@ namespace corsika::sibyll {
public:
HadronInteractionModel();
+ /**
+ * @brief construct hadron interaction model with SIBYLL
+ *
+ * @param list: list of particles that should be stable inside SIBYLL. all other
+ * particles should decay (if they can)
+ */
+ HadronInteractionModel(std::set<Code> list);
~HadronInteractionModel();
/**
@@ -104,6 +112,8 @@ namespace corsika::sibyll {
FourMomentum const& projectileP4, FourMomentum const& targetP4);
private:
+ void setParticleListStable(std::set<Code>);
+ void setAllParticlesUnstable();
HEPEnergyType constexpr getMinEnergyCoM() const { return minEnergyCoM_; }
HEPEnergyType constexpr getMaxEnergyCoM() const { return maxEnergyCoM_; }
@@ -120,6 +130,8 @@ namespace corsika::sibyll {
int count_ = 0;
int nucCount_ = 0;
bool sibyll_listing_;
+ bool const internal_decays_;
+ std::set<Code> stable_particles_;
};
} // namespace corsika::sibyll
diff --git a/corsika/modules/sibyll/InteractionModel.hpp b/corsika/modules/sibyll/InteractionModel.hpp
index c5fc50087e6e2701a4e1e1e996ff8d489c615583..579f12d6d1b4933234ce83c9d5f932f0e2213404 100644
--- a/corsika/modules/sibyll/InteractionModel.hpp
+++ b/corsika/modules/sibyll/InteractionModel.hpp
@@ -25,8 +25,7 @@ namespace corsika::sibyll {
public:
using nuclear_model_type = NuclearInteractionModel<HadronInteractionModel>;
- template <typename TEnvironment>
- InteractionModel(TEnvironment const&);
+ InteractionModel(std::set<Code> const&, std::set<Code> const&);
CrossSectionType getCrossSection(Code, Code, FourMomentum const&,
FourMomentum const&) const;
diff --git a/examples/cascade_examples/em_shower.cpp b/examples/cascade_examples/em_shower.cpp
index 138717edc83b3b451691a40dd47b25b00899dd78..c3083ec139af35150dff14e1f4943e13cbc64326 100644
--- a/examples/cascade_examples/em_shower.cpp
+++ b/examples/cascade_examples/em_shower.cpp
@@ -42,6 +42,7 @@
#include <corsika/setup/SetupStack.hpp>
#include <corsika/setup/SetupTrajectory.hpp>
+#include <corsika/setup/SetupC7trackedParticles.hpp>
#include <iomanip>
#include <iostream>
@@ -148,7 +149,8 @@ int main(int argc, char** argv) {
ParticleCut<SubWriter<decltype(energyloss)>> cut(5_MeV, 5_MeV, 100_GeV, 100_GeV, true,
energyloss);
- corsika::sibyll::Interaction sibyll(corsika::get_all_elements_in_universe(env));
+ corsika::sibyll::Interaction sibyll(corsika::get_all_elements_in_universe(env),
+ corsika::setup::C7trackedParticles);
corsika::sophia::InteractionModel sophia;
HEPEnergyType heThresholdNN = 80_GeV;
corsika::proposal::Interaction emCascade(
diff --git a/examples/cascade_examples/mars.cpp b/examples/cascade_examples/mars.cpp
index 2610f15d00331d5cc7bb9dcb4091c7612993e8da..b19d1de9e1d36e3caf23d62835139b5384e70ce7 100644
--- a/examples/cascade_examples/mars.cpp
+++ b/examples/cascade_examples/mars.cpp
@@ -60,6 +60,7 @@
#include <corsika/setup/SetupStack.hpp>
#include <corsika/setup/SetupTrajectory.hpp>
+#include <corsika/setup/SetupC7trackedParticles.hpp>
#include <CLI/App.hpp>
#include <CLI/Config.hpp>
@@ -325,7 +326,7 @@ int main(int argc, char** argv) {
/* === START: SETUP PROCESS LIST === */
auto const all_elements = corsika::get_all_elements_in_universe(env);
- corsika::sibyll::Interaction sibyll(all_elements);
+ corsika::sibyll::Interaction sibyll(all_elements, corsika::setup::C7trackedParticles);
InteractionCounter sibyllCounted(sibyll);
corsika::pythia8::Decay decayPythia;
diff --git a/examples/cascade_examples/mc_conex.cpp b/examples/cascade_examples/mc_conex.cpp
index 9fca5112fd528dd3cdfdfc2a080730e6e43955d2..2fbd39ea0e125f84d66a0d2835a21115c6d0f5d0 100644
--- a/examples/cascade_examples/mc_conex.cpp
+++ b/examples/cascade_examples/mc_conex.cpp
@@ -47,6 +47,7 @@
#include <corsika/setup/SetupStack.hpp>
#include <corsika/setup/SetupTrajectory.hpp>
+#include <corsika/setup/SetupC7trackedParticles.hpp>
#include <iomanip>
#include <iostream>
@@ -245,7 +246,8 @@ int main(int argc, char** argv) {
// SETUP PROCESSES, DECAYS, INTERACTIONS
- corsika::sibyll::Interaction sibyll(corsika::get_all_elements_in_universe(env));
+ corsika::sibyll::Interaction sibyll(corsika::get_all_elements_in_universe(env),
+ corsika::setup::C7trackedParticles);
InteractionCounter sibyllCounted{sibyll};
corsika::pythia8::Decay decayPythia;
diff --git a/examples/cascade_examples/radio_em_shower.cpp b/examples/cascade_examples/radio_em_shower.cpp
index f4d5d4524ecdde13088a0d1c9c626ab545052ce7..0b89a1f8adc5b19520ddab62f2c00367b237e7d1 100644
--- a/examples/cascade_examples/radio_em_shower.cpp
+++ b/examples/cascade_examples/radio_em_shower.cpp
@@ -55,6 +55,7 @@
#include <corsika/setup/SetupStack.hpp>
#include <corsika/setup/SetupTrajectory.hpp>
+#include <corsika/setup/SetupC7trackedParticles.hpp>
#include <iomanip>
#include <iostream>
@@ -233,7 +234,8 @@ int main(int argc, char** argv) {
ParticleCut<SubWriter<decltype(energyloss)>> cut(5_MeV, 5_MeV, 100_GeV, 100_GeV, true,
energyloss);
- corsika::sibyll::Interaction sibyll(corsika::get_all_elements_in_universe(env));
+ corsika::sibyll::Interaction sibyll(corsika::get_all_elements_in_universe(env),
+ corsika::setup::C7trackedParticles);
corsika::sophia::InteractionModel sophia;
HEPEnergyType heThresholdNN = 80_GeV;
corsika::proposal::Interaction emCascade(
diff --git a/examples/cascade_examples/water.cpp b/examples/cascade_examples/water.cpp
index 6bf603a8058c259e1b691a826d3772c40f3e62a0..26380f0e72a4b318513a7b6add58226b9c088e97 100644
--- a/examples/cascade_examples/water.cpp
+++ b/examples/cascade_examples/water.cpp
@@ -43,6 +43,7 @@
#include <corsika/setup/SetupStack.hpp>
#include <corsika/setup/SetupTrajectory.hpp>
+#include <corsika/setup/SetupC7trackedParticles.hpp>
#include <CLI/App.hpp>
#include <CLI/Config.hpp>
@@ -223,7 +224,8 @@ int main(int argc, char** argv) {
// hadronic interactions
HEPEnergyType heHadronModelThreshold = std::pow(10, 1.9) * 1_GeV;
- corsika::sibyll::Interaction sibyll(corsika::get_all_elements_in_universe(env));
+ corsika::sibyll::Interaction sibyll(corsika::get_all_elements_in_universe(env),
+ corsika::setup::C7trackedParticles);
auto const all_elements = corsika::get_all_elements_in_universe(env);
corsika::fluka::Interaction leIntModel{all_elements};
diff --git a/modules/sibyll/sibyll2.3d.hpp b/modules/sibyll/sibyll2.3d.hpp
index a00efa88930562009d24bd0ac4214d882fa0bfca..4fa7fe8f4da18e0715d87e2210d0a9295f92fab6 100644
--- a/modules/sibyll/sibyll2.3d.hpp
+++ b/modules/sibyll/sibyll2.3d.hpp
@@ -108,6 +108,8 @@ void sibyll_(const int&, const int&, const double&);
// subroutine to initiate sibyll
void sibyll_ini_();
+void decsib_();
+
// print event
void sib_list_(int&);
diff --git a/tests/modules/testCONEX.cpp b/tests/modules/testCONEX.cpp
index 296700bccab695a83b56be716dfbcc08417d14e0..4ad241839b99529e92915889b1ef094f898745e9 100644
--- a/tests/modules/testCONEX.cpp
+++ b/tests/modules/testCONEX.cpp
@@ -28,6 +28,8 @@
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/utility/CorsikaFenv.hpp>
+#include <corsika/setup/SetupC7trackedParticles.hpp>
+
#include <catch2/catch_all.hpp>
using namespace corsika;
@@ -94,7 +96,7 @@ TEST_CASE("CONEX") {
std::set<Code> const nuclearcomp = {Code::Nitrogen, Code::Oxygen};
// need to initialize Sibyll, done in constructor:
- corsika::sibyll::HadronInteractionModel sibyll;
+ corsika::sibyll::HadronInteractionModel sibyll(corsika::setup::C7trackedParticles);
[[maybe_unused]] corsika::sibyll::NuclearInteractionModel sibyllNuc(sibyll,
nuclearcomp);
diff --git a/tests/modules/testSibyll.cpp b/tests/modules/testSibyll.cpp
index 41a3c88ba3ccf15cf39b5ebae5c3fc440379294b..8e2ecd686e3b7a8b1eec905aa6f2058c95191c42 100644
--- a/tests/modules/testSibyll.cpp
+++ b/tests/modules/testSibyll.cpp
@@ -114,9 +114,11 @@ TEST_CASE("SibyllInterface", "modules") {
std::set<Code> const nuclearcomp = {Code::Hydrogen, Code::Nitrogen, Code::Oxygen,
Code::Carbon};
+
+ std::set<Code> stable_particles = {Code::PiPlus, Code::PiMinus};
SECTION("InteractionInterface - valid targets") {
- corsika::sibyll::HadronInteractionModel model;
+ corsika::sibyll::HadronInteractionModel model(stable_particles);
// sibyll only accepts protons or nuclei with 4<=A<=18 as targets
CHECK_FALSE(model.isValid(Code::Proton, Code::Electron, 100_GeV));
CHECK(model.isValid(Code::Proton, Code::Hydrogen, 100_GeV));
@@ -164,7 +166,7 @@ TEST_CASE("SibyllInterface", "modules") {
const HEPEnergyType P0 = 60_GeV;
MomentumVector const plab = MomentumVector(cs, {P0, 0_eV, 0_eV});
// also print particles after sibyll was called
- corsika::sibyll::HadronInteractionModel model;
+ corsika::sibyll::HadronInteractionModel model(stable_particles);
model.setVerbose(true);
HEPEnergyType const Elab = sqrt(static_pow<2>(P0) + static_pow<2>(Proton::mass));
FourMomentum const projectileP4(Elab, plab);
@@ -248,7 +250,7 @@ TEST_CASE("SibyllInterface", "modules") {
HEPMomentumType const P0 = 50_TeV;
MomentumVector const plab = MomentumVector(cs, {P0, 0_eV, 0_eV});
- corsika::sibyll::HadronInteractionModel hmodel;
+ corsika::sibyll::HadronInteractionModel hmodel(stable_particles);
NuclearInteractionModel nuclearModel(hmodel, nuclearcomp);
@@ -280,7 +282,7 @@ TEST_CASE("SibyllInterface", "modules") {
}
SECTION("CombinedInterface") {
- corsika::sibyll::InteractionModel combinedModel{nuclearcomp};
+ corsika::sibyll::InteractionModel combinedModel{nuclearcomp, stable_particles};
corsika::sibyll::HadronInteractionModel const& hmodel =
combinedModel.getHadronInteractionModel();
auto const& nuclearModel = combinedModel.getNuclearInteractionModel();