Merge branch 'sibyll-decay-with-ministack' into 'master'

Overhaul of Sibyll decay

See merge request !392

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

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

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&);

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") {