additional changes

CMakeLists.txt

@@ -17,9 +17,6 @@ set (CMAKE_INSTALL_MESSAGE LAZY)
 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 -Wno-ignored-qualifiers)
-
 set (CMAKE_CXX_STANDARD 17)
 enable_testing ()
 set (CTEST_OUTPUT_ON_FAILURE 1)
@@ -39,9 +36,10 @@ if(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES)
     "Debug" "Release" "MinSizeRel" "RelWithDebInfo")
 endif()
 
-#set(CMAKE_CXX_FLAGS "-Wall -Wextra")
-set(CMAKE_CXX_FLAGS_DEBUG "-O0 -g")
-set(CMAKE_CXX_FLAGS_RELEASE "-O3")
+# enable warnings and disallow non-standard language
+set(CMAKE_CXX_FLAGS "-Wall -pedantic -Wextra -Wno-ignored-qualifiers")
+set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS} -O0 -g")
+set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS} -O3")
 
 # unit testing coverage, does not work yet
 #include (CodeCoverage)

Documentation/Examples/CMakeLists.txt

@@ -25,7 +25,7 @@ add_executable (cascade_example cascade_example.cc)
 target_compile_options(cascade_example PRIVATE -g) # do not skip asserts
 target_link_libraries (cascade_example SuperStupidStack CORSIKAunits CORSIKAlogging
    CORSIKArandom
-  CORSIKAsibyll
+   ProcessSibyll
   CORSIKAcascade
   ProcessStackInspector
   CORSIKAprocesses

Documentation/Examples/cascade_example.cc

@@ -20,14 +20,10 @@
 #include <corsika/environment/HomogeneousMedium.h>
 #include <corsika/environment/NuclearComposition.h>
 
-#include <corsika/random/RNGManager.h>
-
-#include <corsika/cascade/SibStack.h>
-#include <corsika/cascade/sibyll2.3c.h>
 #include <corsika/geometry/Sphere.h>
-#include <corsika/process/sibyll/ParticleConversion.h>
 
-#include <corsika/process/sibyll/ProcessDecay.h>
+#include <corsika/process/sibyll/Decay.h>
+#include <corsika/process/sibyll/Interaction.h>
 
 #include <corsika/units/PhysicalUnits.h>
 
@@ -47,7 +43,6 @@ using namespace corsika::environment;
 using namespace std;
 using namespace corsika::units::si;
 
-static int fCount = 0;
 static EnergyType fEnergy = 0. * 1_GeV;
 
 // FOR NOW: global static variables for ParticleCut process
@@ -121,7 +116,7 @@ public:
   }
 
   template <typename Particle>
-  double MinStepLength(Particle& p, setup::Trajectory&) const {
+  double MaxStepLength(Particle& p, setup::Trajectory&) const {
     const Code pid = p.GetPID();
     if (isEmParticle(pid) || isInvisible(pid)) {
       cout << "ProcessCut: MinStep: next cut: " << 0. << endl;
@@ -173,7 +168,6 @@ public:
     } else if (isBelowEnergyCut(p)) {
       cout << "removing low en. particle..." << endl;
       fEnergy += p.GetEnergy();
-      fCount += 1;
       p.Delete();
     }
   }
@@ -206,296 +200,6 @@ public:
 private:
 };
 
