diff --git a/CMakeLists.txt b/CMakeLists.txt
index 7eb52c6d7c60a351adabb860dab45008f8ba4180..d5f78faea60daf8fa2eb980a0e3ac64eb6d12db7 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -72,6 +72,8 @@ if (Boost_FOUND)
   include_directories(${Boost_INCLUDE_DIRS})
 endif (Boost_FOUND)
 
+find_package (Pythia8) # optional
+  
 find_package (Eigen3 REQUIRED)
 #find_package (HDF5) # not yet needed
 
diff --git a/CMakeModules/CorsikaUtilities.cmake b/CMakeModules/CorsikaUtilities.cmake
index 85105468ad429e013f499fd221b2d657e60fd612..9b00787289f634872b707c208cc899136418d597 100644
--- a/CMakeModules/CorsikaUtilities.cmake
+++ b/CMakeModules/CorsikaUtilities.cmake
@@ -1,4 +1,14 @@
+#
+# (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.
+#
 
+#################################################
 #
 # takes a list of input files and prepends a path
 #
@@ -12,6 +22,7 @@ function (CORSIKA_PREPEND_PATH return prefix)
 endfunction (CORSIKA_PREPEND_PATH)
 
 
+#################################################
 #
 # use: CORSIKA_COPY_HEADERS_TO_NAMESPACE theLib theNamesapce header1.h header2.h ...
 #
@@ -58,6 +69,7 @@ endfunction (CORSIKA_COPY_HEADERS_TO_NAMESPACE)
 
 
 
+#################################################
 #
 # use: CORSIKA_ADD_FILES_ABSOLUTE varname
 #
@@ -81,6 +93,7 @@ endmacro(CORSIKA_ADD_FILES_ABSOLUTE)
 
 
 
+#################################################
 #
 # central macro to activate unit tests in cmake
 #
diff --git a/CMakeModules/FindPythia8.cmake b/CMakeModules/FindPythia8.cmake
new file mode 100644
index 0000000000000000000000000000000000000000..57aea3d0bb86f4a90a17e36cb0e18da5561ba309
--- /dev/null
+++ b/CMakeModules/FindPythia8.cmake
@@ -0,0 +1,73 @@
+#
+# (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.
+#
+
+#################################################
+#
+# run pythia8-config and interpret result
+#
+
+function (_PYTHIA8_CONFIG_ option variable type doc)
+  execute_process(COMMAND ${PYTHIA8_CONFIG} ${option}
+    OUTPUT_VARIABLE _local_out_
+    RESULT_VARIABLE _local_res_)
+  string(REGEX REPLACE "\n$" "" _local_out_ "${_local_out_}")
+  if (NOT ${_local_res_} EQUAL 0)
+    message ("Error in running ${PYTHIA8_CONFIG} ${option}")
+  else ()
+    set(${variable} "${_local_out_}" CACHE ${type} ${doc})
+  endif ()
+endfunction (_PYTHIA8_CONFIG_)
+  
+
+
+#################################################
+# 
+# Searched PYTHIA8 on system. Expect pythia8-config in PATH, or typical installation location
+#
+# This module defines
+# HAVE_PYTHIA8
+# PYTHIA8_INCLUDE_DIR   where to locate Pythia.h file
+# PYTHIA8_LIBRARY       where to find the libpythia8 library
+# PYTHIA8_LIBRARIES     (not cached) the libraries to link against to use Pythia8
+# PYTHIA8_VERSION       version of Pythia8 if found
+#
+
+set(_SEARCH_PYTHIA8_
+  ${PROJECT_BINARY_DIR}/ThirdParty/pythia8-install
+  ${PYTHIA8}
+  $ENV{PYTHIA8}
+  ${PYTHIA8DIR}
+  $ENV{PYTHIA8DIR}
+  ${PYTHIA8_DIR}
+  $ENV{PYTHIA8_DIR}
+  ${PYTHIA8_ROOT}
+  $ENV{PYTHIA8_ROOT}
+  /opt/pythia8)
+
+find_program(PYTHIA8_CONFIG
+  NAME pythia8-config
+  PATHS ${_SEARCH_PYTHIA8_}
+  PATH_SUFFIXES "/bin"
+  DOC "The location of the pythia8-config script")
+
+if (PYTHIA8_CONFIG)
+  set(HAVE_PYTHIA8 1 CACHE BOOL "presence of pythia8, found via pythia8-config")
+
+  _PYTHIA8_CONFIG_("--prefix" PYTHIA8_PREFIX PATH "location of pythia8 installation")
+  _PYTHIA8_CONFIG_("--includedir" PYTHIA8_INCLUDE_DIR PATH "pythia8 include directory")
+  _PYTHIA8_CONFIG_("--libs" PYTHIA8_LIBRARY STRING "the pythia8 libs")
+  _PYTHIA8_CONFIG_("--datadir" PYTHIA8_DATA_DIR PATH "the pythia8 data dir")
+  _PYTHIA8_CONFIG_("--cxxflags" PYTHIA8_CXXFLAGS STRING "the pythia8 cxxflags")
+endif ()
+
+# standard cmake infrastructure:
+include(FindPackageHandleStandardArgs)
+find_package_handle_standard_args(Pythia8 DEFAULT_MSG PYTHIA8_PREFIX PYTHIA8_INCLUDE_DIR PYTHIA8_LIBRARY)
+mark_as_advanced(PYTHIA8_DATA_DIR PYTHIA8_CXXFLAGS PYTHIA8_INCLUDE_DIR PYTHIA8_LIBRARY)
diff --git a/Documentation/Examples/CMakeLists.txt b/Documentation/Examples/CMakeLists.txt
index 621cab01d1facaa687c12f287c12ed2d54a6501f..03eacd0a0206365680151f67d1131d7362e3d31c 100644
--- a/Documentation/Examples/CMakeLists.txt
+++ b/Documentation/Examples/CMakeLists.txt
@@ -18,14 +18,17 @@ CORSIKA_ADD_TEST (logger_example)
 
 add_executable (stack_example stack_example.cc)
 target_compile_options(stack_example PRIVATE -g) # do not skip asserts
