added cascade example with sibyll, primitive

CMakeLists.txt

@@ -21,6 +21,8 @@ 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/geometry/LineTrajectory.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>
+
+using namespace corsika;
+using namespace corsika::process;
+using namespace corsika::units;
+using namespace corsika::particles;
+using namespace corsika::random;
+
+#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
+    int kBeam   = 1;
+    // target nuclei: A < 18
+    int kTarget = 0.4*16 + 0.6*14;
+    double beamEnergy =  p.GetEnergy() / 1_GeV; 
+    std::cout << "ProcessSplit: " << "MinStep: en: " << beamEnergy << " pid:" << kBeam << std::endl;
+    double prodCrossSection,dummy;
+    sib_sigma_hnuc_(kBeam, kTarget, beamEnergy, prodCrossSection, dummy );
+    std::cout << "ProcessSplit: " << "MinStep: sibyll return: " << prodCrossSection << std::endl;
+    CrossSectionType sig = prodCrossSection  / 1000. * barn;
+    std::cout << "ProcessSplit: " << "MinStep: CrossSection= " << sig << std::endl;
+
+    // calculate interaction length in medium
+
+    // pick random step lenth
+    
+    return prodCrossSection;
+  }
+
+  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 {
+    // get energy of particle from stack
+    // stack is in GeV in lab. frame
+    // convert to GeV in cm. frame
+    EnergyType E   = p.GetEnergy();
+    double Ecm = sqrt( 2. * E * 0.93827_GeV ) / 1_GeV ;
+    // FOR NOW: set beam to proton
+    int kBeam   = 13; //p.GetPID();
+    // FOR NOW: set target to proton
+    int kTarget = 1; //p.GetPID();
+    std::cout << "ProcessSplit: " << " DoDiscrete: E(GeV):" << E / 1_GeV << " Ecm(GeV): " << Ecm << std::endl;
+    if (E < 8.5_GeV || Ecm < 10. ) {
+      std::cout << "ProcessSplit: " << " DoDiscrete: dropping particle.." << std::endl;
+      p.Delete();
+      fCount++;
+    } else {
+      //p.SetEnergy(E / 2);
+      //s.NewParticle().SetEnergy(E / 2);
+      
+      // running sibyll
+      sibyll_( kBeam, kTarget, Ecm);
+      int print_unit = 6;
+      sib_list_( print_unit );
+      
+      // delete current particle
+      p.Delete();
+      
+      // add particles from sibyll to stack
+      for(int i=0; i<s_plist_.np; ++i){
+	
+	s.NewParticle().SetEnergy( s_plist_.p[3][i] * 1_GeV );
+      }
+    }
+  }
+
+  void Init() //{ fCount = 0; }
+  {
+    fCount = 0;
+
+    // initialize random numbers for sibyll
+    // FOR NOW USE SIBYLL INTERNAL !!!
+    rnd_ini_();
+    
+    // corsika::random::RNGManager rmng;
+    // 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<5; ++i)
+      std::cout << i << " " << s_rndm_(a) << std::endl;
+    
+    //initialize Sibyll
+    sibyll_ini_();
+  }
+  
+  int GetCount() { return fCount; }
+
+private:
+};
+
+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,7 @@ set (
 set (
   CORSIKAcascade_HEADERS
   Cascade.h
+  sibyll2.3c.h
   )
 
 #set (
@@ -16,9 +17,23 @@ set (
 #  Cascade.cc
 #  )
 
+set (
+  CORSIKAsibyll_NAMESPACE
+  corsika/cascade
+  )
+
+set (
+  CORSIKAsibyll_SOURCES
+  sibyll2.3c.f
+  rndm_dbl.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,7 +66,9 @@ add_executable (
 
 target_link_libraries (
   testCascade
-#  CORSIKAutls
+  #  CORSIKAutls
+  CORSIKArandom
+  CORSIKAsibyll
   CORSIKAcascade
   ProcessStackInspector
   CORSIKAprocesses

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_;
+
+
+  // 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&);
+  
+  double s_rndm_(int&);
+
+  // phojet random generator setup
+  void pho_rndin_(int&,int&,int&,int&);
+}
+#endif
+

Framework/Units/PhysicalConstants.h

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

Framework/Units/PhysicalUnits.h

@@ -38,6 +38,8 @@ namespace corsika::units::si {
   using EnergyType = phys::units::quantity<phys::units::energy_d, double>;
   using MassType = phys::units::quantity<phys::units::mass_d, double>;
 
+  using CrossSectionType = phys::units::quantity<phys::units::area_d, double>;
+  
 } // end namespace corsika::units::si
 
 // we want to call the operator<< without namespace... I think