-class ProcessSplit : public corsika::process::InteractionProcess<ProcessSplit> {
-public:
-  ProcessSplit() {}
-
-  void setTrackedParticlesStable() const {
-    /*
-      Sibyll is hadronic generator
-      only hadrons decay
-     */
-    // set particles unstable
-    corsika::process::sibyll::setHadronsUnstable();
-    // make tracked particles stable
-    std::cout << "ProcessSplit: setting tracked hadrons stable.." << std::endl;
-    setup::Stack ds;
-    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
-      s_csydec_.idb[s_id - 1] = (-1) * abs(s_csydec_.idb[s_id - 1]);
-      p.Delete();
-    }
-  }
-
-  template <typename Particle, typename Track>
-  double GetInteractionLength(Particle& p, Track&) const {
-
-    // coordinate system, get global frame of reference
-    CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
-
-    const Code corsikaBeamId = p.GetPID();
-
-    // beam particles for sibyll : 1, 2, 3 for p, pi, k
-    // read from cross section code table
-    int kBeam = process::sibyll::GetSibyllXSCode(corsikaBeamId);
-
-    bool kInteraction = process::sibyll::CanInteract(corsikaBeamId);
-
-    /*
-       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 = 16;
-
-    EnergyType Etot = p.GetEnergy() + kTarget * corsika::particles::Proton::GetMass();
-    super_stupid::MomentumVector Ptot(
-        rootCS, {0.0_newton_second, 0.0_newton_second, 0.0_newton_second});
-    // FOR NOW: assume target is at rest
-    super_stupid::MomentumVector pTarget(
-        rootCS, {0.0_newton_second, 0.0_newton_second, 0.0_newton_second});
-    Ptot += p.GetMomentum();
-    Ptot += pTarget;
-    // calculate cm. energy
-    EnergyType sqs = sqrt(Etot * Etot - Ptot.squaredNorm() * si::constants::cSquared);
-    double Ecm = sqs / 1_GeV;
-
-    std::cout << "ProcessSplit: "
-              << "MinStep: input en: " << p.GetEnergy() / 1_GeV << endl
-              << " beam can interact:" << kBeam << endl
-              << " beam XS code:" << kBeam << endl
-              << " beam pid:" << p.GetPID() << endl
-              << " target mass number:" << kTarget << std::endl;
-
-    double int_length = 0;
-    if (kInteraction) {
-
-      double prodCrossSection, dummy, dum1, dum2, dum3, dum4;
-      double dumdif[3];
-
-      if (kTarget == 1)
-        sib_sigma_hp_(kBeam, Ecm, dum1, dum2, prodCrossSection, dumdif, dum3, dum4);
-      else
-        sib_sigma_hnuc_(kBeam, kTarget, Ecm, 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
-      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;
-    } else
-      int_length = std::numeric_limits<double>::infinity();
-
-    /*
-      what are the units of the output? slant depth or 3space length?
-
-    */
-    return int_length;
-    //
-    // int a = 0;
-    // const double next_step = -int_length * log(s_rndm_(a));
-    // std::cout << "ProcessSplit: "
-    //        << "next step (g/cm2): " << next_step << std::endl;
-    // return next_step;
-  }
-
-  template <typename Particle, typename Track, typename Stack>
-  EProcessReturn DoContinuous(Particle&, Track&, Stack&) const {
-    return EProcessReturn::eOk;
-  }
-
-  template <typename Particle, typename Stack>
-  void DoInteraction(Particle& p, Stack& s) const {
-    cout << "ProcessSibyll: "
-         << "DoInteraction: " << p.GetPID() << " interaction? "
-         << process::sibyll::CanInteract(p.GetPID()) << endl;
-    if (process::sibyll::CanInteract(p.GetPID())) {
-      cout << "defining coordinates" << endl;
-      // coordinate system, get global frame of reference
-      CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
-
-      QuantityVector<length_d> const coordinates{0_m, 0_m, 0_m};
-      Point pOrig(rootCS, coordinates);
-
-      /*
-         the target should be defined by the Environment,
-         ideally as full particle object so that the four momenta
-         and the boosts can be defined..
-
-         here we need: GetTargetMassNumber() or GetTargetPID()??
-                       GetTargetMomentum() (zero in EAS)
-      */
-      // FOR NOW: set target to proton
-      int kTarget = 1; // env.GetTargetParticle().GetPID();
-
-      cout << "defining target momentum.." << endl;
-      // FOR NOW: target is always at rest
-      const EnergyType Etarget = 0. * 1_GeV + corsika::particles::Proton::GetMass();
-      const auto pTarget = super_stupid::MomentumVector(
-          rootCS, 0. * 1_GeV / si::constants::c, 0. * 1_GeV / si::constants::c,
-          0. * 1_GeV / si::constants::c);
-      cout << "target momentum (GeV/c): "
-           << pTarget.GetComponents() / 1_GeV * si::constants::c << endl;
-      cout << "beam momentum (GeV/c): "
-           << p.GetMomentum().GetComponents() / 1_GeV * si::constants::c << endl;
-
-      // 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)
-      */
-      // total energy: E_beam + E_target
-      // in lab. frame: E_beam + m_target*c**2
-      EnergyType E = p.GetEnergy();
-      EnergyType Etot = E + Etarget;
-      // total momentum
-      super_stupid::MomentumVector Ptot = p.GetMomentum(); // + pTarget;
-      // invariant mass, i.e. cm. energy
-      EnergyType Ecm = sqrt(Etot * Etot -
-                            Ptot.squaredNorm() *
-                                si::constants::cSquared); // sqrt( 2. * E * 0.93827_GeV );
-      /*
-       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
-     */
-      const double gamma = Etot / Ecm;
-      const auto gambet = Ptot / (Ecm / si::constants::c);
-
-      std::cout << "ProcessSplit: "
-                << " DoDiscrete: gamma:" << gamma << endl;
-      std::cout << "ProcessSplit: "
-                << " DoDiscrete: gambet:" << gambet.GetComponents() << endl;
-
-      int kBeam = process::sibyll::ConvertToSibyllRaw(p.GetPID());
-
-      std::cout << "ProcessSplit: "
-                << " DoInteraction: E(GeV):" << E / 1_GeV << " Ecm(GeV): " << Ecm / 1_GeV
-                << std::endl;
-      if (E < 8.5_GeV || Ecm < 10_GeV) {
-        std::cout << "ProcessSplit: "
-                  << " DoInteraction: 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
-        setTrackedParticlesStable();
-        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;
-
-        // SibStack does not know about momentum yet so we need counter to access momentum
-        // array in Sibyll
-        int i = -1;
-        super_stupid::MomentumVector Ptot_final(
-            rootCS, {0.0_newton_second, 0.0_newton_second, 0.0_newton_second});
-        for (auto& psib : ss) {
-          ++i;
-          // skip particles that have decayed in Sibyll
-          if (abs(s_plist_.llist[i]) > 100) continue;
-
-          // transform energy to lab. frame, primitve
-          // compute beta_vec * p_vec
-          // arbitrary Lorentz transformation based on sibyll routines
-          const auto gammaBetaComponents = gambet.GetComponents();
-          const auto pSibyllComponents = psib.GetMomentum().GetComponents();
-          EnergyType en_lab = 0. * 1_GeV;
-          MomentumType p_lab_components[3];
-          en_lab = psib.GetEnergy() * gamma;
-          EnergyType pnorm = 0. * 1_GeV;
-          for (int j = 0; j < 3; ++j)
-            pnorm += (pSibyllComponents[j] * gammaBetaComponents[j] * si::constants::c) /
-                     (gamma + 1.);
-          pnorm += psib.GetEnergy();
-
-          for (int j = 0; j < 3; ++j) {
-            p_lab_components[j] = pSibyllComponents[j] - (-1) * pnorm *
-                                                             gammaBetaComponents[j] /
-                                                             si::constants::c;
-            en_lab -=
-                (-1) * pSibyllComponents[j] * gammaBetaComponents[j] * si::constants::c;
-          }
-
-          // add to corsika stack
-          auto pnew = s.NewParticle();
-          pnew.SetEnergy(en_lab);
-          pnew.SetPID(process::sibyll::ConvertFromSibyll(psib.GetPID()));
-
-          corsika::geometry::QuantityVector<momentum_d> p_lab_c{
-              p_lab_components[0], p_lab_components[1], p_lab_components[2]};
-          pnew.SetMomentum(super_stupid::MomentumVector(rootCS, p_lab_c));
-          Ptot_final += pnew.GetMomentum();
-        }
-        // cout << "tot. momentum final (GeV/c): " << Ptot_final.GetComponents() / 1_GeV *
-        // si::constants::c << endl;
-      }
-    }
-  }
-
-  void Init() {
-    fCount = 0;
-
-    corsika::random::RNGManager& rmng = corsika::random::RNGManager::GetInstance();
-
-    rmng.RegisterRandomStream("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_();
-
-    setTrackedParticlesStable();
-  }
-
-  int GetCount() { return fCount; }
-  EnergyType GetEnergy() { return fEnergy; }
-
-private:
-};
-
-double s_rndm_(int&) {
-  static corsika::random::RNG& rmng =
-      corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");
-  ;
-  return rmng() / (double)rmng.max();
-}
 
 int main() {
   corsika::environment::Environment env; // dummy environment
@@ -520,8 +224,8 @@ int main() {
   tracking_line::TrackingLine<setup::Stack> tracking(env);
   stack_inspector::StackInspector<setup::Stack> p0(true);
 
-  ProcessSplit p1;
-  corsika::process::sibyll::ProcessDecay p2;
+  corsika::process::sibyll::Interaction/*<setup::Stack>,setup::Trajectory>*/ p1;
+  corsika::process::sibyll::Decay p2;
   ProcessEMCut p3;
   const auto sequence = /*p0 +*/ p1 + p2 + p3;
   setup::Stack stack;
@@ -543,8 +247,10 @@ int main() {
   EAS.Init();
   EAS.Run();
   cout << "Result: E0=" << E0 / 1_GeV
-       << "GeV, particles below energy threshold =" << p1.GetCount() << endl;
-  cout << "total energy below threshold (GeV): " << p1.GetEnergy() / 1_GeV << std::endl;
+    //<< "GeV, particles below energy threshold =" << p1.GetCount()
+       << endl;
+  cout << "total energy below threshold (GeV): " //<< p1.GetEnergy() / 1_GeV
+       << std::endl;
   p3.ShowResults();
   cout << "total energy (GeV): "
        << (p3.GetCutEnergy() + p3.GetInvEnergy() + p3.GetEmEnergy()) / 1_GeV << endl;

Framework/Cascade/CMakeLists.txt

@@ -9,33 +9,10 @@ set (
 set (
   CORSIKAcascade_HEADERS
   Cascade.h
-  sibyll2.3c.h
-  SibStack.h
   )
 
-#set (
-#  CORSIKAcascade_SOURCES
-#  TrackingStep.cc
-#  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 (
@@ -70,7 +47,7 @@ target_link_libraries (
   testCascade
   #  CORSIKAutls
   CORSIKArandom
-  CORSIKAsibyll
+  ProcessSibyll
   CORSIKAcascade
   ProcessStackInspector
   CORSIKAstackinterface

Framework/ProcessSequence/ContinuousProcess.h

@@ -33,8 +33,8 @@ namespace corsika::process {
 
     // here starts the interface part
     // -> enforce derived to implement DoContinuous...
-    template <typename D, typename T, typename S>
-    inline EProcessReturn DoContinuous(D&, T&, S&) const;
+    template <typename P, typename T, typename S>
+    inline EProcessReturn DoContinuous(P&, T&, S&) const;
   };
 
 } // namespace corsika::process

Framework/ProcessSequence/InteractionProcess.h

@@ -34,11 +34,11 @@ namespace corsika::process {
 
     /// here starts the interface-definition part
     // -> enforce derived to implement DoInteraction...
-    template <typename Particle, typename Stack>
-    inline EProcessReturn DoInteraction(Particle&, Stack&) const;
+    template <typename P, typename S>
+    inline EProcessReturn DoInteraction(P&, S&) const;
 
-    template <typename Particle, typename Track>
-    inline double GetInteractionLength(Particle& p, Track& t) const;
+    template <typename P, typename T>
+    inline double GetInteractionLength(P&, T&) const;
 
     template <typename Particle, typename Track>
     inline double GetInverseInteractionLength(Particle& p, Track& t) const {

Framework/ProcessSequence/ProcessSequence.h

@@ -348,8 +348,8 @@ namespace corsika::process {
   OPSEQ(DecayProcess, ContinuousProcess)
   OPSEQ(DecayProcess, DecayProcess)
 
-  template <template <typename, typename> class T, typename A, typename B>
-  struct is_process_sequence<T<A, B> > {
+  template <typename A, typename B>
+    struct is_process_sequence<corsika::process::ProcessSequence<A, B> > {
     static const bool value = true;
   };
 

Processes/Sibyll/CMakeLists.txt

@@ -15,12 +15,18 @@ add_custom_command (
 set (
   MODEL_SOURCES
   ParticleConversion.cc
+  sibyll2.3c.f
+  sibyll2.3c.cc
+  gasdev.f
   )
 
 set (
   MODEL_HEADERS
   ParticleConversion.h
-  ProcessDecay.h
+  sibyll2.3c.h
+  SibStack.h
+  Decay.h
+  Interaction.h
   ${PROJECT_BINARY_DIR}/Processes/Sibyll/Generated.inc
   )
 
@@ -61,6 +67,7 @@ target_link_libraries (
   CORSIKAparticles
   CORSIKAunits
   CORSIKAthirdparty
+  CORSIKAgeometry
   )
 
 target_include_directories (

Processes/Sibyll/Decay.h

-#ifndef _include_ProcessDecay_h_
-#define _include_ProcessDecay_h_
+#ifndef _include_corsika_process_sibyll_decay_h_
+#define _include_corsika_process_sibyll_decay_h_
 
-#include <corsika/process/ContinuousProcess.h>
-
-#include <corsika/cascade/SibStack.h>
+#include <corsika/process/sibyll/SibStack.h>
 #include <corsika/process/sibyll/ParticleConversion.h>
+#include <corsika/process/DecayProcess.h>
 #include <corsika/setup/SetupTrajectory.h>
 
-// using namespace corsika::particles;
+#include <corsika/particles/ParticleProperties.h>
 
 namespace corsika::process {
 
@@ -56,11 +55,37 @@ namespace corsika::process {
       }
     }
 
-    class ProcessDecay : public corsika::process::BaseProcess<ProcessDecay> {
+    void setTrackedParticlesStable() {
+      /*
+	Sibyll is hadronic generator
+	only hadrons decay
+      */
+      // set particles unstable
+      setHadronsUnstable();
+      // make tracked particles stable
+      std::cout << "Interaction: setting tracked hadrons stable.." << std::endl;
+      setup::Stack ds;
+      ds.NewParticle().SetPID(particles::Code::PiPlus);
+      ds.NewParticle().SetPID(particles::Code::PiMinus);
+      ds.NewParticle().SetPID(particles::Code::KPlus);
+      ds.NewParticle().SetPID(particles::Code::KMinus);
+      ds.NewParticle().SetPID(particles::Code::K0Long);
+      ds.NewParticle().SetPID(particles::Code::K0Short);
+      
+      for (auto& p : ds) {
+	int s_id = process::sibyll::ConvertToSibyllRaw(p.GetPID());
+	// set particle stable by setting table value negative
+	s_csydec_.idb[s_id - 1] = (-1) * abs(s_csydec_.idb[s_id - 1]);
+	p.Delete();
+      }
+    }
+
+    class Decay : public corsika::process::DecayProcess<Decay> {
     public:
-      ProcessDecay() {}
+      Decay() {}
       void Init() {
-        // setHadronsUnstable();
+	setHadronsUnstable();
+        setTrackedParticlesStable();
       }
 
       void setAllStable() {
@@ -85,31 +110,32 @@ namespace corsika::process {
       friend void setHadronsUnstable();
 
       template <typename Particle>
-      double MinStepLength(Particle& p, setup::Trajectory&) const {
+      double GetLifetime(Particle& p) const {
         corsika::units::hep::EnergyType E = p.GetEnergy();
         corsika::units::hep::MassType m = corsika::particles::GetMass(p.GetPID());
-        // env.GetDensity();
 
-        const MassDensityType density = 1.25e-3 * kilogram / (1_cm * 1_cm * 1_cm);
+        //const MassDensityType density = 1.25e-3 * kilogram / (1_cm * 1_cm * 1_cm);
 
         const double gamma = E / m;
 
-        TimeType t0 = GetLifetime(p.GetPID());
-        cout << "ProcessDecay: code: " << (p.GetPID()) << endl;
-        cout << "ProcessDecay: MinStep: t0: " << t0 << endl;
-        cout << "ProcessDecay: MinStep: gamma: " << gamma << endl;
-        cout << "ProcessDecay: MinStep: density: " << density << endl;
+        const TimeType t0 = particles::GetLifetime(p.GetPID());
+        cout << "Decay: code: " << (p.GetPID()) << endl;
+        cout << "Decay: MinStep: t0: " << t0 << endl;
+        cout << "Decay: MinStep: gamma: " << gamma << endl;
+        // cout << "Decay: MinStep: density: " << density << endl;
         // return as column density
-        const double x0 = density * t0 * gamma * constants::c / kilogram * 1_cm * 1_cm;
-        cout << "ProcessDecay: MinStep: x0: " << x0 << endl;
-        int a = 1;
-        const double x = -x0 * log(s_rndm_(a));
-        cout << "ProcessDecay: next decay: " << x << endl;
-        return x;
+        // const double x0 = density * t0 * gamma * constants::c / kilogram * 1_cm * 1_cm;
+        // cout << "Decay: MinStep: x0: " << x0 << endl;
+	const double lifetime = gamma*t0 / 1_s;
+	cout << "Decay: MinStep: tau: " << lifetime << endl;
+        //int a = 1;
+        //const double x = -x0 * log(s_rndm_(a));
+        //cout << "Decay: next decay: " << x << endl;
+        return lifetime;
       }
 
       template <typename Particle, typename Stack>
-      void DoDiscrete(Particle& p, Stack& s) const {
+      void DoDecay(Particle& p, Stack& s) const {
         SibStack ss;
         ss.Clear();
         // copy particle to sibyll stack
@@ -122,7 +148,7 @@ namespace corsika::process {
         // set all particles/hadrons unstable
         setHadronsUnstable();
         // call sibyll decay
-        std::cout << "ProcessDecay: calling Sibyll decay routine.." << std::endl;
+        std::cout << "Decay: calling Sibyll decay routine.." << std::endl;
         decsib_();
         // print output
         int print_unit = 6;
@@ -144,11 +170,6 @@ namespace corsika::process {
         // empty sibyll stack
         ss.Clear();
       }
-
-      template <typename Particle, typename Stack>
-      EProcessReturn DoContinuous(Particle&, setup::Trajectory&, Stack&) const {
-        return EProcessReturn::eOk;
-      }
     };
   } // namespace sibyll
 } // namespace corsika::process

Processes/Sibyll/ParticleConversion.h

@@ -25,11 +25,11 @@ namespace corsika::process::sibyll {
 
 #include <corsika/process/sibyll/Generated.inc>
 
-  bool KnownBySibyll(corsika::particles::Code pCode) {
+  bool constexpr  KnownBySibyll(corsika::particles::Code pCode) {
     return isKnown[static_cast<corsika::particles::CodeIntType>(pCode)];
   }
 
-  bool CanInteract(corsika::particles::Code pCode) {
+  bool constexpr CanInteract(corsika::particles::Code pCode) {
     return canInteract[static_cast<corsika::particles::CodeIntType>(pCode)];
   }
 
@@ -43,12 +43,12 @@ namespace corsika::process::sibyll {
     return sibyll2corsika[static_cast<SibyllCodeIntType>(pCode) - minSibyll];
   }
 
-  int ConvertToSibyllRaw(corsika::particles::Code pCode) {
+  int constexpr ConvertToSibyllRaw(corsika::particles::Code pCode) {
     return (int)static_cast<corsika::process::sibyll::SibyllCodeIntType>(
         corsika::process::sibyll::ConvertToSibyll(pCode));
   }
 
-  int GetSibyllXSCode(corsika::particles::Code pCode) {
+  int constexpr GetSibyllXSCode(corsika::particles::Code pCode) {
     return corsika2sibyllXStype[static_cast<corsika::particles::CodeIntType>(pCode)];
   }
 

Processes/Sibyll/SibStack.h

 #ifndef _include_sibstack_h_
 #define _include_sibstack_h_
 
-#include <string>
-#include <vector>
-
-#include <corsika/cascade/sibyll2.3c.h>
+#include <corsika/process/sibyll/sibyll2.3c.h>
 #include <corsika/process/sibyll/ParticleConversion.h>
 #include <corsika/stack/Stack.h>
 #include <corsika/units/PhysicalUnits.h>
 
+#include <corsika/stack/super_stupid/SuperStupidStack.h>
+
 using namespace std;
 using namespace corsika::stack;
-using namespace corsika::units;
+using namespace corsika::units::si;
 using namespace corsika::geometry;
 
+namespace corsika::process::sibyll {
+  
 class SibStackData {
 
 public:
@@ -30,7 +31,7 @@ public:
   void SetMomentum(const int i, const super_stupid::MomentumVector& v) {
     auto tmp = v.GetComponents();
     for (int idx = 0; idx < 3; ++idx)
-      s_plist_.p[idx][i] = tmp[idx] / 1_GeV * si::constants::c;
+      s_plist_.p[idx][i] = tmp[idx] / 1_GeV * constants::c;
   }
 
   int GetId(const int i) const { return s_plist_.llist[i]; }
@@ -40,9 +41,9 @@ public:
   super_stupid::MomentumVector GetMomentum(const int i) const {
     CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
     corsika::geometry::QuantityVector<momentum_d> components{
-        s_plist_.p[0][i] * 1_GeV / si::constants::c,
-        s_plist_.p[1][i] * 1_GeV / si::constants::c,
-        s_plist_.p[2][i] * 1_GeV / si::constants::c};
+        s_plist_.p[0][i] * 1_GeV / constants::c,
+        s_plist_.p[1][i] * 1_GeV / constants::c,
+        s_plist_.p[2][i] * 1_GeV / constants::c};
     super_stupid::MomentumVector v1(rootCS, components);
     return v1;
   }
@@ -82,4 +83,6 @@ public:
 
 typedef Stack<SibStackData, ParticleInterface> SibStack;
 
+}
+
 #endif

Processes/Sibyll/code_generator.py

@@ -89,7 +89,7 @@ def generate_sibyll2corsika(pythia_db) :
             pDict[sib_code] = identifier
     
     nPart = max(pDict.keys()) - min(pDict.keys()) + 1
-    string += "std::array<corsika::particles::Code, {:d}> sibyll2corsika = {{\n".format(nPart)
+    string += "std::array<corsika::particles::Code, {:d}> constexpr sibyll2corsika = {{\n".format(nPart)
     
     for iPart in range(nPart) :
         if iPart in pDict: