Merge branch 'sibyll'

CMakeLists.txt

@@ -15,12 +15,14 @@ set(CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/CMakeModules)
 include (CorsikaUtilities) # a few cmake function
 
 # enable warnings and disallow non-standard language
-add_definitions(-Wall -pedantic -Wextra)
+add_definitions(-Wall -pedantic -Wextra -Wno-ignored-qualifiers)
 
 set (CMAKE_CXX_STANDARD 17)
 enable_testing ()
 set (CTEST_OUTPUT_ON_FAILURE 1)
 
+ENABLE_LANGUAGE (Fortran)
+
 # unit testing coverage, does not work yet
 #include (CodeCoverage)
 ##set(COVERAGE_LCOV_EXCLUDES 'Documentation/*')

Documentation/Examples/CMakeLists.txt

@@ -14,6 +14,19 @@ add_executable (stack_example stack_example.cc)
 target_link_libraries (stack_example SuperStupidStack CORSIKAunits CORSIKAlogging)
 install (TARGETS stack_example DESTINATION share/examples)
 
+add_executable (cascade_example cascade_example.cc)
+target_link_libraries (cascade_example SuperStupidStack CORSIKAunits CORSIKAlogging
+   CORSIKArandom
+  CORSIKAsibyll
+  CORSIKAcascade
+  ProcessStackInspector
+  CORSIKAprocesses
+  CORSIKAparticles
+  CORSIKAgeometry
+  CORSIKAprocesssequence
+  )
+install (TARGETS cascade_example DESTINATION share/examples)
+
 add_executable (staticsequence_example staticsequence_example.cc)
 target_link_libraries (staticsequence_example
   CORSIKAprocesssequence

Documentation/Examples/cascade_example.cc

+
+/**
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#include <corsika/cascade/Cascade.h>
+#include <corsika/process/ProcessSequence.h>
+#include <corsika/process/stack_inspector/StackInspector.h>
+
+#include <corsika/setup/SetupStack.h>
+#include <corsika/setup/SetupTrajectory.h>
+
+#include <corsika/random/RNGManager.h>
+
+#include <corsika/cascade/sibyll2.3c.h>
+#include <corsika/cascade/SibStack.h>
+#include <corsika/process/sibyll/ParticleConversion.h>
+
+#include <corsika/units/PhysicalUnits.h>
+using namespace corsika;
+using namespace corsika::process;
+using namespace corsika::units;
+using namespace corsika::particles;
+using namespace corsika::random;
+
+#include <typeinfo>
+#include <iostream>
+using namespace std;
+
+static int fCount = 0;
+
+class ProcessSplit : public corsika::process::BaseProcess<ProcessSplit> {
+public:
+  ProcessSplit() {}
+
+  template <typename Particle>
+  double MinStepLength(Particle& p) const {
+
+    // beam particles for sibyll : 1, 2, 3 for p, pi, k
+    // read from cross section code table
+    int kBeam   =  1;
+    
+    /* 
+       the target should be defined by the Environment,
+       ideally as full particle object so that the four momenta 
+       and the boosts can be defined..
+     */
+    // target nuclei: A < 18
+    // FOR NOW: assume target is oxygen
+    int kTarget = 1;
+    double beamEnergy =  p.GetEnergy() / 1_GeV; 
+    std::cout << "ProcessSplit: " << "MinStep: en: " << beamEnergy << " pid:" << kBeam << std::endl;
+    double prodCrossSection,dummy,dum1,dum2,dum3,dum4;
+    double dumdif[3];
+
+    if(kTarget==1)
+      sib_sigma_hp_(kBeam, beamEnergy, dum1, dum2, prodCrossSection, dumdif,dum3, dum4 );
+    else
+      sib_sigma_hnuc_(kBeam, kTarget, beamEnergy, prodCrossSection, dummy );
+    
+    std::cout << "ProcessSplit: " << "MinStep: sibyll return: " << prodCrossSection << std::endl;
+    CrossSectionType sig = prodCrossSection * 1_mbarn;
+    std::cout << "ProcessSplit: " << "MinStep: CrossSection (mb): " << sig / 1_mbarn << std::endl;
+
+    const MassType nucleon_mass = 0.93827_GeV / corsika::units::si::constants::cSquared;
+    std::cout << "ProcessSplit: " << "nucleon mass " << nucleon_mass << std::endl;
+    // calculate interaction length in medium
+    double int_length = kTarget * ( nucleon_mass / 1_g ) / ( sig / 1_cmeter / 1_cmeter );
+    // pick random step lenth
+    std::cout << "ProcessSplit: " << "interaction length (g/cm2): " << int_length << std::endl;
+    // add exponential sampling
+    int a = 0;
+    const double next_step = -int_length * log(s_rndm_(a));
+    /*
+      what are the units of the output? slant depth or 3space length?
+      
+    */
+    std::cout << "ProcessSplit: " << "next step (g/cm2): " << next_step << std::endl;
+    return next_step;
+  }
+
+  template <typename Particle, typename Trajectory, typename Stack>
+  EProcessReturn DoContinuous(Particle&, Trajectory&, Stack&) const {
+    // corsika::utls::ignore(p);
+    return EProcessReturn::eOk;
+  }
+
+  template <typename Particle, typename Stack>
+  void DoDiscrete(Particle& p, Stack& s) const {
+    cout << "DoDiscrete: " << p.GetPID() << " interaction? " << process::sibyll::CanInteract( p.GetPID() )  << endl; 
+    if( process::sibyll::CanInteract( p.GetPID() ) ){
+    
+    // get energy of particle from stack
+    /*
+      stack is in GeV in lab. frame
+      convert to GeV in cm. frame 
+      (assuming proton at rest as target AND 
+      assuming no pT, i.e. shower frame-z is aligned with hadron-int-frame-z)
+    */
+    EnergyType E   = p.GetEnergy();
+    EnergyType Ecm = sqrt( 2. * E * 0.93827_GeV );
+
+    int kBeam   = process::sibyll::ConvertToSibyllRaw( p.GetPID() );
+    
+    /* 
+       the target should be defined by the Environment,
+       ideally as full particle object so that the four momenta 
+       and the boosts can be defined..
+     */
+    // FOR NOW: set target to proton
+    int kTarget = 1; //p.GetPID();
+  
+    std::cout << "ProcessSplit: " << " DoDiscrete: E(GeV):" << E / 1_GeV << " Ecm(GeV): " << Ecm / 1_GeV << std::endl;
+    if (E < 8.5_GeV || Ecm < 10_GeV ) {
+      std::cout << "ProcessSplit: " << " DoDiscrete: dropping particle.." << std::endl;
+      p.Delete();
+      fCount++;
+    } else {
+      // Sibyll does not know about units..
+      double sqs = Ecm / 1_GeV;
+      // running sibyll, filling stack
+      sibyll_( kBeam, kTarget, sqs);
+      // running decays
+      //decsib_();
+      // print final state
+      int print_unit = 6;
+      sib_list_( print_unit );
+
+      // delete current particle
+      p.Delete();
+
+      // add particles from sibyll to stack
+      // link to sibyll stack
+      SibStack ss;
+      /*
+	get transformation between Stack-frame and SibStack-frame
+	for EAS Stack-frame is lab. frame, could be different for CRMC-mode
+	the transformation should be derived from the input momenta
+	in general transformation is rotation + boost
+      */
+      const EnergyType proton_mass = 0.93827_GeV;
+      const double gamma  = ( E + proton_mass ) / ( Ecm );
+      const double gambet =  sqrt( E * E - proton_mass * proton_mass ) / Ecm;
+
+      // SibStack does not know about momentum yet so we need counter to access momentum array in Sibyll
+      int i = -1;
+      for (auto &p: ss){
+	++i;
+	//transform to lab. frame, primitve
+	const double en_lab = gambet * s_plist_.p[2][i] + gamma * p.GetEnergy();	
+	// add to corsika stack
+	auto pnew = s.NewParticle();
+	pnew.SetEnergy( en_lab * 1_GeV );
+	pnew.SetPID( process::sibyll::ConvertFromSibyll( p.GetPID() ) );
+      }
+    }
+    }else
+      p.Delete();
+  }
+  
+  
+  void Init() 
+  {
+    fCount = 0;
+
+    // define reference frame? --> defines boosts between corsika stack and model stack.
+    
+    // initialize random numbers for sibyll
+    // FOR NOW USE SIBYLL INTERNAL !!!
+    //    rnd_ini_();
+    
+    corsika::random::RNGManager & rmng = corsika::random::RNGManager::GetInstance();;
+    const std::string str_name = "s_rndm";
+    rmng.RegisterRandomStream(str_name);
+    
+    // //    corsika::random::RNG srng;
+    // auto srng = rmng.GetRandomStream("s_rndm");
+
+    // test random number generator
+    std::cout << "ProcessSplit: " << " test sequence of random numbers."  << std::endl;
+    int a = 0;
+    for(int i=0; i<8; ++i)
+      std::cout << i << " " << s_rndm_(a) << std::endl;
+    
+    //initialize Sibyll
+    sibyll_ini_();
+
+    // set particles stable / unstable
+    // use stack to loop over particles
+    setup::Stack ds;
+    ds.NewParticle().SetPID( Code::Proton  );
+    ds.NewParticle().SetPID( Code::Neutron );
+    ds.NewParticle().SetPID( Code::PiPlus  );
+    ds.NewParticle().SetPID( Code::PiMinus );
+    ds.NewParticle().SetPID( Code::KPlus   );
+    ds.NewParticle().SetPID( Code::KMinus  );
+    ds.NewParticle().SetPID( Code::K0Long  );
+    ds.NewParticle().SetPID( Code::K0Short );
+
+    for( auto &p: ds){
+      int s_id = process::sibyll::ConvertToSibyllRaw( p.GetPID() );
+      // set particle stable by setting table value negative
+      cout << "ProcessSplit: Init: setting " << p.GetPID() << "(" << s_id << ")" << " stable in Sibyll .." << endl;
+      s_csydec_.idb[ s_id ] = -s_csydec_.idb[ s_id-1 ];
+      p.Delete();
+    }
+  }
+  
+  int GetCount() { return fCount; }
+
+private:
+};
+
+double s_rndm_(int &)
+{
+  static corsika::random::RNG& rmng = corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");;
+  return rmng()/(double)rmng.max();
+}
+
+
+int main(){
+
+  stack_inspector::StackInspector<setup::Stack, setup::Trajectory> p0(true);
+  ProcessSplit p1;
+  const auto sequence = p0 + p1;
+  setup::Stack stack;
+
+  
+  corsika::cascade::Cascade EAS(sequence, stack);
+
+  stack.Clear();
+  auto particle = stack.NewParticle();
+  EnergyType E0 = 100_GeV;
+  particle.SetEnergy(E0);
+  particle.SetPID( Code::Proton );
+  EAS.Init();
+  EAS.Run();
+  cout << "Result: E0=" << E0 / 1_GeV << "GeV, count=" << p1.GetCount() << endl;
+  
+}

