diff --git a/corsika/detail/modules/sibyll/Decay.inl b/corsika/detail/modules/sibyll/Decay.inl
index 465190f6c04a55a99bcfd13bbca4fbeaf052124e..6b8ae281d4a4be07408785e1a846c376c8b4bc0f 100644
--- a/corsika/detail/modules/sibyll/Decay.inl
+++ b/corsika/detail/modules/sibyll/Decay.inl
@@ -173,46 +173,58 @@ namespace corsika::sibyll {
throw std::runtime_error("STOP! Sibyll not configured to execute this decay!");
count_++;
- SibStack ss;
- ss.clear();
- // copy particle to sibyll stack
- ss.addParticle(sibyll::convertToSibyllRaw(pCode), projectile.getEnergy(),
- projectile.getMomentum(),
- // setting particle mass with Corsika values, may be inconsistent
- // with sibyll internal values
- get_mass(pCode));
// remember position
Point const decayPoint = projectile.getPosition();
TimeType const t0 = projectile.getTime();
- // remember if particles is unstable
- // auto const priorIsUnstable = isUnstable(pCode);
// switch on decay for this particle
setUnstable(pCode);
printDecayConfig(pCode);
// call sibyll decay
CORSIKA_LOG_DEBUG("Decay: calling Sibyll decay routine..");
- decsib_();
- if (sibyll_listing_) {
- // print output
- int print_unit = 6;
- sib_list_(print_unit);
- }
+ // particle to pass to sibyll decay
+ int inputSibPID = sibyll::convertToSibyllRaw(pCode);
+
+ // particle momentum format: px, py, pz, e, mass. units: GeV
+ double inputMomentum[5];
+ QuantityVector<hepmomentum_d> input_components =
+ projectile.getMomentum().getComponents();
+ for (int idx = 0; idx < 3; ++idx) inputMomentum[idx] = input_components[idx] / 1_GeV;
+
+ inputMomentum[3] = projectile.getEnergy() / 1_GeV;
+ inputMomentum[4] = get_mass(pCode) / 1_GeV;
+ int nFinalParticles;
+
+ double outputMomentum[5 * 10];
+ int outputSibPID[10];
+
+ // run decay routine
+ decpar_(inputSibPID, inputMomentum, nFinalParticles, outputSibPID,
+ &outputMomentum[0]);
+
+ CORSIKA_LOG_TRACE("Sibyll::Decay: number of final state particles: {}",
+ nFinalParticles);
// reset to stable
setStable(pCode);
- // copy particles from sibyll stack to corsika
- for (auto const& psib : ss) {
- // FOR NOW: skip particles that have decayed in Sibyll, move to iterator?
- if (psib.hasDecayed()) continue;
- // add to corsika stack
- projectile.addSecondary(std::make_tuple(sibyll::convertFromSibyll(psib.getPID()),
- psib.getMomentum(), decayPoint, t0));
+ CoordinateSystemPtr const& rootCS = get_root_CoordinateSystem();
+
+ // copy particles from sibyll ministack to corsika
+ for (int i = 0; i < nFinalParticles; ++i) {
+ QuantityVector<hepmomentum_d> components = {outputMomentum[10 * 0 + i] * 1_GeV,
+ outputMomentum[10 * 1 + i] * 1_GeV,
+ outputMomentum[10 * 2 + i] * 1_GeV};
+
+ auto const pid = sibyll::convertFromSibyll(
+ static_cast<corsika::sibyll::SibyllCode>(outputSibPID[i]));
+
+ CORSIKA_LOG_TRACE("Sibyll::Decay: i={} id={} p={} GeV", i, pid, components / 1_GeV);
+
+ projectile.addSecondary(
+ std::make_tuple(pid, MomentumVector(rootCS, components), decayPoint, t0));
}
- // empty sibyll stack
- ss.clear();
}
} // namespace corsika::sibyll
diff --git a/modules/sibyll/sibyll2.3d.f b/modules/sibyll/sibyll2.3d.f
index 8fccb7a42026eef3621cd80ef20797e27a65a4e1..9d2ff3b7fe07ab7dcaa8145249fa696a38ee2a1a 100644
--- a/modules/sibyll/sibyll2.3d.f
+++ b/modules/sibyll/sibyll2.3d.f
@@ -7,7 +7,7 @@ C SSSSSS IIIIIII BBBBB YY LLLLLLL LLLLLLL
C=======================================================================
C Code for SIBYLL: hadronic interaction Monte Carlo event generator
C=======================================================================
-C Version 2.3d (Jun-01-2017, modified May-20-2020)
+C Version 2.3d01 (Jun-01-2017, modified Oct-15-2021)
C
C with CHARM production
C
@@ -32,7 +32,10 @@ C gaisser@bartol.udel.edu
C paolo.lipari@roma1.infn.it
C friehn@lip.pt
C stanev@bartol.udel.edu
-C
+C
+C last changes relative to Sibyll 2.3d:
+C * stop decays with more than 10 particles
+C
C last changes relative to Sibyll 2.3c:
C * no pi0 suppression in minijets
C * added cross section tables for hadron-nitrogen and hadron-oxygen
@@ -471,8 +474,8 @@ C-----------------------------------------------------------------------
* /,' ','| |',
* /,' ','| Publication to be cited when using this program: |',
* /,' ','| Eun-Joo AHN et al., Phys.Rev. D80 (2009) 094003 |',
- * /,' ','| F. RIEHN et al., hep-ph: 1912.03300 |',
- * /,' ','| last modifications: F. Riehn (05/20/2020) |',
+ * /,' ','| F. RIEHN et al., Phys.Rev.D 102 (2020) 6, 063002 |',
+ * /,' ','| last modifications: F. Riehn (08/15/2021) |',
* /,' ','====================================================',
* /)
@@ -6616,6 +6619,10 @@ C...Choose decay channel
KD =6*(IDC-1)+1
ND = KDEC(KD)
+ IF(ND.GT.10) THEN
+ WRITE(LUN,*) 'DECPAR: too many final state particles in decay!'
+ STOP
+ ENDIF
MAT= KDEC(KD+1)
MBST=0
IF (MAT .GT.0 .AND. P0(4) .GT. 20.D0*P0(5)) MBST=1
diff --git a/modules/sibyll/sibyll2.3d.hpp b/modules/sibyll/sibyll2.3d.hpp
index 741d45775bf00a4503e242cb4245605b8448d89a..b773d0f99789e482c7dba276cddd04be31b440a4 100644
--- a/modules/sibyll/sibyll2.3d.hpp
+++ b/modules/sibyll/sibyll2.3d.hpp
@@ -55,12 +55,12 @@ extern struct {
// number of wounded nucleons, number of hard and soft scatterings etc.
extern struct { int nnsof[20], nnjet[20], jdif[20], nwd, njet, nsof; } s_chist_;
-extern struct {
- double cbr[223 + 16 + 12 + 8];
- int kdec[1338 + 6 * (16 + 12 + 8)];
- int lbarp[99];
- int idb[99];
-} s_csydec_;
+ extern struct {
+ double cbr[223 + 16 + 12 + 8];
+ int kdec[1338 + 6 * (16 + 12 + 8)];
+ int lbarp[99];
+ int idb[99];
+ } s_csydec_;
// additional particle stack for the mother particles of unstable particles
// stable particles have entry zero
@@ -100,17 +100,11 @@ void sibyll_(const int&, const int&, const double&);
// subroutine to initiate sibyll
void sibyll_ini_();
-// subroutine to SET DECAYS
-void dec_ini_();
-
-// subroutine to initiate random number generator
-// void rnd_ini_();
-
// print event
void sib_list_(int&);
-// decay routine
-void decsib_();
+// decay routine (LA,P0,ND,LL,P)
+void decpar_(const int&, const double*, int&, int*, double*);
// interaction length
// double fpni_(double&, int&);
diff --git a/tests/modules/testSibyll.cpp b/tests/modules/testSibyll.cpp
index c5dca766faab79acc981d78c2ef81b7751ef68e6..14051538285f203739b84ae4a3b7eaebad29ada1 100644
--- a/tests/modules/testSibyll.cpp
+++ b/tests/modules/testSibyll.cpp
@@ -100,7 +100,7 @@ auto sumMomentum(TStackView const& view, CoordinateSystemPtr const& vCS) {
return sum;
}
-TEST_CASE("SibyllInterface", "modules") {
+TEST_CASE("SibyllInteractionInterface", "modules") {
logging::set_level(logging::level::info);
@@ -304,6 +304,32 @@ TEST_CASE("SibyllInterface", "modules") {
// CHECK(view.getSize() == 20); // also sibyll not stable wrt. to compiler changes
CHECK(view.getSize() == Approx(100).margin(90)); // this is not physics validation
}
+}
+
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/geometry/RootCoordinateSystem.hpp>
+#include <corsika/framework/geometry/Vector.hpp>
+
+#include <corsika/framework/core/PhysicalUnits.hpp>
+
+#include <corsika/framework/core/ParticleProperties.hpp>
+
+#include <SetupTestEnvironment.hpp>
+#include <SetupTestStack.hpp>
+
+#include <corsika/media/Environment.hpp>
+#include <corsika/media/HomogeneousMedium.hpp>
+#include <corsika/media/NuclearComposition.hpp>
+#include <corsika/media/UniformMagneticField.hpp>
+
+TEST_CASE("SibyllDecayInterface", "modules") {
+ logging::set_level(logging::level::info);
+
+ auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Oxygen);
+ auto const& cs = *csPtr;
+ { [[maybe_unused]] auto const& env_dummy = env; }
+
+ RNGManager<>::getInstance().registerRandomStream("sibyll");
SECTION("DecayInterface") {