diff --git a/Processes/Sibyll/ProcessDecay.h b/Processes/Sibyll/Decay.h
similarity index 100%
rename from Processes/Sibyll/ProcessDecay.h
rename to Processes/Sibyll/Decay.h
diff --git a/Processes/Sibyll/Interaction.h b/Processes/Sibyll/Interaction.h
new file mode 100644
index 0000000000000000000000000000000000000000..c0ceb87dd54c49660803df78725be7aa1d574a05
--- /dev/null
+++ b/Processes/Sibyll/Interaction.h
@@ -0,0 +1,283 @@
+#ifndef _corsika_process_sibyll_interaction_h_
+#define _corsika_process_sibyll_interaction_h_
+
+#include <corsika/process/InteractionProcess.h>
+
+//#include <corsika/setup/SetupStack.h>
+//#include <corsika/setup/SetupTrajectory.h>
+
+#include <corsika/process/sibyll/SibStack.h>
+#include <corsika/process/sibyll/sibyll2.3c.h>
+#include <corsika/process/sibyll/ParticleConversion.h>
+
+#include <corsika/particles/ParticleProperties.h>
+#include <corsika/units/PhysicalUnits.h>
+#include <corsika/random/RNGManager.h>
+
+using namespace corsika;
+using namespace corsika::process::sibyll;
+using namespace corsika::units::si;
+
+
+namespace corsika::process::sibyll {
+
+ // template <typename Stack, typename Track>
+ //template <typename Stack>
+ class Interaction : public corsika::process::InteractionProcess<Interaction> { // <Stack,Track>> {
+
+ //typedef typename Stack::ParticleType Particle;
+ //typedef typename corsika::setup::Stack::ParticleType Particle;
+ //typedef corsika::setup::Trajectory Track;
+
+ public:
+
+ Interaction() {}
+ ~Interaction() {}
+
+ void Init() {
+ corsika::random::RNGManager& rmng = corsika::random::RNGManager::GetInstance();
+ rmng.RegisterRandomStream("s_rndm");
+
+ // test random number generator
+ std::cout << "Interaction: "
+ << " 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_();
+ }
+
+ //void setTrackedParticlesStable();
+
+ template <typename Particle, typename Track>
+ double GetInteractionLength(Particle& p, Track&) const {
+
+ // coordinate system, get global frame of reference
+ CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
+
+ const particles::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() * constants::cSquared);
+ double Ecm = sqs / 1_GeV;
+
+ std::cout << "Interaction: "
+ << "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 << "Interaction: "
+ << "MinStep: sibyll return: " << prodCrossSection << std::endl;
+ CrossSectionType sig = prodCrossSection * 1_mbarn;
+ std::cout << "Interaction: "
+ << "MinStep: CrossSection (mb): " << sig / 1_mbarn << std::endl;
+
+ const MassType nucleon_mass = 0.93827_GeV / corsika::units::si::constants::cSquared;
+ std::cout << "Interaction: "
+ << "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 << "Interaction: "
+ << "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 << "Interaction: "
+ // << "next step (g/cm2): " << next_step << std::endl;
+ // return next_step;
+ }
+
+
+ template <typename Particle, typename Stack>
+ corsika::process::EProcessReturn 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 / constants::c, 0. * 1_GeV / constants::c,
+ 0. * 1_GeV / constants::c);
+ cout << "target momentum (GeV/c): "
+ << pTarget.GetComponents() / 1_GeV * constants::c << endl;
+ cout << "beam momentum (GeV/c): "
+ << p.GetMomentum().GetComponents() / 1_GeV * 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() *
+ 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 / constants::c);
+
+ std::cout << "Interaction: "
+ << " DoDiscrete: gamma:" << gamma << endl;
+ std::cout << "Interaction: "
+ << " DoDiscrete: gambet:" << gambet.GetComponents() << endl;
+
+ int kBeam = process::sibyll::ConvertToSibyllRaw(p.GetPID());
+
+ std::cout << "Interaction: "
+ << " DoInteraction: E(GeV):" << E / 1_GeV << " Ecm(GeV): " << Ecm / 1_GeV
+ << std::endl;
+ if (E < 8.5_GeV || Ecm < 10_GeV) {
+ std::cout << "Interaction: "
+ << " DoInteraction: dropping particle.." << std::endl;
+ p.Delete();
+ } 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] * constants::c) /
+ (gamma + 1.);
+ pnorm += psib.GetEnergy();
+
+ for (int j = 0; j < 3; ++j) {
+ p_lab_components[j] = pSibyllComponents[j] - (-1) * pnorm *
+ gammaBetaComponents[j] /
+ constants::c;
+ en_lab -=
+ (-1) * pSibyllComponents[j] * gammaBetaComponents[j] * 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 *
+ // constants::c << endl;
+ }
+ }
+ return process::EProcessReturn::eOk;
+ }
+
+ };
+
+}
+
+#endif
diff --git a/Framework/Cascade/SibStack.h b/Processes/Sibyll/SibStack.h
similarity index 100%
rename from Framework/Cascade/SibStack.h
rename to Processes/Sibyll/SibStack.h
diff --git a/Framework/Cascade/gasdev.f b/Processes/Sibyll/gasdev.f
similarity index 100%
rename from Framework/Cascade/gasdev.f
rename to Processes/Sibyll/gasdev.f
diff --git a/Framework/Cascade/rndm_dbl.f b/Processes/Sibyll/rndm_dbl.f
similarity index 100%
rename from Framework/Cascade/rndm_dbl.f
rename to Processes/Sibyll/rndm_dbl.f
diff --git a/Processes/Sibyll/sibyll2.3c.cc b/Processes/Sibyll/sibyll2.3c.cc
new file mode 100644
index 0000000000000000000000000000000000000000..336617fbfe42e31dcc00b9d8215e8a6f5e24c485
--- /dev/null
+++ b/Processes/Sibyll/sibyll2.3c.cc
@@ -0,0 +1,11 @@
+#include <corsika/process/sibyll/sibyll2.3c.h>
+
+#include <corsika/random/RNGManager.h>
+
+
+double s_rndm_(int&) {
+ static corsika::random::RNG& rmng =
+ corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");
+ ;
+ return rmng() / (double)rmng.max();
+}
diff --git a/Framework/Cascade/sibyll2.3c.f b/Processes/Sibyll/sibyll2.3c.f
similarity index 100%
rename from Framework/Cascade/sibyll2.3c.f
rename to Processes/Sibyll/sibyll2.3c.f
diff --git a/Framework/Cascade/sibyll2.3c.h b/Processes/Sibyll/sibyll2.3c.h
similarity index 100%
rename from Framework/Cascade/sibyll2.3c.h
rename to Processes/Sibyll/sibyll2.3c.h