-target_link_libraries (stack_example SuperStupidStack CORSIKAunits CORSIKAlogging)
+target_link_libraries (stack_example SuperStupidStack CORSIKAunits
+  CORSIKAlogging)
 CORSIKA_ADD_TEST (stack_example)
 
 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
-   ProcessSibyll
+target_link_libraries (cascade_example SuperStupidStack CORSIKAunits
+  CORSIKAlogging
+  CORSIKArandom
+  ProcessSibyll
+  ProcessPythia
   CORSIKAcascade
   ProcessStackInspector
   ProcessTrackWriter
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index f9d1e95a83aa38ed0475c06850698d2b867f6f99..62e2347f4f3f096f7314cd9087ffcef6ff0fdf40 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -28,6 +28,8 @@
 #include <corsika/process/sibyll/Interaction.h>
 #include <corsika/process/sibyll/NuclearInteraction.h>
 
+#include <corsika/process/pythia/Decay.h>
+
 #include <corsika/process/track_writer/TrackWriter.h>
 
 #include <corsika/units/PhysicalUnits.h>
@@ -251,10 +253,17 @@ int main() {
   tracking_line::TrackingLine<setup::Stack, setup::Trajectory> tracking(env);
   stack_inspector::StackInspector<setup::Stack> p0(true);
 
+  const std::vector<particles::Code> trackedHadrons = {
+        particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
+        particles::Code::KMinus, particles::Code::K0Long,  particles::Code::K0Short};
+
+  
   random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
   process::sibyll::Interaction sibyll(env);
   process::sibyll::NuclearInteraction sibyllNuc(env, sibyll);
-  process::sibyll::Decay decay;
+  process::sibyll::Decay decay(trackedHadrons);
+  //random::RNGManager::GetInstance().RegisterRandomStream("pythia");
+  //process::pythia::Decay decay(trackedHadrons);
   ProcessCut cut(20_GeV);
 
   // random::RNGManager::GetInstance().RegisterRandomStream("HadronicElasticModel");
@@ -274,7 +283,7 @@ int main() {
   setup::Stack stack;
   stack.Clear();
   const Code beamCode = Code::Nucleus;
-  const int nuclA = 56;
+  const int nuclA = 4;
   const int nuclZ = int(nuclA / 2.15 + 0.7);
   const HEPMassType mass = particles::Proton::GetMass() * nuclZ +
                            (nuclA - nuclZ) * particles::Neutron::GetMass();
diff --git a/Framework/Geometry/Trajectory.h b/Framework/Geometry/Trajectory.h
index 32063061058c446430cd1ba52e13924c20685eff..e93a9d3029121c8183fd0e416e53d47194050eeb 100644
--- a/Framework/Geometry/Trajectory.h
+++ b/Framework/Geometry/Trajectory.h
@@ -47,10 +47,10 @@ namespace corsika::geometry {
     void LimitEndTo(corsika::units::si::LengthType limit) {
       fTimeLength = T::TimeFromArclength(limit);
     }
-    
+
     auto NormalizedDirection() const {
-        static_assert(std::is_same_v<T, corsika::geometry::Line>);
-        return T::GetV0().normalized();
+      static_assert(std::is_same_v<T, corsika::geometry::Line>);
+      return T::GetV0().normalized();
     }
   };
 
diff --git a/Framework/Geometry/testGeometry.cc b/Framework/Geometry/testGeometry.cc
index d73ae99b2f21590984aae1e90c683777549b7b32..dae272be1bc052b3174ad42ae0b2a3f50777322b 100644
--- a/Framework/Geometry/testGeometry.cc
+++ b/Framework/Geometry/testGeometry.cc
@@ -177,8 +177,10 @@ TEST_CASE("Trajectories") {
     CHECK(line.GetPosition(t).GetCoordinates() == base.GetPosition(1.).GetCoordinates());
 
     CHECK(base.ArcLength(1_s, 2_s) / 1_m == Approx(3));
-    
-    CHECK((base.NormalizedDirection().GetComponents(rootCS) - QuantityVector<dimensionless_d>{1, 0, 0}).eVector.norm() == Approx(0).margin(absMargin));
+
+    CHECK((base.NormalizedDirection().GetComponents(rootCS) -
+           QuantityVector<dimensionless_d>{1, 0, 0})
+              .eVector.norm() == Approx(0).margin(absMargin));
   }
 
   SECTION("Helix") {
diff --git a/Processes/CMakeLists.txt b/Processes/CMakeLists.txt
index 26a44fe6e1ea2ddf806a4f9fd1e2b0552153a64d..eb78018365128dde2308ba1730ed2920a20ae800 100644
--- a/Processes/CMakeLists.txt
+++ b/Processes/CMakeLists.txt
@@ -1,6 +1,9 @@
 
 add_subdirectory (NullModel)
 add_subdirectory (Sibyll)
+if (PYTHIA8_FOUND)
+  add_subdirectory (Pythia)
+endif (PYTHIA8_FOUND)
 add_subdirectory (StackInspector)
 add_subdirectory (TrackWriter)
 add_subdirectory (TrackingLine)
@@ -10,5 +13,8 @@ add_subdirectory (HadronicElasticModel)
 add_library (CORSIKAprocesses INTERFACE)
 add_dependencies(CORSIKAprocesses ProcessNullModel)
 add_dependencies(CORSIKAprocesses ProcessSibyll)
+if (PYTHIA8_FOUND)
+  add_dependencies(CORSIKAprocesses ProcessPythia)
+endif (PYTHIA8_FOUND)
 add_dependencies(CORSIKAprocesses ProcessStackInspector)
 add_dependencies(CORSIKAprocesses ProcessTrackingLine)
diff --git a/Processes/Pythia/CMakeLists.txt b/Processes/Pythia/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..3af574c892dd423cdef3b7f018265c9f6f3912af
--- /dev/null
+++ b/Processes/Pythia/CMakeLists.txt
@@ -0,0 +1,91 @@
+set(Python_ADDITIONAL_VERSIONS 3)
+find_package(PythonInterp 3 REQUIRED)
+  
+set (
+  MODEL_SOURCES
+  Decay.cc
+  Random.cc
+  )
+
+set (
+  MODEL_HEADERS
+  Decay.h
+  Random.h
+  )
+
+set (
+  MODEL_NAMESPACE
+  corsika/process/pythia
+  )
+
+add_library (ProcessPythia STATIC ${MODEL_SOURCES})
+CORSIKA_COPY_HEADERS_TO_NAMESPACE (ProcessPythia ${MODEL_NAMESPACE} ${MODEL_HEADERS})
+
+
+set_target_properties (
+  ProcessPythia
+  PROPERTIES
+  VERSION ${PROJECT_VERSION}
+  SOVERSION 1
+  PUBLIC_HEADER "${MODEL_HEADERS}"
+  )
+
+# target dependencies on other libraries (also the header onlys)
+target_link_libraries (
+  ProcessPythia
+  CORSIKAparticles
+  CORSIKAutilities
+  CORSIKAunits
+  CORSIKAthirdparty
+  CORSIKAgeometry
+  CORSIKAenvironment
+  ${PYTHIA8_LIBRARY}
+  )
+
+target_include_directories (
+  ProcessPythia
+  INTERFACE
+  $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
+  $<INSTALL_INTERFACE:include/include>
+  )
+# do it again, since we want -isystem for pythia exteranls
+target_include_directories (
+  ProcessPythia
+  SYSTEM PUBLIC
+  ${PYTHIA8_INCLUDE_DIR}
+  )
+
+target_include_directories (
+  ProcessPythia
+  PRIVATE
+  ${PYTHIA8_INCLUDE_DIR}
+  )
+
+install (
+  TARGETS ProcessPythia
+  LIBRARY DESTINATION lib
+  ARCHIVE DESTINATION lib
+#  PUBLIC_HEADER DESTINATION include/${MODEL_NAMESPACE}
+  )
+
+
+# --------------------
+# code unit testing
+add_executable (testPythia
+  testPythia.cc
+  ${MODEL_HEADERS}
+  )
+
+  
+target_link_libraries (
+  testPythia
+  ProcessPythia
+  CORSIKAsetup
+  CORSIKArandom
+  CORSIKAgeometry
+  CORSIKAunits
+  CORSIKAthirdparty # for catch2
+  ${PYTHIA8_LIBRARY}
+  )
+
+CORSIKA_ADD_TEST(testPythia)
diff --git a/Processes/Pythia/Decay.cc b/Processes/Pythia/Decay.cc
new file mode 100644
index 0000000000000000000000000000000000000000..8bdfb554fb44538c816e76440acb8e0cf0845609
--- /dev/null
+++ b/Processes/Pythia/Decay.cc
@@ -0,0 +1,156 @@
+/*
+ * (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/process/pythia/Decay.h>
+#include <corsika/process/pythia/Random.h>
+#include <Pythia8/Pythia.h>
+
+#include <corsika/setup/SetupStack.h>
+#include <corsika/setup/SetupTrajectory.h>
+
+using std::cout;
+using std::endl;
+using std::tuple;
+using std::vector;
+
+using namespace corsika;
+using namespace corsika::setup;
+using Particle = Stack::StackIterator; // ParticleType;
+using Track = Trajectory;
+
+namespace corsika::process::pythia {
+
+  Decay::Decay(vector<particles::Code> pParticles)
+	: fTrackedParticles(pParticles) {}
+  
+  Decay::~Decay() { cout << "Pythia::Decay n=" << fCount << endl; }
+  
+  void Decay::Init() {
+    
+    Decay::SetParticleListStable(fTrackedParticles);
+
+    // set random number generator in pythia
+    Pythia8::RndmEngine* rndm = new corsika::process::pythia::Random();
+    fPythia.setRndmEnginePtr( rndm );
+
+    fPythia.readString("Next:numberShowInfo = 0");
+    fPythia.readString("Next:numberShowProcess = 0");
+    fPythia.readString("Next:numberShowEvent = 0");
+
+    fPythia.readString("Print:quiet = on");
+    
+    fPythia.readString("ProcessLevel:all = off");
+    fPythia.readString("ProcessLevel:resonanceDecays = off");
+
+    fPythia.particleData.readString("59:m0 = 101.00");
+    
+    fPythia.init();
+  }
+
+  void Decay::SetParticleListStable(const vector<particles::Code> particleList) {
+    for (auto p : particleList)
+      Decay::SetStable( p );
+  }
+
+  void Decay::SetUnstable(const particles::Code pCode) {
+    cout << "Pythia::Decay: setting " << pCode << " unstable.." << endl;
+    fPythia.particleData.mayDecay( static_cast<int>( particles::GetPDG(pCode) ) , true);
+  }
+
+  void Decay::SetStable(const particles::Code pCode) {
+    cout << "Pythia::Decay: setting " << pCode << " stable.." << endl;
+    fPythia.particleData.mayDecay( static_cast<int>( particles::GetPDG(pCode) ) , false);
+  }
+
+  template <>
+  units::si::TimeType Decay::GetLifetime(Particle const& p) {
+    using namespace units::si;
+
+    HEPEnergyType E = p.GetEnergy();
+    HEPMassType m = particles::GetMass(p.GetPID());
+
+    const double gamma = E / m;
+
+    const TimeType t0 = particles::GetLifetime(p.GetPID());
+    auto const lifetime = gamma * t0;
+
+    return lifetime;
+  }
+
+  template <>
+  void Decay::DoDecay(Particle& p, Stack&) {
+    using geometry::Point;
+    using namespace units;
+    using namespace units::si;
+
+    auto const decayPoint = p.GetPosition();
+    auto const t0 = p.GetTime();
+
+    // coordinate system, get global frame of reference
+    geometry::CoordinateSystem& rootCS =
+      geometry::RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
+    
+    fCount++;
+
+    // pythia stack
+    Pythia8::Event& event = fPythia.event;
+    event.reset();
+    
+    // set particle unstable
+    Decay::SetUnstable( p.GetPID() );
+
+    // input particle PDG
+    auto const pdgCode = static_cast<int>( particles::GetPDG( p.GetPID() ));
+
+    auto const pcomp =     p.GetMomentum().GetComponents();
+    double px = pcomp[0] / 1_GeV ;
+    double py = pcomp[1] / 1_GeV ;
+    double pz = pcomp[2] / 1_GeV ;
+    double en = p.GetEnergy() / 1_GeV;
+    double m = particles::GetMass( p.GetPID() ) / 1_GeV;
+
+    // add particle to pythia stack
+    event.append( pdgCode, 1, 0, 0, px, py, pz, en, m);
+
+    if (!fPythia.next())
+      cout << "Pythia::Decay: decay failed!" << endl;
+    else
+      cout << "Pythia::Decay: particles after decay: " << event.size() << endl;
+
+    // list final state
+    event.list();
+
+    // loop over final state
+    for (int i = 0; i < event.size(); ++i)
+      if (event[i].isFinal()) {
+	auto const pyId = particles::ConvertFromPDG(static_cast<particles::PDGCode>(event[i].id()));
+	HEPEnergyType pyEn = event[i].e() * 1_GeV;
+	MomentumVector pyP(rootCS, { event[i].px() * 1_GeV,
+				     event[i].py() * 1_GeV,
+				     event[i].pz() * 1_GeV});
+	
+	cout << "particle: id=" << pyId << " momentum=" << pyP.GetComponents() / 1_GeV
+	     << " energy=" << pyEn << endl;
+	
+	p.AddSecondary( tuple<particles::Code, units::si::HEPEnergyType,
+			corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
+			  pyId, pyEn, pyP, decayPoint, t0} );
+      }
+    
+    // set particle stable
+    Decay::SetStable( p.GetPID() );
+
+    // remove original particle from corsika stack
+    p.Delete();
+    //    if (fCount>10) throw std::runtime_error("stop here");
+  }
+
+  
+} // namespace corsika::process::pythia
diff --git a/Processes/Pythia/Decay.h b/Processes/Pythia/Decay.h
new file mode 100644
index 0000000000000000000000000000000000000000..06f1838867a9be760128eefadaf3e653e8c89598
--- /dev/null
+++ b/Processes/Pythia/Decay.h
@@ -0,0 +1,51 @@
+
+/*
+ * (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.
+ */
+
+#ifndef _include_corsika_process_pythia_decay_h_
+#define _include_corsika_process_pythia_decay_h_
+
+#include <corsika/particles/ParticleProperties.h>
+#include <corsika/process/DecayProcess.h>
+#include <Pythia8/Pythia.h>
+
+namespace corsika::process {
+
+  namespace pythia {
+
+    typedef corsika::geometry::Vector<corsika::units::si::hepmomentum_d> MomentumVector;
+    
+    class Decay : public corsika::process::DecayProcess<Decay> {
+      const std::vector<particles::Code> fTrackedParticles;
+      int fCount = 0;
+
+    public:
+      Decay(std::vector<corsika::particles::Code> );
+      ~Decay();
+      void Init();
+
+      void SetParticleListStable(const std::vector<particles::Code>);
+      void SetUnstable(const corsika::particles::Code );
+      void SetStable(const corsika::particles::Code );
+
+      template <typename Particle>
+      corsika::units::si::TimeType GetLifetime(Particle const& p);
+
+      template <typename Particle, typename Stack>
+      void DoDecay(Particle& p, Stack&);
+
+    private:
+      Pythia8::Pythia fPythia;
+    };
+    
+  } // namespace pythia
+} // namespace corsika::process
+
+#endif
diff --git a/Processes/Pythia/Random.cc b/Processes/Pythia/Random.cc
new file mode 100644
index 0000000000000000000000000000000000000000..1432cc3cb4267a641118baf8c22e1e37f45ba051
--- /dev/null
+++ b/Processes/Pythia/Random.cc
@@ -0,0 +1,21 @@
+/*
+ * (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/process/pythia/Random.h>
+
+namespace corsika::process::pythia {
+
+  double Random::flat(){
+    std::uniform_real_distribution<double> dist;
+    return dist(fRNG);
+  }
+
+  
+} // namespace corsika::process::pythia
diff --git a/Processes/Pythia/Random.h b/Processes/Pythia/Random.h
new file mode 100644
index 0000000000000000000000000000000000000000..885bb598a983079b9b0fe438a78762f109a88067
--- /dev/null
+++ b/Processes/Pythia/Random.h
@@ -0,0 +1,32 @@
+
+/*
+ * (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.
+ */
+
+#ifndef _include_corsika_process_pythia_random_h_
+#define _include_corsika_process_pythia_random_h_
+
+#include <corsika/random/RNGManager.h>
+#include <Pythia8/Pythia.h>
+
+namespace corsika::process {
+
+  namespace pythia {
+
+    class Random : public Pythia8::RndmEngine {
+      double flat();
+    private:
+      corsika::random::RNG& fRNG =
+        corsika::random::RNGManager::GetInstance().GetRandomStream("pythia");
+    };
+    
+  } // namespace pythia
+} // namespace corsika::process
+
+#endif
diff --git a/Processes/Pythia/testPythia.cc b/Processes/Pythia/testPythia.cc
new file mode 100644
index 0000000000000000000000000000000000000000..05895edaca71bc5a54ad2dccebd19fc7cdd883dc
--- /dev/null
+++ b/Processes/Pythia/testPythia.cc
@@ -0,0 +1,157 @@
+
+/*
+ * (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/process/pythia/Decay.h>
+#include <Pythia8/Pythia.h>
+
+#include <corsika/random/RNGManager.h>
+
+#include <corsika/particles/ParticleProperties.h>
+
+#include <corsika/geometry/Point.h>
+#include <corsika/units/PhysicalUnits.h>
+
+#define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one
+                          // cpp file
+#include <catch2/catch.hpp>
+
+TEST_CASE("Pythia", "[processes]") {
+
+  SECTION("linking pythia") {
+    using namespace Pythia8;
+    using std::cout;
+    using std::endl;
+
+    // Generator. Process selection. LHC initialization. Histogram.
+    Pythia pythia;
+
+    pythia.readString("Next:numberShowInfo = 0");
+    pythia.readString("Next:numberShowProcess = 0");
+    pythia.readString("Next:numberShowEvent = 0");
+
+    pythia.readString("ProcessLevel:all = off");
+
+    pythia.init();
+
+    Event& event = pythia.event;
+    event.reset();
+
+    pythia.particleData.mayDecay(321, true);
+    double pz = 100.;
+    double m = 0.49368;
+    event.append(321, 1, 0, 0, 0., 0., 100., sqrt(pz * pz + m * m), m);
+
+    if (!pythia.next())
+      cout << "decay failed!" << endl;
+    else
+      cout << "particles after decay: " << event.size() << endl;
+    event.list();
+
+    // loop over final state
+    for (int i = 0; i < pythia.event.size(); ++i)
+      if (pythia.event[i].isFinal()) {
+	cout << "particle: id=" << pythia.event[i].id() << endl;
+      }
+  }
+
+
+  SECTION("pythia interface") {
+    using namespace corsika;
+    
+    const std::vector<particles::Code> particleList = {
+        particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
+        particles::Code::KMinus, particles::Code::K0Long,  particles::Code::K0Short};
+
+    random::RNGManager::GetInstance().RegisterRandomStream("pythia");
+    
+    process::pythia::Decay model(particleList);
+    
+    model.Init();
+
+  }
+
+}
+
+#include <corsika/geometry/Point.h>
+#include <corsika/geometry/RootCoordinateSystem.h>
+#include <corsika/geometry/Vector.h>
+
+#include <corsika/units/PhysicalUnits.h>
+
+#include <corsika/particles/ParticleProperties.h>
+#include <corsika/setup/SetupStack.h>
+#include <corsika/setup/SetupTrajectory.h>
+
+#include <corsika/environment/Environment.h>
+#include <corsika/environment/HomogeneousMedium.h>
+#include <corsika/environment/NuclearComposition.h>
+
+using namespace corsika;
+using namespace corsika::units::si;
+
+TEST_CASE("pytia decay"){  
+
+  // setup environment, geometry
+  environment::Environment env;
+  auto& universe = *(env.GetUniverse());
+
+  auto theMedium = environment::Environment::CreateNode<geometry::Sphere>(
+      geometry::Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
+      1_km * std::numeric_limits<double>::infinity());
+
+  using MyHomogeneousModel = environment::HomogeneousMedium<environment::IMediumModel>;
+  theMedium->SetModelProperties<MyHomogeneousModel>(
+      1_kg / (1_m * 1_m * 1_m),
+      environment::NuclearComposition(
+          std::vector<particles::Code>{particles::Code::Oxygen}, std::vector<float>{1.}));
+
+  universe.AddChild(std::move(theMedium));
+
+  const geometry::CoordinateSystem& cs = env.GetCoordinateSystem();
+
+  geometry::Point const origin(cs, {0_m, 0_m, 0_m});
+  geometry::Vector<units::si::SpeedType::dimension_type> v(cs, 0_m / second, 0_m / second,
+                                                           1_m / second);
+  geometry::Line line(origin, v);
+  geometry::Trajectory<geometry::Line> track(line, 10_s);
+
+  random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
+
+  SECTION("pythia decay") {
+    
+    setup::Stack stack;
+    const HEPEnergyType E0 = 10_GeV;
+    HEPMomentumType P0 =
+        sqrt(E0 * E0 - particles::PiPlus::GetMass() * particles::PiPlus::GetMass());
+    auto plab = corsika::stack::MomentumVector(cs, {0_GeV, 0_GeV, -P0});
+    geometry::Point pos(cs, 0_m, 0_m, 0_m);
+    auto particle = stack.AddParticle(
+        std::tuple<particles::Code, units::si::HEPEnergyType,
+                   corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
+            particles::Code::PiPlus, E0, plab, pos, 0_ns});
+
+    
+    const std::vector<particles::Code> particleList = {
+        particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
+        particles::Code::KMinus, particles::Code::K0Long,  particles::Code::K0Short};
+
+    random::RNGManager::GetInstance().RegisterRandomStream("pythia");
+    
+    process::pythia::Decay model(particleList);
+    
+    model.Init();
+    /*[[maybe_unused]] const process::EProcessReturn ret =*/model.DoDecay(particle,
+                                                                          stack);
+    [[maybe_unused]] const TimeType time = model.GetLifetime(particle);
+    
+  }
+
+}
diff --git a/Processes/Sibyll/Decay.cc b/Processes/Sibyll/Decay.cc
index 1e0d1ca2b64be8db9444c94894fc6983a01e62ba..470ea1819c70b0cd8cfad5b5ddb38b9f4b2f8b0f 100644
--- a/Processes/Sibyll/Decay.cc
+++ b/Processes/Sibyll/Decay.cc
@@ -29,44 +29,37 @@ using Track = Trajectory;
 
 namespace corsika::process::sibyll {
 
-  Decay::Decay() {}
+  Decay::Decay(vector<particles::Code>pParticles)
+    : fTrackedParticles(pParticles) {}
   Decay::~Decay() { cout << "Sibyll::Decay n=" << fCount << endl; }
   void Decay::Init() {
-    setHadronsUnstable();
-    setTrackedParticlesStable();
+    SetHadronsUnstable();
+    SetParticleListStable( fTrackedParticles );
   }
 
-  void Decay::setTrackedParticlesStable() {
+  void Decay::SetParticleListStable(const vector<particles::Code> particleList) {
     /*
        Sibyll is hadronic generator
        only hadrons decay
      */
     // set particles unstable
-    setHadronsUnstable();
-    // make tracked particles stable
+    SetHadronsUnstable();
     cout << "Interaction: setting tracked hadrons stable.." << endl;
-    const vector<particles::Code> particleList = {
-        particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
-        particles::Code::KMinus, particles::Code::K0Long,  particles::Code::K0Short};
-
-    for (auto p : particleList) {
-      // set particle stable by setting table value negative
-      const int sibid = process::sibyll::ConvertToSibyllRaw(p);
-      s_csydec_.idb[sibid - 1] = (-1) * abs(s_csydec_.idb[sibid - 1]);
-    }
+    for (auto p : particleList) 
+      Decay::SetStable( p );
   }
 
-  void Decay::setUnstable(const particles::Code pCode) {
+  void Decay::SetUnstable(const particles::Code pCode) {
     int s_id = process::sibyll::ConvertToSibyllRaw(pCode);
     s_csydec_.idb[s_id - 1] = abs(s_csydec_.idb[s_id - 1]);
   }
 
-  void Decay::setStable(const particles::Code pCode) {
+  void Decay::SetStable(const particles::Code pCode) {
     int s_id = process::sibyll::ConvertToSibyllRaw(pCode);
     s_csydec_.idb[s_id - 1] = (-1) * abs(s_csydec_.idb[s_id - 1]);
   }
 
-  void Decay::setAllStable() {
+  void Decay::SetAllStable() {
     // name? also makes EM particles stable
 
     cout << "Decay: setting all particles stable.." << endl;
@@ -83,7 +76,7 @@ namespace corsika::process::sibyll {
     }
   }
 
-  void Decay::setHadronsUnstable() {
+  void Decay::SetHadronsUnstable() {
 
     // name? also makes EM particles stable
 
@@ -157,16 +150,14 @@ namespace corsika::process::sibyll {
     // remember position
     Point const decayPoint = p.GetPosition();
     TimeType const t0 = p.GetTime();
-    // remove original particle from corsika stack
-    p.Delete();
     // set all particles/hadrons unstable
     // setHadronsUnstable();
-    setUnstable(pCode);
+    SetUnstable(pCode);
     // call sibyll decay
     cout << "Decay: calling Sibyll decay routine.." << endl;
     decsib_();
     // reset to stable
-    setStable(pCode);
+    SetStable(pCode);
     // print output
     int print_unit = 6;
     sib_list_(print_unit);
@@ -184,6 +175,9 @@ namespace corsika::process::sibyll {
     }
     // empty sibyll stack
     ss.Clear();
+    // remove original particle from corsika stack
+    p.Delete();
+
   }
 
 } // namespace corsika::process::sibyll
diff --git a/Processes/Sibyll/Decay.h b/Processes/Sibyll/Decay.h
index c078af7b8763c4d6df9015177ed764f6a663b460..5a1f6323b0cb0dde9d73395d9ecc77f460860c53 100644
--- a/Processes/Sibyll/Decay.h
+++ b/Processes/Sibyll/Decay.h
@@ -20,18 +20,19 @@ namespace corsika::process {
   namespace sibyll {
 
     class Decay : public corsika::process::DecayProcess<Decay> {
+      std::vector<particles::Code> fTrackedParticles;
       int fCount = 0;
 
     public:
-      Decay();
+      Decay(std::vector<particles::Code>);
       ~Decay();
       void Init();
 
-      void setTrackedParticlesStable();
-      void setUnstable(const corsika::particles::Code pCode);
-      void setStable(const corsika::particles::Code pCode);
-      void setAllStable();
-      void setHadronsUnstable();
+      void SetParticleListStable(const std::vector<particles::Code>);
+      void SetUnstable(const corsika::particles::Code );
+      void SetStable(const corsika::particles::Code );
+      void SetAllStable();
+      void SetHadronsUnstable();
 
       template <typename Particle>
       corsika::units::si::TimeType GetLifetime(Particle const& p);
diff --git a/Processes/Sibyll/Interaction.cc b/Processes/Sibyll/Interaction.cc
index 0e65c074bba270dcd4db98866e079893d80bf622..5594ec39de31270b0e1014f5aa1e4b64c917c99f 100644
--- a/Processes/Sibyll/Interaction.cc
+++ b/Processes/Sibyll/Interaction.cc
@@ -35,7 +35,7 @@ using Track = Trajectory;
 namespace corsika::process::sibyll {
 
   Interaction::Interaction(environment::Environment const& env)
-      : fEnvironment(env) {}
+    : fEnvironment(env) {}
 
   Interaction::~Interaction() {
     cout << "Sibyll::Interaction n=" << fCount << " Nnuc=" << fNucCount << endl;
@@ -49,10 +49,46 @@ namespace corsika::process::sibyll {
     if (!fInitialized) {
       sibyll_ini_();
 
+      // any decays in sibyll? if yes need to define which particles
+      if( fInternalDecays){
+	// define which particles are passed to corsika, i.e. which particles make it into history
+	// even very shortlived particles like charm or pi0 are of interest here
+	const std::vector<particles::Code> HadronsWeWantTrackedByCorsika = {
+        particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::Pi0,
+        particles::Code::KMinus, particles::Code::KPlus,
+	particles::Code::K0Long,  particles::Code::K0Short,
+	particles::Code::SigmaPlus, particles::Code::SigmaMinus,
+	particles::Code::Lambda0,
+	particles::Code::Xi0, particles::Code::XiMinus,
+	particles::Code::OmegaMinus,
+	particles::Code::DPlus, particles::Code::DMinus, particles::Code::D0, particles::Code::D0Bar};
+	
+	Interaction::SetParticleListStable(HadronsWeWantTrackedByCorsika);
+      }
+      
       fInitialized = true;
     }
   }
 
+  void Interaction::SetParticleListStable(const std::vector<particles::Code> particleList) {
+    for (auto p : particleList)
+      Interaction::SetStable( p );
+  }
+
+  void Interaction::SetUnstable(const particles::Code pCode) {
+    cout << "Sibyll::Interaction: setting " << pCode << " unstable.." << endl;
+    int s_id = process::sibyll::ConvertToSibyllRaw(pCode);
+    s_csydec_.idb[s_id - 1] = abs(s_csydec_.idb[s_id - 1]);
+  }
+
+  void Interaction::SetStable(const particles::Code pCode) {
+    cout << "Sibyll::Interaction: setting " << pCode << " stable.." << endl;
+    int s_id = process::sibyll::ConvertToSibyllRaw(pCode);
+    s_csydec_.idb[s_id - 1] = (-1) * abs(s_csydec_.idb[s_id - 1]);
+  }
+
+  
+  
   tuple<units::si::CrossSectionType, units::si::CrossSectionType>
   Interaction::GetCrossSection(const particles::Code BeamId,
                                const particles::Code TargetId,
@@ -270,8 +306,8 @@ namespace corsika::process::sibyll {
             sigEla; // to avoid not used warning in array binding
       }
 
-      const auto targetCode = mediumComposition.SampleTarget(
-          cross_section_of_components, fRNG);
+      const auto targetCode =
+          mediumComposition.SampleTarget(cross_section_of_components, fRNG);
       cout << "Interaction: target selected: " << targetCode << endl;
       /*
         FOR NOW: allow nuclei with A<18 or protons only.
@@ -309,8 +345,7 @@ namespace corsika::process::sibyll {
         // running sibyll, filling stack
         sibyll_(kBeam, targetSibCode, sqs);
         // running decays
-        // setTrackedParticlesStable();
-        decsib_();
+        if (fInternalDecays) decsib_();
         // print final state
         int print_unit = 6;
         sib_list_(print_unit);
diff --git a/Processes/Sibyll/Interaction.h b/Processes/Sibyll/Interaction.h
index 509f50a61861d30e0797e86d5190abcf469cb5b8..0f2fc58efc2d27aa337501c9b7d5174c6ef3f460 100644
--- a/Processes/Sibyll/Interaction.h
+++ b/Processes/Sibyll/Interaction.h
@@ -36,6 +36,10 @@ namespace corsika::process::sibyll {
 
     void Init();
 
+    void SetParticleListStable(const std::vector<particles::Code>);
+    void SetUnstable(const corsika::particles::Code );
+    void SetStable(const corsika::particles::Code );
+
     bool WasInitialized() { return fInitialized; }
     bool IsValidCoMEnergy(corsika::units::si::HEPEnergyType ecm) {
       return (fMinEnergyCoM <= ecm) && (ecm <= fMaxEnergyCoM);
@@ -70,6 +74,8 @@ namespace corsika::process::sibyll {
     corsika::environment::Environment const& fEnvironment;
     corsika::random::RNG& fRNG =
         corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");
+
+    const bool fInternalDecays = true;
     const corsika::units::si::HEPEnergyType fMinEnergyCoM = 
       10. * 1e9 * corsika::units::si::electronvolt;
     const corsika::units::si::HEPEnergyType fMaxEnergyCoM = 
diff --git a/Processes/Sibyll/NuclearInteraction.cc b/Processes/Sibyll/NuclearInteraction.cc
index 0e65121d478942af665824f1348777a3764e21da..6c543c3581538976d0c6e84189dba272d773538b 100644
--- a/Processes/Sibyll/NuclearInteraction.cc
+++ b/Processes/Sibyll/NuclearInteraction.cc
@@ -448,7 +448,8 @@ namespace corsika::process::sibyll {
       [[maybe_unused]] auto sigElaCopy = sigEla; // ONLY TO AVOID COMPILER WARNINGS
     }
 
-    const auto targetCode = mediumComposition.SampleTarget(cross_section_of_components, fRNG);
+    const auto targetCode =
+        mediumComposition.SampleTarget(cross_section_of_components, fRNG);
     cout << "Interaction: target selected: " << targetCode << endl;
     /*
       FOR NOW: allow nuclei with A<18 or protons only.
diff --git a/Processes/Sibyll/testSibyll.cc b/Processes/Sibyll/testSibyll.cc
index 93a2840e04b220081b9f307854e56e02b65a6241..c9b2cddfa30489fe19904716ec7d990c6384c52c 100644
--- a/Processes/Sibyll/testSibyll.cc
+++ b/Processes/Sibyll/testSibyll.cc
@@ -170,7 +170,11 @@ TEST_CASE("SibyllInterface", "[processes]") {
                    corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
             particles::Code::Proton, E0, plab, pos, 0_ns});
 
-    Decay model;
+    const std::vector<particles::Code> particleList = {
+        particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
+        particles::Code::KMinus, particles::Code::K0Long,  particles::Code::K0Short};
+
+    Decay model(particleList);
 
     model.Init();
     /*[[maybe_unused]] const process::EProcessReturn ret =*/model.DoDecay(particle,