Framework/Cascade/CMakeLists.txt

@@ -9,6 +9,8 @@ set (
 set (
   CORSIKAcascade_HEADERS
   Cascade.h
+  sibyll2.3c.h
+  SibStack.h
   )
 
 #set (
@@ -16,9 +18,23 @@ set (
 #  Cascade.cc
 #  )
 
+set (
+  CORSIKAsibyll_NAMESPACE
+  corsika/cascade
+  )
+
+set (
+  CORSIKAsibyll_SOURCES
+  sibyll2.3c.f
+  gasdev.f
+  )
+add_library (CORSIKAsibyll STATIC ${CORSIKAsibyll_SOURCES})
+
 #add_library (CORSIKAcascade STATIC ${CORSIKAcascade_SOURCES})
 add_library (CORSIKAcascade INTERFACE)
 
+
+
 CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKAcascade ${CORSIKAcascade_NAMESPACE} ${CORSIKAcascade_HEADERS})
 
 #target_link_libraries (
@@ -51,9 +67,12 @@ add_executable (
 
 target_link_libraries (
   testCascade
-#  CORSIKAutls
+  #  CORSIKAutls
+  CORSIKArandom
+  CORSIKAsibyll
   CORSIKAcascade
   ProcessStackInspector
+  CORSIKAstackinterface
   CORSIKAprocesses
   CORSIKAparticles
   CORSIKAgeometry

Framework/Cascade/Cascade.cc

@@ -38,5 +38,6 @@ void Cascade::Step(Particle& particle) {
   corsika::geometry::LineTrajectory trajectory =
       fProcesseList.Transport(particle, nextStep);
   sequence.DoContinuous(particle, trajectory);
+  // whats going on here? Everywhere else DoDiscrete is passed a Stack reference as well
   sequence.DoDiscrete(particle);
 }

Framework/Cascade/SibStack.h

+#ifndef _include_sibstack_h_
+#define _include_sibstack_h_
+
+#include <vector>
+#include <string>
+
+#include <corsika/stack/Stack.h>
+#include <corsika/cascade/sibyll2.3c.h>
+#include <corsika/process/sibyll/ParticleConversion.h>
+
+using namespace std;
+using namespace corsika::stack;
+
+class SibStackData {
+  
+ public:
+  void Init();
+  
+  void Clear() { s_plist_.np = 0; }
+  
+  int GetSize() const { return s_plist_.np;  }
+  int GetCapacity() const { return 8000; }
+
+  
+  void SetId(const int i, const int v) { s_plist_.llist[i]=v; }
+  void SetEnergy(const int i, const double v) { s_plist_.p[3][i]=v;  }
+
+  int GetId(const int i) const { return s_plist_.llist[i]; }
+  double GetEnergy(const int i) const { return s_plist_.p[3][i]; }
+  
+  void Copy(const int i1, const int i2) {
+    s_plist_.llist[i1] = s_plist_.llist[i2];
+    s_plist_.p[3][i1] = s_plist_.p[3][i2];
+  }
+  
+ protected:
+  void IncrementSize() { s_plist_.np++; }
+  void DecrementSize() { if ( s_plist_.np>0) { s_plist_.np--; } }
+};
+
+
+template<typename StackIteratorInterface>
+class ParticleInterface : public ParticleBase<StackIteratorInterface> {
+  using ParticleBase<StackIteratorInterface>::GetStackData;
+  using ParticleBase<StackIteratorInterface>::GetIndex;
+ public:
+  void SetEnergy(const double v) { GetStackData().SetEnergy(GetIndex(), v); }
+  double GetEnergy() const { return GetStackData().GetEnergy(GetIndex()); }
+  void SetPID(const int v) { GetStackData().SetId(GetIndex(), v); }
+  corsika::process::sibyll::SibyllCode GetPID() const { return static_cast<corsika::process::sibyll::SibyllCode> (GetStackData().GetId(GetIndex())); }
+  
+};
+
+
+typedef Stack<SibStackData, ParticleInterface> SibStack;
+
+#endif

Framework/Cascade/gasdev.f

+      DOUBLE PRECISION FUNCTION GASDEV(Idum)
+C***********************************************************************
+C     Gaussian deviation
+C***********************************************************************
+      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+      IMPLICIT INTEGER(I-N)
+      COMMON /RNDMGAS/ ISET
+      SAVE
+      DATA ISET/0/      
+      gasdev=idum
+      IF (ISET.EQ.0) THEN
+1       V1=2.D0*S_RNDM(0)-1.D0
+        V2=2.D0*S_RNDM(1)-1.D0
+        R=V1**2+V2**2
+        IF(R.GE.1.D0)GO TO 1
+        FAC=SQRT(-2.D0*LOG(R)/R)
+        GSET=V1*FAC
+        GASDEV=V2*FAC
+        ISET=1
+      ELSE
+        GASDEV=GSET
+        ISET=0
+      ENDIF
+      RETURN
+      END
+C***********************************************************************
+ 

Framework/Cascade/rndm_dbl.f

+C***********************************************************************
+C
+C    interface to PHOJET double precision random number generator 
+C    for SIBYLL \FR'14
+C
+C***********************************************************************
+      DOUBLE PRECISION FUNCTION S_RNDM(IDUMMY)
+      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+      DUMMY = dble(IDUMMY)
+      S_RNDM= PHO_RNDM(DUMMY)
+      END
+
+C***********************************************************************
+C
+C    initialization routine for double precision random number generator
+C    calls PHO_RNDIN \FR'14
+C
+C***********************************************************************
+      SUBROUTINE RND_INI
+      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+      COMMON /RNDMGAS/ ISET
+      ISET = 0
+      CALL PHO_RNDIN(12,34,56,78)
+      END
+
+
+      DOUBLE PRECISION FUNCTION GASDEV(Idum)
+C***********************************************************************
+C     Gaussian deviation
+C***********************************************************************
+      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+      IMPLICIT INTEGER(I-N)
+      COMMON /RNDMGAS/ ISET
+      SAVE
+      DATA ISET/0/      
+      gasdev=idum
+      IF (ISET.EQ.0) THEN
+1       V1=2.D0*S_RNDM(0)-1.D0
+        V2=2.D0*S_RNDM(1)-1.D0
+        R=V1**2+V2**2
+        IF(R.GE.1.D0)GO TO 1
+        FAC=SQRT(-2.D0*LOG(R)/R)
+        GSET=V1*FAC
+        GASDEV=V2*FAC
+        ISET=1
+      ELSE
+        GASDEV=GSET
+        ISET=0
+      ENDIF
+      RETURN
+      END
+C***********************************************************************
+      
+
+      DOUBLE PRECISION FUNCTION PHO_RNDM(DUMMY)
+C***********************************************************************
+C
+C    random number generator
+C
+C    initialization by call to PHO_RNDIN needed!
+C     
+C    the algorithm is taken from
+C      G.Marsaglia, A.Zaman: 'Toward a unversal random number generator'
+C      Florida State Univ. preprint FSU-SCRI-87-70
+C
+C    implementation by K. Hahn (Dec. 88), changed to include possibility
+C    of saving / reading generator registers to / from file (R.E. 10/98)
+C
+C    generator should not depend on the hardware (if a real has
+C    at least 24 significant bits in internal representation),
+C    the period is about 2**144,
+C
+C    internal registers:
+C       U(97),C,CD,CM,I,J  - seed values as initialized in PHO_RNDIN
+C
+C
+C***********************************************************************
+      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+      SAVE
+
+      COMMON /PORAND/ U(97),C,CD,CM,I,J
+
+ 100  CONTINUE
+      RNDMI = DUMMY
+      RNDMI = U(I)-U(J)
+      IF ( RNDMI.LT.0.D0 ) RNDMI = RNDMI+1.D0
+      U(I) = RNDMI
+      I    = I-1
+      IF ( I.EQ.0 ) I = 97
+      J    = J-1
+      IF ( J.EQ.0 ) J = 97
+      C    = C-CD
+      IF ( C.LT.0.D0 ) C = C+CM
+      RNDMI = RNDMI-C
+      IF ( RNDMI.LT.0.D0 ) RNDMI = RNDMI+1.D0
+
+      IF((ABS(RNDMI).LT.0.D0).OR.(ABS(RNDMI-1.D0).LT.1.D-10)) GOTO 100
+      PHO_RNDM = RNDMI
+
+      END
+
+
+CDECK  ID>, PHO_RNDIN
+      SUBROUTINE PHO_RNDIN(NA1,NA2,NA3,NB1)
+C***********************************************************************
+C
+C     initialization of PHO_RNDM, has to be called before using PHO_RNDM
+C
+C     input:
+C       NA1,NA2,NA3,NB1  - values for initializing the generator
+C                          NA? must be in 1..178 and not all 1;
+C                          12,34,56  are the standard values
+C                          NB1 must be in 1..168;
+C                          78  is the standard value
+C
+C***********************************************************************
+      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+      SAVE
+
+      COMMON /PORAND/ U(97),C,CD,CM,I,J
+      MA1 = NA1
+      MA2 = NA2
+      MA3 = NA3
+      MB1 = NB1
+      I   = 97
+      J   = 33
+      DO 20 II2 = 1,97
+        S = 0.D0
+        T = 0.5D0
+        DO 10 II1 = 1,24
+          MAT  = MOD(MOD(MA1*MA2,179)*MA3,179)
+          MA1  = MA2
+          MA2  = MA3
+          MA3  = MAT
+          MB1  = MOD(53*MB1+1,169)
+          IF ( MOD(MB1*MAT,64).GE.32 ) S = S+T
+          T    = 0.5D0*T
+ 10     CONTINUE
+        U(II2) = S
+ 20   CONTINUE
+      C  =   362436.D0/16777216.D0
+      CD =  7654321.D0/16777216.D0
+      CM = 16777213.D0/16777216.D0
+
+      END
+
+
+CDECK  ID>, PHO_RNDSI
+      SUBROUTINE PHO_RNDSI(UIN,CIN,CDIN,CMIN,IIN,JIN)
+C***********************************************************************
+C
+C     updates internal random number generator registers using
+C     registers given as arguments
+C
+C***********************************************************************
+      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+      SAVE
+
+      DIMENSION UIN(97)
+      COMMON /PORAND/ U(97),C,CD,CM,I,J
+      DO 10 KKK = 1,97
+        U(KKK) = UIN(KKK)
+ 10   CONTINUE
+      C  = CIN
+      CD = CDIN
+      CM = CMIN
+      I  = IIN
+      J  = JIN
+
+      END
+
+
+CDECK  ID>, PHO_RNDSO
+      SUBROUTINE PHO_RNDSO(UOUT,COUT,CDOUT,CMOUT,IOUT,JOUT)
+C***********************************************************************
+C
+C     copies internal registers from randon number generator
+C     to arguments
+C
+C***********************************************************************
+      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+      SAVE
+
+      DIMENSION UOUT(97)
+      COMMON /PORAND/ U(97),C,CD,CM,I,J
+      DO 10 KKK = 1,97
+        UOUT(KKK) = U(KKK)
+ 10   CONTINUE
+      COUT  = C
+      CDOUT = CD
+      CMOUT = CM
+      IOUT  = I
+      JOUT  = J
+
+      END
+
+
+CDECK  ID>, PHO_RNDTE
+      SUBROUTINE PHO_RNDTE(IO)
+C***********************************************************************
+C
+C     test of random number generator PHO_RNDM
+C
+C     input:
+C       IO defines output
+C           0  output only if an error is detected
+C           1  output independend on an error
+C
+C     uses PHO_RNDSI and PHO_RNDSO to bring the random number generator
+C     to same status as it had before the test run
+C
+C***********************************************************************
+      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+      SAVE
+
+
+C  input/output channels
+      INTEGER LI,LO
+      COMMON /POINOU/ LI,LO
+
+      DIMENSION UU(97)
+      DIMENSION U(6),X(6),D(6)
+      DATA U / 6533892.D0 , 14220222.D0 ,  7275067.D0 ,
+     &         6172232.D0 ,  8354498.D0 , 10633180.D0 /
+
+      CALL PHO_RNDSO(UU,CC,CCD,CCM,II,JJ)
+
+      CALL PHO_RNDIN(12,34,56,78)
+      DO 10 II1 = 1,20000
+        XX      = PHO_RNDM(SD)
+ 10   CONTINUE
+
+      SD        = 0.D0
+      DO 20 II2 = 1,6
+        X(II2)  = 4096.D0*(4096.D0*PHO_RNDM(XX))
+        D(II2)  = X(II2)-U(II2)
+        SD      = SD+ABS(D(II2))
+ 20   CONTINUE
+
+      CALL PHO_RNDSI(UU,CC,CCD,CCM,II,JJ)
+
+      IF ((IO.EQ.1).OR.(ABS(SD).GT.0.D-10)) THEN
+        WRITE(LO,50) (U(I),X(I),D(I),I=1,6)
+      ENDIF
+
+ 50   FORMAT(/,' PHO_RNDTE: test of the random number generator:',/,
+     &  '    expected value    calculated value     difference',/,
+     &  6(F17.1,F20.1,F15.3,/),
+     &  ' generator has the same status as before calling PHO_RNDTE',/)
+
+      END
+
+
+CDECK  ID>, PHO_RNDST
+      SUBROUTINE PHO_RNDST(MODE,FILENA)
+C***********************************************************************
+C
+C     read / write random number generator status from / to file
+C
+C     input:    MODE        1   read registers from file
+C                           2   dump registers to file
+C
+C               FILENA      file name
+C
+C***********************************************************************
+
+      IMPLICIT NONE
+
+
+
+      SAVE
+
+      INTEGER       MODE
+      CHARACTER*(*) FILENA
+
+
+C  input/output channels
+      INTEGER LI,LO
+      COMMON /POINOU/ LI,LO
+
+
+      DOUBLE PRECISION UU,CC,CCD,CCM
+      DIMENSION UU(97)
+
+      INTEGER I,II,JJ
+
+      CHARACTER*80 CH_DUMMY
+
+      IF(MODE.EQ.1) THEN
+
+        WRITE(LO,'(/,1X,2A,A,/)') 'PHO_RNDST: ',
+     &    'reading random number registers from file ',FILENA
+
+        OPEN(12,FILE=FILENA,ERR=1010,STATUS='OLD')
+        READ(12,*,ERR=1010) CH_DUMMY
+        DO I=1,97
+          READ(12,*,ERR=1010) UU(I)
+        ENDDO
+        READ(12,*,ERR=1010) CC
+        READ(12,*,ERR=1010) CCD
+        READ(12,*,ERR=1010) CCM
+        READ(12,*,ERR=1010) II,JJ
+        CLOSE(12)
+        CALL PHO_RNDSI(UU,CC,CCD,CCM,II,JJ)
+
+      ELSE IF(MODE.EQ.2) THEN
+
+        WRITE(LO,'(/,1X,2A,A,/)') 'PHO_RNDST: ',
+     &    'dumping random number registers to file ',FILENA
+
+        OPEN(12,FILE=FILENA,ERR=1010,STATUS='UNKNOWN')
+        CALL PHO_RNDSO(UU,CC,CCD,CCM,II,JJ)
+        WRITE(12,'(1X,A)',ERR=1020) 'random number status registers:'
+        DO I=1,97
+          WRITE(12,'(1X,1P,E28.20)',ERR=1020) UU(I)
+        ENDDO
+        WRITE(12,'(1X,1P,E28.20)',ERR=1020) CC
+        WRITE(12,'(1X,1P,E28.20)',ERR=1020) CCD
+        WRITE(12,'(1X,1P,E28.20)',ERR=1020) CCM
+        WRITE(12,'(1X,2I4)',ERR=1020) II,JJ
+        CLOSE(12)
+
+      ELSE
+
+        WRITE(LO,'(/,1X,2A,I6,/)') 'PHO_RNDST: ',
+     &    'called with invalid mode, nothing done (mode)',MODE
+
+      ENDIF
+
+      RETURN
+
+ 1010 CONTINUE
+      WRITE(LO,'(1X,2A,A,/)') 'PHO_RNDST: ',
+     &  'cannot open or read file ',FILENA
+      RETURN
+
+ 1020 CONTINUE
+      WRITE(LO,'(1X,A,A,/)') 'PHO_RNDST: ',
+     &  'cannot open or write file ',FILENA
+      RETURN
+    
+      END
+
+C----------------------------------------
+C standard generator
+C----------------------------------------
+      REAL FUNCTION S_RNDM_std(IDUMMY)
+C...Generator  from the LUND montecarlo
+C...Purpose: to generate random numbers uniformly distributed between
+C...0 and 1, excluding the endpoints.
+      COMMON/LUDATR/MRLU(6),RRLU(100)
+      SAVE /LUDATR/
+      EQUIVALENCE (MRLU1,MRLU(1)),(MRLU2,MRLU(2)),(MRLU3,MRLU(3)),
+     &(MRLU4,MRLU(4)),(MRLU5,MRLU(5)),(MRLU6,MRLU(6)),
+     &(RRLU98,RRLU(98)),(RRLU99,RRLU(99)),(RRLU00,RRLU(100))
+ 
+C...  Initialize generation from given seed.
+      S_RNDM_std = real(idummy)
+      IF(MRLU2.EQ.0) THEN
+        IF (MRLU1 .EQ. 0)  MRLU1 = 19780503    ! initial seed
+        IJ=MOD(MRLU1/30082,31329)
+        KL=MOD(MRLU1,30082)
+        I=MOD(IJ/177,177)+2
+        J=MOD(IJ,177)+2
+        K=MOD(KL/169,178)+1
+        L=MOD(KL,169)
+        DO 110 II=1,97
+        S=0.
+        T=0.5
+        DO 100 JJ=1,24
+        M=MOD(MOD(I*J,179)*K,179)
+        I=J
+        J=K
+        K=M
+        L=MOD(53*L+1,169)
+        IF(MOD(L*M,64).GE.32) S=S+T
+        T=0.5*T
+  100   CONTINUE
+        RRLU(II)=S
+  110   CONTINUE
+        TWOM24=1.
+        DO 120 I24=1,24
+        TWOM24=0.5*TWOM24
+  120   CONTINUE
+        RRLU98=362436.*TWOM24
+        RRLU99=7654321.*TWOM24
+        RRLU00=16777213.*TWOM24
+        MRLU2=1
+        MRLU3=0
+        MRLU4=97
+        MRLU5=33
+      ENDIF
+ 
+C...Generate next random number.
+  130 RUNI=RRLU(MRLU4)-RRLU(MRLU5)
+      IF(RUNI.LT.0.) RUNI=RUNI+1.
+      RRLU(MRLU4)=RUNI
+      MRLU4=MRLU4-1
+      IF(MRLU4.EQ.0) MRLU4=97
+      MRLU5=MRLU5-1
+      IF(MRLU5.EQ.0) MRLU5=97
+      RRLU98=RRLU98-RRLU99
+      IF(RRLU98.LT.0.) RRLU98=RRLU98+RRLU00
+      RUNI=RUNI-RRLU98
+      IF(RUNI.LT.0.) RUNI=RUNI+1.
+      IF(RUNI.LE.0.OR.RUNI.GE.1.) GOTO 130
+ 
+C...Update counters. Random number to output.
+      MRLU3=MRLU3+1
+      IF(MRLU3.EQ.1000000000) THEN
+        MRLU2=MRLU2+1
+        MRLU3=0
+      ENDIF
+      S_RNDM_std=RUNI
+
+      END

Framework/Cascade/sibyll2.3c.h

+#ifndef _include_sib23c_interface_h_
+#define _include_sib23c_interface_h_
+//----------------------------------------------
+//  C++ interface for the SIBYLL event generator
+//----------------------------------------------
+//wrapper
+
+
+
+
+extern"C"{
+
+  typedef char s_name[6];
+
+  // SIBYLL particle stack (FORTRAN COMMON)
+  // variables are: np : numer of particles on stack
+  //                 p : 4momentum + mass of particles on stack
+  //             llist : id of particles on stack
+  extern struct { 
+    double p[5][8000]; 
+    int llist[8000]; 
+    int np; 
+  }s_plist_;
+
+
+  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
+  extern struct { 
+    int llist1[8000]; 
+  }s_plist1_;
+
+
+  // tables with particle properties
+  // charge, strangeness and baryon number
+  extern struct {
+    int ichp[99]; 
+    int istr[99]; 
+    int ibar[99]; 
+  }s_chp_;
+
+  // tables with particle properties
+  // mass and mass squared
+  extern struct {
+    double am[99]; 
+    double am2[99]; 
+  }s_mass1_;
+
+  // table with particle names
+  extern struct {char namp[6][99];}s_cnam_;
+
+  // debug info
+  extern struct {int ncall; int ndebug; int lun; }s_debug_;
+
+  // lund random generator setup
+  //extern struct {int mrlu[6]; float rrlu[100]; }ludatr_;
+
+
+  // sibyll main subroutine
+  void sibyll_(int&,int&,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_();
+
+  // interaction length
+  //double fpni_(double&, int&);
+
+  void sib_sigma_hnuc_(int&,int&,double&,double&,double&);
+  void sib_sigma_hp_(int&,double&,double&,double&,double&,double*,double&,double&);
+  
+  double s_rndm_(int&);
+
+  // phojet random generator setup
+  void pho_rndin_(int&,int&,int&,int&);
+}
+#endif
+

Framework/Units/PhysicalConstants.h

@@ -54,6 +54,10 @@ namespace corsika::units::si::constants {
   // unified atomic mass unit
   constexpr quantity<mass_d> u{Rep(1.6605402e-27L) * kilogram};
 
+  // barn moved to PhysicalUnits
+  //  constexpr quantity<area_d> barn{Rep(1.e-28L) * meter * meter};
+  
+  
   // etc.
 
 } // namespace corsika::units::si::constants

Framework/Units/PhysicalUnits.h

@@ -9,26 +9,25 @@
 /**
  * @file PhysicalUnits
  *
- * Define _XeV literals, alowing 10_GeV in the code.
+ * Define new units and unit-types
  */
 
-namespace phys {
-  namespace units {
-    namespace literals {
-      QUANTITY_DEFINE_SCALING_LITERALS(eV, energy_d,
-                                       magnitude(corsika::units::si::constants::eV))
-
-      // phys::units::quantity<energy_d/mass_d> Joule2Kg = c2; // 1_Joule / 1_kg;
-
-    } // namespace literals
-  }   // namespace units
-} // namespace phys
-
 namespace corsika::units::si {
   using namespace phys::units;
   using namespace phys::units::literals;
   // namespace literals = phys::units::literals;
+  
+  /// defining momentum you suckers
+  /// dimensions, i.e. composition in base SI dimensions
+  using momentum_d = phys::units::dimensions<1, 1, -1>;
+  // defining the unit of momentum, so far newton-meter, maybe go to HEP?
+  constexpr phys::units::quantity<momentum_d> newton_second{meter * kilogram / second};
 
+  /// defining cross section
+  using sigma_d = phys::units::dimensions<2, 0, 0>;
+  constexpr phys::units::quantity<sigma_d> barn{Rep(1.e-28L) * meter * meter};
+
+  /// add the unit-types
   using LengthType = phys::units::quantity<phys::units::length_d, double>;
   using TimeType = phys::units::quantity<phys::units::time_interval_d, double>;
   using SpeedType = phys::units::quantity<phys::units::speed_d, double>;
@@ -36,10 +35,44 @@ namespace corsika::units::si {
   using ElectricChargeType =
       phys::units::quantity<phys::units::electric_charge_d, double>;
   using EnergyType = phys::units::quantity<phys::units::energy_d, double>;
-  using MassType = phys::units::quantity<phys::units::mass_d, double>;
+  using MassType = phys::units::quantity<phys::units::mass_d, double>;  
+  using MomentumType = phys::units::quantity<momentum_d, double>;
+  using CrossSectionType = phys::units::quantity<sigma_d, double>;
 
 } // end namespace corsika::units::si
 
+
+
+/**
+ * @file PhysicalUnits
+ *
+ * Define _XeV literals, alowing 10_GeV in the code.
+ * Define _meter literal
+ * Define _barn literal
+ * Define _newton_second literal for SI momenta
+ */
+
+namespace phys {
+  namespace units {
+    namespace literals {
+      QUANTITY_DEFINE_SCALING_LITERALS(eV, energy_d,
+                                       magnitude(corsika::units::si::constants::eV))
+
+      QUANTITY_DEFINE_SCALING_LITERALS(barn, corsika::units::si::sigma_d,
+                                       magnitude(corsika::units::si::constants::barn))
+
+      // phys::units::quantity<energy_d/mass_d> Joule2Kg = c2; // 1_Joule / 1_kg;
+
+      QUANTITY_DEFINE_SCALING_LITERALS(meter, length_d,
+                                       magnitude(corsika::units::si::constants::meter))
+
+      QUANTITY_DEFINE_SCALING_LITERALS(newton_second, corsika::units::si::momentum_d,
+                                       magnitude(corsika::units::si::newton_second))
+
+    } // namespace literals
+  }   // namespace units
+} // namespace phys
+
 // we want to call the operator<< without namespace... I think
 using namespace phys::units::io;
 

Framework/Units/testUnits.cc

@@ -32,7 +32,7 @@ TEST_CASE("PhysicalUnits", "[Units]") {
 
     LengthType l1 = 10_nm;
     l1 = l1;
-
+    
     LengthType arr0[5];
     arr0[0] = 5_m;
 
@@ -41,6 +41,9 @@ TEST_CASE("PhysicalUnits", "[Units]") {
     std::array<EnergyType, 4> arr2; // empty array
 
     [[maybe_unused]] std::array<EnergyType, 4> arr3 = {1_GeV, 1_eV, 5_MeV};
+
+    [[maybe_unused]] auto p1 = 10_newton_second;
+    REQUIRE(p1 == 10_newton_second);
   }
 
   SECTION("Powers in literal units") {
@@ -54,7 +57,8 @@ TEST_CASE("PhysicalUnits", "[Units]") {
     REQUIRE(1_mol / 1_amol == Approx(1e18));
     REQUIRE(1_K / 1_zK == Approx(1e21));
     REQUIRE(1_K / 1_yK == Approx(1e24));
-
+    REQUIRE(1_barn / 1_mbarn == Approx(1e3));
+    
     REQUIRE(1_A / 1_hA == Approx(1e-2));
     REQUIRE(1_m / 1_km == Approx(1e-3));
     REQUIRE(1_m / 1_Mm == Approx(1e-6));
@@ -76,6 +80,10 @@ TEST_CASE("PhysicalUnits", "[Units]") {
     REQUIRE(E2 == 40_GeV);
     REQUIRE(E2 / 1_GeV == Approx(40));
 
+    const MassType m = 1_kg;
+    const SpeedType v = 1_m / 1_s;
+    REQUIRE( m*v == 1_newton_second);
+    
     const double lgE = log10(E2 / 1_GeV);
     REQUIRE(lgE == Approx(log10(40.)));
 

Processes/CMakeLists.txt

 
 add_subdirectory (NullModel)
-#add_subdirectory (Sibyll)
+add_subdirectory (Sibyll)
 add_subdirectory (StackInspector)
 
 #add_custom_target(CORSIKAprocesses)

Processes/Sibyll/CMakeLists.txt

@@ -7,7 +7,7 @@ add_custom_command (
           sibyll_codes.dat          
           ${PROJECT_BINARY_DIR}/Framework/Particles/pythia_db.pkl
   WORKING_DIRECTORY
-          ${PROJECT_BINARY_DIR}/Processes/Sibyll/Particles/
+          ${PROJECT_BINARY_DIR}/Processes/Sibyll/
   COMMENT "Generate conversion tables for particle codes SIBYLL <-> CORSIKA"
   VERBATIM
   )
@@ -31,6 +31,21 @@ set (
 add_library (ProcessSibyll STATIC ${MODEL_SOURCES})
 CORSIKA_COPY_HEADERS_TO_NAMESPACE (ProcessSibyll ${MODEL_NAMESPACE} ${MODEL_HEADERS})
 
+# ....................................................
+# since Generated.inc is an automatically produced file in the build directory,
+# create a symbolic link into the source tree, so that it can be found and edited more easily
+# this is not needed for the build to succeed! .......
+add_custom_command (
+  OUTPUT  ${CMAKE_CURRENT_SOURCE_DIR}/Generated.inc
+  COMMAND ${CMAKE_COMMAND} -E create_symlink ${PROJECT_BINARY_DIR}/include/corsika/process/sibyll/Generated.inc ${CMAKE_CURRENT_SOURCE_DIR}/Generated.inc
+  COMMENT "Generate link in source-dir: ${CMAKE_CURRENT_SOURCE_DIR}/Generated.inc"
+  )
+add_custom_target (SourceDirLink2 DEPENDS ${PROJECT_BINARY_DIR}/Processes/Sibyll/Generated.inc)
+add_dependencies (ProcessSibyll SourceDirLink2)
+# .....................................................
+
+
+
 set_target_properties (
   ProcessSibyll
   PROPERTIES
@@ -42,7 +57,9 @@ set_target_properties (
 # target dependencies on other libraries (also the header onlys)
 target_link_libraries (
   ProcessSibyll
+  CORSIKAparticles
   CORSIKAunits
+  CORSIKAthirdparty
   )
 
 target_include_directories (

Processes/Sibyll/ParticleConversion.h

@@ -14,27 +14,42 @@
 
 #include <corsika/particles/ParticleProperties.h>
 
+#include <bitset2/bitset2.hpp>
+
 #include <map>
 
 namespace corsika::process::sibyll {
-  enum class Code : int8_t;
-  using PIDIntType = std::underlying_type<Code>::type;
 
-#include <corsika/processes/sibyll/Generated.inc>
+  enum class SibyllCode : int8_t;
+  using SibyllCodeIntType = std::underlying_type<SibyllCode>::type;
+
+#include <corsika/process/sibyll/Generated.inc>
+
+  bool KnownBySibyll(corsika::particles::Code pCode) {
+    return isKnown[static_cast<corsika::particles::CodeIntType>(pCode)];
+  }
 
-  bool handledBySibyll(corsika::particles::Code pCode) {
-    return handleable[static_cast<CodeIntType>(pCode)];
+  bool CanInteract(corsika::particles::Code pCode) {
+    return canInteract[static_cast<corsika::particles::CodeIntType>(pCode)];
   }
 
-  Code constexpr ConvertToSibyll(corsika::particles::Code pCode) {
+  SibyllCode constexpr ConvertToSibyll(corsika::particles::Code pCode) {
     //~ assert(handledBySibyll(pCode));
-    return static_cast<Code>(corsika2sibyll[static_cast<CodeIntType>(pCode)]);
+    return static_cast<SibyllCode>(corsika2sibyll[static_cast<corsika::particles::CodeIntType>(pCode)]);
+  }
+  
+  corsika::particles::Code constexpr ConvertFromSibyll(SibyllCode pCode) {
+    return sibyll2corsika[static_cast<SibyllCodeIntType>(pCode) - minSibyll];
   }
 
-  corsika::particles::Code constexpr ConvertFromSibyll(Code pCode) {
-    return sibyll2corsika[static_cast<PIDIntType>(pCode) - minSibyll];
+  int ConvertToSibyllRaw(corsika::particles::Code pCode){
+    return  (int)static_cast<corsika::process::sibyll::SibyllCodeIntType>( corsika::process::sibyll::ConvertToSibyll( pCode ) );
   }
 
+  int GetSibyllXSCode(corsika::particles::Code pCode) {
+    return corsika2sibyllXStype[static_cast<corsika::particles::CodeIntType>(pCode)];
+  }
+  
 } // namespace corsika::process::sibyll
 
 #endif

Processes/Sibyll/code_generator.py

@@ -12,12 +12,14 @@
 import pickle, sys, itertools
 
 
+
 # loads the pickled pythia_db (which is an OrderedDict)
 def load_pythiadb(filename):
     with open(filename, "rb") as f:
         pythia_db = pickle.load(f)
     return pythia_db
-    
+
+
 
 # 
 def read_sibyll_codes(filename, pythia_db):
@@ -25,68 +27,91 @@ def read_sibyll_codes(filename, pythia_db):
         for line in f:
             line = line.strip()
             if line[0] == '#':
-                continue
-            identifier, sib_code = line.split()
+                continue            
+            identifier, sib_code, canInteractFlag, xsType = line.split()
             try:
                 pythia_db[identifier]["sibyll_code"] = int(sib_code)
+                pythia_db[identifier]["sibyll_canInteract"] = int(canInteractFlag)
+                pythia_db[identifier]["sibyll_xsType"] = int(xsType)
             except KeyError as e:
                 raise Exception("Identifier '{:s}' not found in pythia_db".format(identifier))
+
+
             
 
 # generates the enum to access sibyll particles by readable names
 def generate_sibyll_enum(pythia_db):
-    output = "enum class PID : int16_t {\n"
-    for identifier, d in pythia_db.items():
-        if d.get('sibyll_code') != None:
-            output += "  {:s} = {:d},\n".format(identifier, d['sibyll_code'])
+    output = "enum class SibyllCode : int8_t {\n"
+    for identifier, pData in pythia_db.items():
+        if pData.get('sibyll_code') != None:
+            output += "  {:s} = {:d},\n".format(identifier, pData['sibyll_code'])
     output += "};\n"
     return output
-            
+
+
 
 # generates the look-up table to convert corsika codes to sibyll codes
 def generate_corsika2sibyll(pythia_db):    
-    string = "std::array<PIDIntType, {:d}> constexpr corsika2sibyll = {{".format(len(pythia_db))
-    for identifier, d in pythia_db.items():
-        sibCode = d.get("sibyll_code", 0)
+    string = "std::array<SibyllCodeIntType, {:d}> constexpr corsika2sibyll = {{\n".format(len(pythia_db))
+    for identifier, pData in pythia_db.items():
+        sibCode = pData.get("sibyll_code", 0)
         string += "  {:d}, // {:s}\n".format(sibCode, identifier if sibCode else identifier + " (not implemented in SIBYLL)")
     string += "};\n"
     return string
     
 
+
+# generates the look-up table to convert corsika codes to sibyll codes
+def generate_corsika2sibyll_xsType(pythia_db):    
+    string = "std::array<int, {:d}> constexpr corsika2sibyllXStype = {{\n".format(len(pythia_db))
+    for identifier, pData in pythia_db.items():
+        sibCodeXS = pData.get("sibyll_xsType", -1)
+        string += "  {:d}, // {:s}\n".format(sibCodeXS, identifier if sibCodeXS else identifier + " (not implemented in SIBYLL)")
+    string += "};\n"
+    return string
+
+
 # generates the look-up table to convert sibyll codes to corsika codes    
-def generate_sibyll2corsika(pythia_db):
-    d = {}
-    for identifier, p in pythia_db.items():
-        if 'sibyll_code' in p:
-            sib_code = p['sibyll_code']
-            corsika_code = p['ngc_code']
-            d[sib_code] = (corsika_code, identifier)
-            
-    string = "std::array<corsika::particles::CodeIntType, {:d}> sibyll2corsika = {{\n".format(len(d))
+def generate_sibyll2corsika(pythia_db) :
+    string = ""
+    
+    minID = 0
+    for identifier, pData in pythia_db.items() :
+        if 'sibyll_code' in pData:
+            minID = min(minID, pData['sibyll_code'])
+
+    string += "SibyllCodeIntType constexpr minSibyll = {:d};\n\n".format(minID)
+
+    pDict = {}
+    for identifier, pData in pythia_db.items() :
+        if 'sibyll_code' in pData:
+            sib_code = pData['sibyll_code'] - minID
+            pDict[sib_code] = identifier
+    
+    nPart = max(pDict.keys()) - min(pDict.keys()) + 1
+    string += "std::array<corsika::particles::Code, {:d}> sibyll2corsika = {{\n".format(nPart)
     
-    for k in range(min(d.keys()), max(d.keys())+1):
-        if k in d:
-            corsika_code = d[k][0]
-            identifier = d[k][1]
+    for iPart in range(nPart) :
+        if iPart in pDict:
+            identifier = pDict[iPart]
         else:
-            corsika_code = 0
-            identifier = ""
-        string += "  {:d}, // {:s}\n".format(corsika_code, identifier)
+            identifier = "Unknown"
+        string += "  corsika::particles::Code::{:s}, \n".format(identifier)
     
     string += "};\n"
-    string += "PIDIntType constexpr minSibyll = {:d};\n".format(min(d.keys()))
     return string
-    
+
+
 
 # generates the bitset for the flag whether Sibyll knows the particle
-def generate_handles_particle(pythia_db):
+def generate_known_particle(pythia_db):
     num_particles = len(pythia_db)
     num_bytes = num_particles // 32 + 1
-    string = "Bitset2::bitset2<{:d}> constexpr handleable{{ std::array<uint32_t, {:d}>{{{{\n".format(num_particles, num_bytes)
+    string = "Bitset2::bitset2<{:d}> constexpr isKnown{{ std::array<uint32_t, {:d}>{{{{\n".format(num_particles, num_bytes)
     
     numeric = 0
-    for identifier, d in reversed(pythia_db.items()):
-        handledBySibyll = int("sibyll_code" in d) & 0x1
+    for identifier, pData in reversed(pythia_db.items()):
+        handledBySibyll = int("sibyll_code" in pData) & 0x1
         numeric = (numeric << 1) | handledBySibyll
     
     while numeric != 0:
@@ -96,8 +121,34 @@ def generate_handles_particle(pythia_db):
         
     string += "}}};\n"
     return string
+
+
+
+# generates the bitset for the flag whether Sibyll can use particle as projectile
+def generate_interacting_particle(pythia_db):
+    num_particles = len(pythia_db)
+    num_bytes = num_particles // 32 + 1
+    string = "Bitset2::bitset2<{:d}> constexpr canInteract{{ std::array<uint32_t, {:d}>{{{{\n".format(num_particles, num_bytes)
+    #string = "std::array<bool, {:d}> constexpr corsika2sibyll = {{\n".format(num_particles)
+    
+    numeric = 0    
+    for identifier, pData in reversed(pythia_db.items()):
+        can = 0
+        if 'sibyll_canInteract' in pData:
+            if pData['sibyll_canInteract'] > 0:
+                can = 0x1
+        numeric = (numeric << 1) | can
+    
+    while numeric != 0:
+        low = numeric & 0xFFFFFFFF
+        numeric = numeric >> 32
+        string += "  0x{:0x},\n".format(low)
+    
+    string += "}}};\n"
+    return string
     
 
+
 if __name__ == "__main__":
     if len(sys.argv) != 3:
         print("usage: {:s} <pythia_db.pkl> <sibyll_codes.dat>".format(sys.argv[0]), file=sys.stderr)
@@ -106,12 +157,15 @@ if __name__ == "__main__":
     print("code_generator.py for SIBYLL")
     
     pythia_db = load_pythiadb(sys.argv[1])
-
     read_sibyll_codes(sys.argv[2], pythia_db)
     
+    print (str(pythia_db))
+    
     with open("Generated.inc", "w") as f:
         print("// this file is automatically generated\n// edit at your own risk!\n", file=f)
         print(generate_sibyll_enum(pythia_db), file=f)
         print(generate_corsika2sibyll(pythia_db), file=f)
-        print(generate_handles_particle(pythia_db), file=f)
+        print(generate_known_particle(pythia_db), file=f)
         print(generate_sibyll2corsika(pythia_db), file=f)
+        print(generate_interacting_particle(pythia_db), file=f)
+        print(generate_corsika2sibyll_xsType(pythia_db), file=f)

Processes/Sibyll/sibyll_codes.dat

-Electron 3
-Positron 2
-NuE 15
-NuEBar 16
-MuMinus 5
-MuPlus 4
-NuMu 17
-NuMuBar 18
-TauMinus 91
-TauPlus 90
-NuTau 92
-NuTauBar 93
-Gamma 1
-Pi0 6
-Rho0 27
-K0Long 11
-PiPlus 7
-PiMinus 8
-RhoPlus 25
-RhoMinus 26
-Eta 23
-Omega 32
-K0Short 12
-KStar0 30
-KStar0Bar 31
-KPlus 9
-KMinus 10
-KStarPlus 28
-KStarMinus 29
-DPlus 59
-DMinus 60
-DStarPlus 78
-DStarMinus 79
-D0 71
-D0Bar 72
-DStar0 80
-DStar0Bar 81
-DsPlus 74
-DsMinus 75
-DStarSPlus 76
-DStarSMinus 77
-EtaC 73
-Neutron 14
-AntiNeutron -14
-Delta0 42
-Delta0Bar -42
-Proton 13
-AntiProton -13
-DeltaPlus 41
-DeltaMinusBar -41
-DeltaPlusPlus 40
-DeltaMinusMinusBar -40
-SigmaMinus 36
-SigmaPlusBar -36
-Lambda0 39
-Lambda0Bar -39
-Sigma0 35
-Sigma0Bar -35
-SigmaPlus 34
-SigmaMinusBar -34
-XiMinus 38
-XiPlusBar -38
-Xi0 37
-Xi0Bar -37
-OmegaMinus 49
-OmegaPlusBar -49
-SigmaC0 86
-SigmaC0Bar -86
-SigmaStarC0 96
-SigmaStarC0Bar -96
-LambdaCPlus 89
-LambdaCMinusBar -89
-XiC0 88
-XiC0Bar -88
-SigmaCPlus 85
-SigmaCMinusBar -85
-SigmaStarCPlus 95
-SigmaStarCMinusBar -95
-SigmaCPlusPlus 84
-SigmaCMinusMinusBar -84
-SigmaStarCPlusPlus 94
-SigmaStarCMinusMinusBar -94
-XiCPlus 87
-XiCMinusBar -87
-OmegaC0 99
-OmegaC0Bar -99
-Jpsi 83
-#Unknown 0
+# input file for particle conversion to/from SIBYLL
+# the format of this file is: "corsika-identifier" "sibyll-id" "can-interact-in-sibyll" "cross-section-type"
+Electron 3 0 0
+Positron 2 0 0
+NuE 15 0 0
+NuEBar 16 0 0
+MuMinus 5 0 0
+MuPlus 4 0 0
+NuMu 17 0 0
+NuMuBar 18 0 0
+TauMinus 91 0 0
+TauPlus 90 0 0
+NuTau 92 0 0
+NuTauBar 93 0 0
+Gamma 1 0 0
+Pi0 6 1 2
+# rho0 could interact but sibyll has no cross section/interaction length. was used for gamma had int
+Rho0 27 0 0
+K0Long 11 1 3
+PiPlus 7 1 2
+PiMinus 8 1 2
+RhoPlus 25 0 0
+RhoMinus 26 0 0
+Eta 23 0 0
+Omega 32 0 0
+K0Short 12 1 3
+KStar0 30 0 0
+KStar0Bar 31 0 0
+KPlus 9 1 3
+KMinus 10 1 3
+KStarPlus 28 0 0
+KStarMinus 29 0 0
+DPlus 59 1 0
+DMinus 60 1 0
+DStarPlus 78 0 0
+DStarMinus 79 0 0
+D0 71 1 0
+D0Bar 72 1 0
+DStar0 80 0 0
+DStar0Bar 81 0 0
+DsPlus 74 1 0
+DsMinus 75 1 0
+DStarSPlus 76 0 0
+DStarSMinus 77 0 0
+EtaC 73 0 0
+Neutron 14 1 1
+AntiNeutron -14 1 1
+Delta0 42 0 0
+Delta0Bar -42 0 0
+Proton 13 1 1
+AntiProton -13 1 1
+DeltaPlus 41 0 0
+DeltaMinusBar -41 0 0
+DeltaPlusPlus 40 0 0
+DeltaMinusMinusBar -40 0 0
+SigmaMinus 36 1 1
+SigmaPlusBar -36 1 1
+Lambda0 39 1 1
+Lambda0Bar -39 1 1
+Sigma0 35 1 1
+Sigma0Bar -35 1 1
+SigmaPlus 34 1 1
+SigmaMinusBar -34 1 1
+XiMinus 38 1 1
+XiPlusBar -38 1 1
+Xi0 37 1 1
+Xi0Bar -37 1 1
+OmegaMinus 49 0 0
+OmegaPlusBar -49 0 0
+SigmaC0 86 0 0
+SigmaC0Bar -86 0 0
+SigmaStarC0 96 0 0
+SigmaStarC0Bar -96 0 0
+LambdaCPlus 89 1 1
+LambdaCMinusBar -89 1 1
+XiC0 88 1 1
+XiC0Bar -88 1 1
+SigmaCPlus 85 0 0
+SigmaCMinusBar -85 0 0
+SigmaStarCPlus 95 0 0
+SigmaStarCMinusBar -95 0 0
+SigmaCPlusPlus 84 0 0
+SigmaCMinusMinusBar -84 0 0
+SigmaStarCPlusPlus 94 0 0
+SigmaStarCMinusMinusBar -94 0 0
+XiCPlus 87 1 1
+XiCMinusBar -87 1 1
+OmegaC0 99 0 0
+OmegaC0Bar -99 0 0
+Jpsi 83 0 0
+#Unknown 0 0 0

Processes/Sibyll/testSibyll.cc

@@ -17,24 +17,45 @@
                           // cpp file
 #include <catch2/catch.hpp>
 
+#include <iostream>
+using namespace std;
 using namespace corsika;
 
 TEST_CASE("Sibyll", "[processes]") {
 
   SECTION("Sibyll -> Corsika") {
     REQUIRE(corsika::particles::Electron::GetCode() ==
-            process::sibyll::ConvertFromSibyll(process::sibyll::Code::Electron));
+            process::sibyll::ConvertFromSibyll(process::sibyll::SibyllCode::Electron));
   }
 
   SECTION("Corsika -> Sibyll") {
     REQUIRE(process::sibyll::ConvertToSibyll(corsika::particles::Electron::GetCode()) ==
-            process::sibyll::Code::Electron);
+            process::sibyll::SibyllCode::Electron);
+    REQUIRE(process::sibyll::ConvertToSibyllRaw(corsika::particles::Proton::GetCode()) ==
+            13);
   }
 
-  SECTION("handledBySibyll") {
-    REQUIRE(process::sibyll::handledBySibyll(corsika::particles::Electron::GetCode()));
+  SECTION("KnownBySibyll") {
+    REQUIRE(process::sibyll::KnownBySibyll(corsika::particles::Electron::GetCode()));
 
     REQUIRE_FALSE(
-        process::sibyll::handledBySibyll(corsika::particles::XiPrimeC0::GetCode()));
+        process::sibyll::KnownBySibyll(corsika::particles::XiPrimeC0::GetCode()));
+  }
+
+  SECTION("canInteractInSibyll") {
+
+    REQUIRE(process::sibyll::CanInteract(corsika::particles::Proton::GetCode()));
+    REQUIRE(process::sibyll::CanInteract(corsika::particles::Code::XiCPlus));
+
+    REQUIRE_FALSE(process::sibyll::CanInteract(corsika::particles::Electron::GetCode()));
+    REQUIRE_FALSE(process::sibyll::CanInteract(corsika::particles::SigmaC0::GetCode()));
+  }
+  
+  SECTION("cross-section type") {
+
+    REQUIRE(process::sibyll::GetSibyllXSCode(corsika::particles::Code::Electron) == 0);
+    REQUIRE(process::sibyll::GetSibyllXSCode(corsika::particles::Code::K0Long) == 3);
+    REQUIRE(process::sibyll::GetSibyllXSCode(corsika::particles::Code::SigmaPlus) == 1);
+    REQUIRE(process::sibyll::GetSibyllXSCode(corsika::particles::Code::PiMinus) == 2);
   }
 }

Processes/StackInspector/StackInspector.cc

@@ -37,12 +37,15 @@ process::EProcessReturn StackInspector<Stack, Trajectory>::DoContinuous(Particle
   if (!fReport) return EProcessReturn::eOk;
   [[maybe_unused]] int i = 0;
   EnergyType Etot = 0_GeV;
+
   for (auto& iterP : s) {
     EnergyType E = iterP.GetEnergy();
     Etot += E;
+    std::cout << " particle data: " << i++ << ", id=" << iterP.GetPID()<< ", en=" << iterP.GetEnergy() / 1_GeV << " | " << std::endl;
   }
   countStep++;
-  cout << countStep << " " << s.GetSize() << " " << Etot / 1_GeV << " " << endl;
+  // cout << "#=" << countStep << " " << s.GetSize() << " " << Etot/1_GeV << endl;
+  cout << "call: " << countStep << " stack size: " << s.GetSize() << " energy: " << Etot / 1_GeV << endl;
   return EProcessReturn::eOk;
 }