diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 2ef93a896e2a7343746c270b20ea6e32834087c8..fe2865b43182d3afac7a1513212bcd1cde988e13 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -364,13 +364,16 @@ build_test_example-clang-8:
     - set -o pipefail
     - make -j2 install
   rules:
+    - if: '$CI_MERGE_REQUEST_LABELS =~ /Ready for code review/' # run on merge requests, if label 'Ready for code review' is set
+    - if: $CI_COMMIT_BRANCH == $CI_DEFAULT_BRANCH
+      when: manual
+      allow_failure: true
     - if: $CI_MERGE_REQUEST_ID
       when: manual
+      allow_failure: true
     - if: $CI_COMMIT_TAG
       when: manual
-    - if: $CI_COMMIT_BRANCH == $CI_DEFAULT_BRANCH
-    - if: $CI_COMMIT_BRANCH
-      when: manual
+      allow_failure: true
   cache:
     paths:
       - ${CI_PROJECT_DIR}/build/
@@ -549,6 +552,7 @@ sanity:
     - corsika
   before_script:
     - apt-get update && apt-get install -y python3-pip
+    - pip3 install --upgrade cython
   script:
     - cd ${CI_PROJECT_DIR}/python  # change into the Python directory
     - pip3 install --user -e '.[test]'  # install the package + test deps
diff --git a/corsika/detail/framework/core/Cascade.inl b/corsika/detail/framework/core/Cascade.inl
index 654dfc6d1aad0b603baa76e701f5357fa85ffc13..eb09a2c763ccf49015ce2d5e22539a1342ee5b8c 100644
--- a/corsika/detail/framework/core/Cascade.inl
+++ b/corsika/detail/framework/core/Cascade.inl
@@ -31,9 +31,6 @@ namespace corsika {
   inline void Cascade<TTracking, TProcessList, TOutput, TStack, TStackView>::run() {
     setNodes(); // put each particle on stack in correct environment volume
 
-    // start this event (i.e. this shower)
-    output_.startOfShower();
-
     while (!stack_.isEmpty()) {
       while (!stack_.isEmpty()) {
         CORSIKA_LOG_TRACE("Stack: {}", stack_.asString());
@@ -54,9 +51,6 @@ namespace corsika {
       // thus, the double loop
       // doCascadeEquations();
     }
-
-    // end this event (i.e. this shower)
-    output_.endOfShower();
   }
 
   template <typename TTracking, typename TProcessList, typename TOutput, typename TStack,
diff --git a/corsika/detail/framework/process/BoundaryCrossingProcess.hpp b/corsika/detail/framework/process/BoundaryCrossingProcess.hpp
index 1efd36454579d35b31d9f97421c42ca42a31674c..bf0e64dfa3643d7f7cce6cca61884ff1ab2c5c86 100644
--- a/corsika/detail/framework/process/BoundaryCrossingProcess.hpp
+++ b/corsika/detail/framework/process/BoundaryCrossingProcess.hpp
@@ -9,6 +9,7 @@
 #pragma once
 
 #include <corsika/framework/process/ProcessTraits.hpp>
+#include <corsika/framework/utility/HasMethodSignature.hpp>
 
 namespace corsika {
 
diff --git a/corsika/detail/framework/process/ContinuousProcess.hpp b/corsika/detail/framework/process/ContinuousProcess.hpp
index 451f45a68be3c592e3ac4d68e58f78cab7f682d3..68e5185bbdd55115d6e5174ee208ef4b6e3c9c70 100644
--- a/corsika/detail/framework/process/ContinuousProcess.hpp
+++ b/corsika/detail/framework/process/ContinuousProcess.hpp
@@ -9,6 +9,7 @@
 #pragma once
 
 #include <corsika/framework/process/ProcessTraits.hpp>
+#include <corsika/framework/utility/HasMethodSignature.hpp>
 
 namespace corsika {
 
diff --git a/corsika/detail/framework/process/DecayProcess.hpp b/corsika/detail/framework/process/DecayProcess.hpp
index b5664418a4d6db1692be00294c82d7e224782dbc..048466e1513e140e67de2ba685928ba70a537635 100644
--- a/corsika/detail/framework/process/DecayProcess.hpp
+++ b/corsika/detail/framework/process/DecayProcess.hpp
@@ -9,6 +9,7 @@
 #pragma once
 
 #include <corsika/framework/process/ProcessTraits.hpp>
+#include <corsika/framework/utility/HasMethodSignature.hpp>
 
 namespace corsika {
 
diff --git a/corsika/detail/framework/process/InteractionProcess.hpp b/corsika/detail/framework/process/InteractionProcess.hpp
index 8fb157a688c8c4828ea5c0b3eff410af9ad9aa76..2446a385d7beec2424da598585dd89bd70466a5c 100644
--- a/corsika/detail/framework/process/InteractionProcess.hpp
+++ b/corsika/detail/framework/process/InteractionProcess.hpp
@@ -9,6 +9,7 @@
 #pragma once
 
 #include <corsika/framework/process/ProcessTraits.hpp>
+#include <corsika/framework/utility/HasMethodSignature.hpp>
 
 namespace corsika {
 
@@ -45,7 +46,7 @@ namespace corsika {
     //! @}
   };
 
-  //! @file BoundaryCrossingProcess.hpp
+  //! @file InteractionProcess.hpp
   //! value traits type
   template <class TProcess, typename TReturn, typename... TArgs>
   bool constexpr has_method_doInteract_v =
@@ -85,7 +86,7 @@ namespace corsika {
     //! @}
   };
 
-  //! @file BoundaryCrossingProcess.hpp
+  //! @file InteractionProcess.hpp
   //! value traits type
 
   template <class TProcess, typename TReturn, typename... TArgs>
diff --git a/corsika/detail/framework/process/SecondariesProcess.hpp b/corsika/detail/framework/process/SecondariesProcess.hpp
index 0310bf3dfa2b73826c10c2ad25bf4eb2d4ab5a21..ee76cce853d29f7007d9a82ed88b12e0826a31f6 100644
--- a/corsika/detail/framework/process/SecondariesProcess.hpp
+++ b/corsika/detail/framework/process/SecondariesProcess.hpp
@@ -9,6 +9,7 @@
 #pragma once
 
 #include <corsika/framework/process/ProcessTraits.hpp>
+#include <corsika/framework/utility/HasMethodSignature.hpp>
 
 namespace corsika {
 
diff --git a/corsika/detail/framework/process/StackProcess.hpp b/corsika/detail/framework/process/StackProcess.hpp
index 68f350adf1bd1db6cc2be74b329e6a4824b2df07..f06126ed410fef6412b529f0fcc516490212bc18 100644
--- a/corsika/detail/framework/process/StackProcess.hpp
+++ b/corsika/detail/framework/process/StackProcess.hpp
@@ -9,6 +9,7 @@
 #pragma once
 
 #include <corsika/framework/process/ProcessTraits.hpp>
+#include <corsika/framework/utility/HasMethodSignature.hpp>
 
 namespace corsika {
 
diff --git a/corsika/detail/modules/proposal/ContinuousProcess.inl b/corsika/detail/modules/proposal/ContinuousProcess.inl
index 91308390ba4f9b5b8e87bb0d94e79646e9a6f101..b8403dc145a85f212ceb95699443e56916c7336a 100644
--- a/corsika/detail/modules/proposal/ContinuousProcess.inl
+++ b/corsika/detail/modules/proposal/ContinuousProcess.inl
@@ -1,3 +1,4 @@
+
 /*
  * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
  *
diff --git a/corsika/detail/setup/SetupStack.hpp b/corsika/detail/setup/SetupStack.hpp
index e58ff5d9c7a7e88450147d8761e22585c1de8512..03016a47ebabcab2d2fd7b087e4f1a8f401b19a0 100644
--- a/corsika/detail/setup/SetupStack.hpp
+++ b/corsika/detail/setup/SetupStack.hpp
@@ -36,9 +36,9 @@ namespace corsika {
         CombinedParticleInterface<nuclear_stack::ParticleDataStack::pi_type,
                                   SetupGeometryDataInterface, TStackIter>;
 
-    using StackWithGeometry = CombinedStack<
-        typename nuclear_stack::ParticleDataStack::stack_implementation_type,
-        node::GeometryData<setup::Environment>, StackWithGeometryInterface>;
+    using StackWithGeometry =
+        CombinedStack<typename nuclear_stack::ParticleDataStack::stack_data_type,
+                      node::GeometryData<setup::Environment>, StackWithGeometryInterface>;
 
     // ------------------------------------------
     // Add [optional] history data to stack, too:
@@ -50,7 +50,7 @@ namespace corsika {
                                   history::HistoryEventDataInterface, TStackIter>;
 
     using StackWithHistory =
-        CombinedStack<typename StackWithGeometry::stack_implementation_type,
+        CombinedStack<typename StackWithGeometry::stack_data_type,
                       history::HistoryEventData, StackWithHistoryInterface>;
 
   } // namespace setup::detail
diff --git a/corsika/framework/process/ProcessTraits.hpp b/corsika/framework/process/ProcessTraits.hpp
index e675fb1717e0c848b91b78ea6b3e0730d2063cf4..cccd80abaf1a11436befc913eaaf785112c5b123 100644
--- a/corsika/framework/process/ProcessTraits.hpp
+++ b/corsika/framework/process/ProcessTraits.hpp
@@ -96,33 +96,4 @@ namespace corsika {
     static unsigned int constexpr count = N;
   };
 
-  namespace detail {
-
-    /**
-       Helper traits class (partial) for static compile time checking.
-
-       Note, this is a poor replacement for C++20 concepts... they are
-       eagerly awaited!
-
-       It defines the default body of a generic test function returning
-       std::false_type.
-
-       In addition it defines the pattern for class-method matching with a
-       return type TReturn and function arguments TArgs... . Right now
-       both method signatures, "const" and "not const", are matched.
-     */
-    template <typename TReturn, typename... TArgs>
-    struct has_method_signature {
-
-      // the non-const version
-      template <class T>
-      static std::true_type testSignature(TReturn (T::*)(TArgs...));
-
-      // the const version
-      template <class T>
-      static std::true_type testSignature(TReturn (T::*)(TArgs...) const);
-    };
-
-  } // namespace detail
-
 } // namespace corsika
diff --git a/corsika/framework/stack/ParticleBase.hpp b/corsika/framework/stack/ParticleBase.hpp
index 2c62285d5fbcd26d4b06841809cd2a0f8bf8a602..3ddb7fba4a5be42499899faed4bb87b657ed0bdf 100644
--- a/corsika/framework/stack/ParticleBase.hpp
+++ b/corsika/framework/stack/ParticleBase.hpp
@@ -8,6 +8,7 @@
 
 #pragma once
 
+#include <type_traits>
 #include <cstdlib> // for size_t
 
 namespace corsika {
@@ -17,8 +18,8 @@ namespace corsika {
    particle stack. Every stack must implement this readout via the
    ParticleBase class.
 
-   The StackIterator template argument is derived from StackIteratorInterface, which is of
-   type <code> template <typename StackData, template <typename> typename
+   The TStackIterator template argument is derived from StackIteratorInterface, which is
+   of type <code> template <typename StackData, template <typename> typename
    ParticleInterface> class StackIteratorInterface : public
    ParticleInterface<StackIteratorInterface<StackData, ParticleInterface>>
    </code>
@@ -38,11 +39,11 @@ namespace corsika {
 
   */
 
-  template <typename StackIterator>
+  template <typename TStackIterator>
   class ParticleBase {
 
   public:
-    typedef StackIterator stack_iterator_type;
+    typedef TStackIterator stack_iterator_type;
 
     ParticleBase() = default;
 
@@ -75,14 +76,13 @@ namespace corsika {
      */
     template <typename... TArgs>
     stack_iterator_type addSecondary(const TArgs... args) {
-
       return this->getStack().addSecondary(this->getIterator(), args...);
     }
 
     // protected: // todo should [MAY]be proteced, but don't now how to 'friend Stack'
     // Function to provide CRTP access to inheriting class (type)
     /**
-     * return the corresponding StackIterator for this particle
+     * return the corresponding TStackIterator for this particle
      */
     stack_iterator_type& getIterator() {
       return static_cast<stack_iterator_type&>(*this);
diff --git a/corsika/framework/stack/SecondaryView.hpp b/corsika/framework/stack/SecondaryView.hpp
index 9e5f60324c17f7a4feb144a5a8b638b74fd88a8c..6ba881b9262be89f403f58aed3338dd15fd9f218 100644
--- a/corsika/framework/stack/SecondaryView.hpp
+++ b/corsika/framework/stack/SecondaryView.hpp
@@ -430,8 +430,8 @@ namespace corsika {
             class MSecondaryProducer = corsika::DefaultSecondaryProducer,
             template <typename> typename pi_type_ = TStack::template pi_type>
   struct MakeView {
-    using type = corsika::SecondaryView<typename TStack::stack_implementation_type,
-                                        pi_type_, MSecondaryProducer>;
+    using type = corsika::SecondaryView<typename TStack::stack_data_type, pi_type_,
+                                        MSecondaryProducer>;
   };
 #endif
 
diff --git a/corsika/framework/stack/Stack.hpp b/corsika/framework/stack/Stack.hpp
index caaf9d065b709fc0161c87f41bd2ac3661f8b2ac..81da792924c4c9d4e3197a44562e93ef90db5ed4 100644
--- a/corsika/framework/stack/Stack.hpp
+++ b/corsika/framework/stack/Stack.hpp
@@ -72,17 +72,17 @@ namespace corsika {
      loops, ranges, etc.
    */
 
-  template <typename StackData, template <typename> typename MParticleInterface>
+  template <typename TStackData, template <typename> typename MParticleInterface>
   class Stack {
 
-    typedef typename std::remove_reference<StackData>::type value_type;
+    typedef typename std::remove_reference<TStackData>::type value_type;
 
   public:
-    typedef StackData stack_implementation_type; ///< this is the type of the
-                                                 ///< user-provided data structure
+    typedef TStackData stack_data_type; ///< this is the type of the
+                                        ///< user-provided data structure
 
-    template <typename TSI>
-    using pi_type = MParticleInterface<TSI>;
+    template <typename TStackIterator>
+    using pi_type = MParticleInterface<TStackIterator>; // @todo pi_type -> pi_template
 
     /**
      * Via the StackIteratorInterface and ConstStackIteratorInterface
@@ -122,26 +122,26 @@ namespace corsika {
         delete; ///< since Stack can be very big, we don't want to copy it
 
     /**
-     * if StackData is a reference member we *HAVE* to initialize
+     * if TStackData is a reference member we *HAVE* to initialize
      * it in the constructor, this is typically needed for SecondaryView
      */
-    template <typename UType = StackData,
+    template <typename UType = TStackData,
               typename = typename std::enable_if<std::is_reference<UType>::value>::type>
-    Stack(StackData vD)
+    Stack(TStackData vD)
         : nDeleted_(0)
         , data_(vD)
         , deleted_(std::vector<bool>(data_.getSize(), false)) {}
 
     /**
      * This constructor takes any argument and passes it on to the
-     * StackData user class. If the user did not provide a suited
+     * TStackData user class. If the user did not provide a suited
      * constructor this will fail with an error message.
      *
      * Furthermore, this is disabled with enable_if for SecondaryView
      * stacks, where the inner data container is always a reference
      * and cannot be initialized here.
      */
-    template <typename... TArgs, typename UType = StackData,
+    template <typename... TArgs, typename UType = TStackData,
               typename = typename std::enable_if<std::is_reference<UType>::value>::type>
     Stack(TArgs... args)
         : nDeleted_(0)
@@ -149,7 +149,7 @@ namespace corsika {
         , deleted_(std::vector<bool>(data_.getSize(), false)) {}
 
     /**
-     * @name Most generic proxy methods for StackData data_
+     * @name Most generic proxy methods for TStackData data_
      * @{
      */
     unsigned int getCapacity() const { return data_.getCapacity(); }
@@ -197,7 +197,10 @@ namespace corsika {
     /**
      * increase stack size, create new particle at end of stack
      */
-    template <typename... TArgs>
+    template <typename... TArgs> //,
+    //        typename = std::enable_if_t<std::is_same_v<
+    //      void, std::invoke_result_t<decltype((**stack_iterator_type).setParticleData),
+    //                                 TArgs...>>>>
     stack_iterator_type addParticle(const TArgs... v);
 
     void swap(stack_iterator_type a, stack_iterator_type b);
@@ -277,12 +280,12 @@ namespace corsika {
     /**
      * Function to perform eventual transformation from
      * StackIterator::getIndex() to index in data stored in
-     * StackData data_. By default (and in almost all cases) this
+     * TStackData data_. By default (and in almost all cases) this
      * should just be identiy. See class SecondaryView for an alternative implementation.
      */
     unsigned int getIndexFromIterator(const unsigned int vI) const;
     /**
-     * @name Return reference to StackData object data_ for data access
+     * @name Return reference to TStackData object data_ for data access
      * @{
      */
 
@@ -306,7 +309,7 @@ namespace corsika {
     unsigned int nDeleted_ = 0;
 
   private:
-    StackData data_; ///< this in general holds all the data and can be quite big
+    TStackData data_; ///< this in general holds all the data and can be quite big
     std::vector<bool> deleted_; ///< bit field to flag deleted entries
   };
 
diff --git a/corsika/framework/stack/StackIteratorInterface.hpp b/corsika/framework/stack/StackIteratorInterface.hpp
index 02ef7e5423c99e64248d4a5135f2d8262c221eb8..d6a267014def6b44d46011cf5c46f0791df5e76e 100644
--- a/corsika/framework/stack/StackIteratorInterface.hpp
+++ b/corsika/framework/stack/StackIteratorInterface.hpp
@@ -10,6 +10,8 @@
 
 #include <corsika/framework/stack/ParticleBase.hpp>
 
+#include <boost/type_index.hpp>
+
 namespace corsika::history {
   template <typename T, template <typename> typename TParticleInterface>
   class HistorySecondaryProducer; // forward decl.
@@ -77,6 +79,9 @@ namespace corsika {
         corsika::StackIteratorInterface<TStackData, TParticleInterface, TStackType>>
         particle_interface_type;
 
+    typedef TStackType stack_type;
+    typedef TStackData stack_data_type;
+
     // it is not allowed to create a "dangling" stack iterator
     StackIteratorInterface() = delete; //! \todo check rule of five
 
@@ -100,7 +105,7 @@ namespace corsika {
           @param data reference to the stack [rw]
           @param index index on stack
        */
-    StackIteratorInterface(TStackType& data, unsigned int const index)
+    StackIteratorInterface(stack_type& data, unsigned int const index)
         : index_(index)
         , data_(&data) {}
 
@@ -112,10 +117,11 @@ namespace corsika {
        particle_interface_type::setParticleData(...) function
      */
     template <typename... TArgs>
-    StackIteratorInterface(TStackType& data, unsigned int const index,
+    StackIteratorInterface(stack_type& data, unsigned int const index,
                            const TArgs... args)
         : index_(index)
         , data_(&data) {
+
       (**this).setParticleData(args...);
     }
 
@@ -129,11 +135,12 @@ namespace corsika {
        consistent with the definition of the user-provided
        particle_interface_type::setParticleData(...) function
     */
-    template <typename... Args>
-    StackIteratorInterface(TStackType& data, unsigned int const index,
-                           StackIteratorInterface& parent, const Args... args)
+    template <typename... TArgs>
+    StackIteratorInterface(stack_type& data, unsigned int const index,
+                           StackIteratorInterface& parent, const TArgs... args)
         : index_(index)
         , data_(&data) {
+
       (**this).setParticleData(*parent, args...);
     }
 
@@ -168,10 +175,10 @@ namespace corsika {
       return index_ != rhs.index_;
     }
     bool operator==(
-        const ConstStackIteratorInterface<TStackData, TParticleInterface, TStackType>&
+        const ConstStackIteratorInterface<TStackData, TParticleInterface, stack_type>&
             rhs) const; // implemented below
     bool operator!=(
-        const ConstStackIteratorInterface<TStackData, TParticleInterface, TStackType>&
+        const ConstStackIteratorInterface<TStackData, TParticleInterface, stack_type>&
             rhs) const; // implemented below
 
     /**
@@ -199,9 +206,9 @@ namespace corsika {
     /// Get current particle index
     unsigned int getIndex() const { return index_; }
     /// Get current particle Stack object
-    TStackType& getStack() { return *data_; }
+    stack_type& getStack() { return *data_; }
     /// Get current particle const Stack object
-    TStackType const& getStack() const { return *data_; }
+    stack_type const& getStack() const { return *data_; }
     /// Get current user particle TStackData object
     TStackData& getStackData() { return data_->getStackData(); }
     /// Get current const user particle TStackData object
@@ -227,13 +234,13 @@ namespace corsika {
     template <typename T, template <typename> typename TParticleInterface_>
     friend class corsika::history::HistorySecondaryProducer;
 
-    friend class ConstStackIteratorInterface<TStackData, TParticleInterface, TStackType>;
+    friend class ConstStackIteratorInterface<TStackData, TParticleInterface, stack_type>;
 
   protected:
     unsigned int index_ = 0;
 
   private:
-    TStackType* data_ = 0; // info: Particles and StackIterators become invalid when
+    stack_type* data_ = 0; // info: Particles and StackIterators become invalid when
                            // parent Stack is copied or deleted!
 
   }; // end class StackIterator
@@ -263,6 +270,9 @@ namespace corsika {
         ConstStackIteratorInterface<TStackData, TParticleInterface, TStackType>>
         particle_interface_type;
 
+    typedef TStackType stack_type;
+    typedef TStackData stack_data_type;
+
     // we don't want to allow dangling iterators to exist
     ConstStackIteratorInterface() = delete; //! \todo check rule of five
 
@@ -271,7 +281,7 @@ namespace corsika {
         : index_(std::move(rhs.index_))
         , data_(std::move(rhs.data_)) {}
 
-    ConstStackIteratorInterface(TStackType const& data, unsigned int const index)
+    ConstStackIteratorInterface(stack_type const& data, unsigned int const index)
         : index_(index)
         , data_(&data) {}
 
@@ -304,12 +314,12 @@ namespace corsika {
     bool operator!=(ConstStackIteratorInterface const& rhs) const {
       return index_ != rhs.index_;
     }
-    bool operator==(StackIteratorInterface<TStackData, TParticleInterface,
-                                           TStackType> const& rhs) const {
+    bool operator==(StackIteratorInterface<stack_data_type, TParticleInterface,
+                                           stack_type> const& rhs) const {
       return index_ == rhs.index_;
     }
-    bool operator!=(StackIteratorInterface<TStackData, TParticleInterface,
-                                           TStackType> const& rhs) const {
+    bool operator!=(StackIteratorInterface<stack_data_type, TParticleInterface,
+                                           stack_type> const& rhs) const {
       return index_ != rhs.index_;
     }
 
@@ -324,8 +334,8 @@ namespace corsika {
         @{
      */
     unsigned int getIndex() const { return index_; }
-    TStackType const& getStack() const { return *data_; }
-    TStackData const& getStackData() const { return data_->getStackData(); }
+    stack_type const& getStack() const { return *data_; }
+    stack_data_type const& getStackData() const { return data_->getStackData(); }
     /// Get data index as mapped in Stack class
     unsigned int getIndexFromIterator() const {
       return data_->getIndexFromIterator(index_);
@@ -333,18 +343,18 @@ namespace corsika {
     ///@}
 
     // friends are needed for access to protected methods
-    friend class Stack<TStackData,
+    friend class Stack<stack_data_type,
                        TParticleInterface>; // for access to GetIndex for Stack
-    friend class Stack<TStackData&, TParticleInterface>; // for access to GetIndex
+    friend class Stack<stack_data_type&, TParticleInterface>; // for access to GetIndex
 
     friend class ParticleBase<ConstStackIteratorInterface>; // for access to GetStackData
 
-    template <typename T1,                     // best fix to: TStackData,
+    template <typename T1,                     // best fix to: stack_data_type,
               template <typename> typename M1, // best fix to: TParticleInterface,
               template <class T2, template <class> class T3> class MSecondaryProducer>
     friend class SecondaryView; // access for SecondaryView
 
-    friend class StackIteratorInterface<TStackData, TParticleInterface, TStackType>;
+    friend class StackIteratorInterface<stack_data_type, TParticleInterface, stack_type>;
 
     template <typename T, template <typename> typename TParticleInterface_>
     friend class corsika::history::HistorySecondaryProducer;
@@ -353,7 +363,7 @@ namespace corsika {
     unsigned int index_ = 0;
 
   private:
-    TStackType const* data_ = 0; // info: Particles and StackIterators become invalid when
+    stack_type const* data_ = 0; // info: Particles and StackIterators become invalid when
                                  // parent Stack is copied or deleted!
 
   }; // end class ConstStackIterator
diff --git a/corsika/framework/utility/HasMethodSignature.hpp b/corsika/framework/utility/HasMethodSignature.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..a3a06567c8078594e0f308fbbb58ed995e3bf33a
--- /dev/null
+++ b/corsika/framework/utility/HasMethodSignature.hpp
@@ -0,0 +1,47 @@
+#pragma once
+
+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * 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.
+ */
+
+/**
+ * @file HasMethodSignature.hpp
+ */
+
+#include <type_traits>
+
+namespace corsika {
+
+  namespace detail {
+
+    /**
+       Helper traits class (partial) for static compile time checking.
+
+       Note, this is a poor replacement for C++20 concepts... they are
+       eagerly awaited!
+
+       It defines the default body of a generic test function returning
+       std::false_type.
+
+       In addition it defines the pattern for class-method matching with a
+       return type TReturn and function arguments TArgs... . Right now
+       both method signatures, "const" and "not const", are matched.
+    */
+    template <typename TReturn, typename... TArgs>
+    struct has_method_signature {
+
+      // the non-const version
+      template <class T>
+      static std::true_type testSignature(TReturn (T::*)(TArgs...));
+
+      // the const version
+      template <class T>
+      static std::true_type testSignature(TReturn (T::*)(TArgs...) const);
+    };
+
+  } // namespace detail
+} // namespace corsika
diff --git a/corsika/modules/qgsjetII/QGSJetIIFragmentsStack.hpp b/corsika/modules/qgsjetII/QGSJetIIFragmentsStack.hpp
index 5539ec807a99804c5c66fe343c16bc7948a0c3a4..8d1279382c244b7c0265a31dbf5ebb7413fe837e 100644
--- a/corsika/modules/qgsjetII/QGSJetIIFragmentsStack.hpp
+++ b/corsika/modules/qgsjetII/QGSJetIIFragmentsStack.hpp
@@ -50,16 +50,16 @@ namespace corsika::qgsjetII {
     }
   };
 
-  template <typename StackIteratorInterface>
-  class FragmentsInterface : public corsika::ParticleBase<StackIteratorInterface> {
+  template <typename TStackIterator>
+  class FragmentsInterface : public corsika::ParticleBase<TStackIterator> {
 
-    using corsika::ParticleBase<StackIteratorInterface>::getStackData;
-    using corsika::ParticleBase<StackIteratorInterface>::getIndex;
+    using corsika::ParticleBase<TStackIterator>::getStackData;
+    using corsika::ParticleBase<TStackIterator>::getIndex;
 
   public:
     void setParticleData(const int vSize) { setFragmentSize(vSize); }
 
-    void setParticleData(FragmentsInterface<StackIteratorInterface>& /*parent*/,
+    void setParticleData(FragmentsInterface<TStackIterator>& /*parent*/,
                          const int vSize) {
       setFragmentSize(vSize);
     }
diff --git a/corsika/modules/qgsjetII/QGSJetIIStack.hpp b/corsika/modules/qgsjetII/QGSJetIIStack.hpp
index bde5bdeacdc219f8ff2080f4ab23c2a59ebe2eb4..fa78938422c504be7a2a3695d434905ee074ee24 100644
--- a/corsika/modules/qgsjetII/QGSJetIIStack.hpp
+++ b/corsika/modules/qgsjetII/QGSJetIIStack.hpp
@@ -44,17 +44,17 @@ namespace corsika::qgsjetII {
     void decrementSize();
   };
 
-  template <typename StackIteratorInterface>
-  class ParticleInterface : public corsika::ParticleBase<StackIteratorInterface> {
+  template <typename TStackIterator>
+  class ParticleInterface : public corsika::ParticleBase<TStackIterator> {
 
-    using corsika::ParticleBase<StackIteratorInterface>::getStackData;
-    using corsika::ParticleBase<StackIteratorInterface>::getIndex;
+    using corsika::ParticleBase<TStackIterator>::getStackData;
+    using corsika::ParticleBase<TStackIterator>::getIndex;
 
   public:
     void setParticleData(const int vID, const HEPEnergyType vE, const MomentumVector& vP,
                          const HEPMassType);
-    void setParticleData(ParticleInterface<StackIteratorInterface>& /*parent*/,
-                         const int vID, const HEPEnergyType vE, const MomentumVector& vP,
+    void setParticleData(ParticleInterface<TStackIterator>& /*parent*/, const int vID,
+                         const HEPEnergyType vE, const MomentumVector& vP,
                          const HEPMassType);
 
     void setEnergy(const HEPEnergyType v);
diff --git a/corsika/modules/sibyll/SibStack.hpp b/corsika/modules/sibyll/SibStack.hpp
index 68e2d63d54c2314d03888276c123cc01129c15aa..636180d206f3486b3276ae25c89f47ba44800ac6 100644
--- a/corsika/modules/sibyll/SibStack.hpp
+++ b/corsika/modules/sibyll/SibStack.hpp
@@ -68,15 +68,14 @@ namespace corsika::sibyll {
     }
   };
 
-  template <typename StackIteratorInterface>
-  class ParticleInterface : public corsika::ParticleBase<StackIteratorInterface> {
+  template <typename TStackIterator>
+  class ParticleInterface : public corsika::ParticleBase<TStackIterator> {
 
-    using corsika::ParticleBase<StackIteratorInterface>::getStackData;
-    using corsika::ParticleBase<StackIteratorInterface>::getIndex;
+    using corsika::ParticleBase<TStackIterator>::getStackData;
+    using corsika::ParticleBase<TStackIterator>::getIndex;
 
   public:
-    void setParticleData(const int vID, // corsika::sibyll::SibyllCode vID,
-                         const HEPEnergyType vE, const MomentumVector& vP,
+    void setParticleData(const int vID, const HEPEnergyType vE, const MomentumVector& vP,
                          const HEPMassType vM) {
       setPID(vID);
       setEnergy(vE);
@@ -84,8 +83,7 @@ namespace corsika::sibyll {
       setMass(vM);
     }
 
-    void setParticleData(ParticleInterface<StackIteratorInterface>& /*parent*/,
-                         const int vID, //  corsika::sibyll::SibyllCode vID,
+    void setParticleData(ParticleInterface<TStackIterator>& /*parent*/, const int vID,
                          const HEPEnergyType vE, const MomentumVector& vP,
                          const HEPMassType vM) {
       setPID(vID);
diff --git a/corsika/setup/SetupStack.hpp b/corsika/setup/SetupStack.hpp
index 847f83834e57d1afa757151df9e92cce7ee8dca4..a051e4db177cd36e087f42612652ffe0737c585d 100644
--- a/corsika/setup/SetupStack.hpp
+++ b/corsika/setup/SetupStack.hpp
@@ -57,7 +57,7 @@ namespace corsika::setup {
 #ifdef WITH_HISTORY
 
 #if defined(__clang__)
-  using StackView = SecondaryView<typename Stack::stack_implementation_type,
+  using StackView = SecondaryView<typename Stack::stack_data_type,
                                   // CHECK with CLANG: setup::Stack::MPIType>;
                                   detail::StackWithHistoryInterface, StackViewProducer>;
 #elif defined(__GNUC__) || defined(__GNUG__)
@@ -67,7 +67,7 @@ namespace corsika::setup {
 #else // WITH_HISTORY
 
 #if defined(__clang__)
-  using StackView = SecondaryView<typename setup::Stack::stack_implementation_type,
+  using StackView = SecondaryView<typename setup::Stack::stack_data_type,
                                   // CHECK with CLANG:
                                   // setup::Stack::MPIType>;
                                   setup::detail::StackWithGeometryInterface>;
diff --git a/corsika/stack/NuclearStackExtension.hpp b/corsika/stack/NuclearStackExtension.hpp
index 6c928656642959d93bf07a469a98e04637884bc8..03870ccf5de4b0225be86add758eaec8ef0c8a15 100644
--- a/corsika/stack/NuclearStackExtension.hpp
+++ b/corsika/stack/NuclearStackExtension.hpp
@@ -291,8 +291,7 @@ namespace corsika::nuclear_stack {
 
   template <typename TInnerStack, template <typename> typename PI_>
   using NuclearStackExtension =
-      Stack<NuclearStackExtensionImpl<typename TInnerStack::stack_implementation_type>,
-            PI_>;
+      Stack<NuclearStackExtensionImpl<typename TInnerStack::stack_data_type>, PI_>;
 
   //
   template <typename TStackIter>
diff --git a/corsika/stack/VectorStack.hpp b/corsika/stack/VectorStack.hpp
index 223d651f90a6a3d37c5d45e362fa8702e21e9c54..154a9850bcb11e8157c8c2e9bc29c862f907184c 100644
--- a/corsika/stack/VectorStack.hpp
+++ b/corsika/stack/VectorStack.hpp
@@ -26,11 +26,11 @@ namespace corsika {
    * Example of a particle object on the stack.
    */
 
-  template <typename StackIteratorInterface>
-  class ParticleInterface : public ParticleBase<StackIteratorInterface> {
+  template <typename TStackIterator>
+  class ParticleInterface : public ParticleBase<TStackIterator> {
 
   private:
-    typedef ParticleBase<StackIteratorInterface> super_type;
+    typedef ParticleBase<TStackIterator> super_type;
 
   public:
     typedef std::tuple<Code, HEPEnergyType, DirectionVector, Point, TimeType>
@@ -59,7 +59,7 @@ namespace corsika {
      *  MomentumVector is only used to determine the DirectionVector, the normalization
      * is lost.
      */
-    void setParticleData(ParticleInterface<StackIteratorInterface> const& p,
+    void setParticleData(ParticleInterface<TStackIterator> const& p,
                          particle_data_type const& v);
 
     /**
@@ -77,7 +77,7 @@ namespace corsika {
      * @param v tuple containing: PID, kinetic Energy, Direction Vector, Position, Time
      *
      */
-    void setParticleData(ParticleInterface<StackIteratorInterface> const& p,
+    void setParticleData(ParticleInterface<TStackIterator> const& p,
                          particle_data_momentum_type const& v);
 
     ///! Set particle corsika::Code
diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index 3240a632cb15c390454e99e0515c52f7945336e7..2d20738b31793611db793fa90bf1cc1bca1b1dd5 100644
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -35,7 +35,7 @@ CORSIKA_REGISTER_EXAMPLE (cascade_proton_example)
 
 add_executable (vertical_EAS vertical_EAS.cpp)
 target_link_libraries (vertical_EAS CORSIKA8::CORSIKA8)
-CORSIKA_REGISTER_EXAMPLE (vertical_EAS RUN_OPTIONS 4 2 10000. 1)
+CORSIKA_REGISTER_EXAMPLE (vertical_EAS RUN_OPTIONS 4 2 10000.)
 
 add_executable (stopping_power stopping_power.cpp)
 target_link_libraries (stopping_power CORSIKA8::CORSIKA8)
@@ -56,3 +56,7 @@ CORSIKA_REGISTER_EXAMPLE (em_shower RUN_OPTIONS "100.")
 add_executable (hybrid_MC hybrid_MC.cpp)
 target_link_libraries (hybrid_MC  CORSIKA8::CORSIKA8)
 CORSIKA_REGISTER_EXAMPLE (hybrid_MC RUN_OPTIONS 4 2 10000.)
+
+add_executable (corsika corsika.cpp)
+target_link_libraries (corsika CORSIKA8::CORSIKA8)
+CORSIKA_REGISTER_EXAMPLE (corsika RUN_OPTIONS 4 2 10000. 1)
diff --git a/examples/boundary_example.cpp b/examples/boundary_example.cpp
index 3ce8df868981a77af9d4df509f414264a2c3138d..e5f707bdef84d99c6a553ad32b58ff0e949f8b86 100644
--- a/examples/boundary_example.cpp
+++ b/examples/boundary_example.cpp
@@ -26,7 +26,6 @@
 #include <corsika/media/UniformMagneticField.hpp>
 #include <corsika/media/MediumPropertyModel.hpp>
 
-#include <corsika/modules/Sibyll.hpp>
 #include <corsika/modules/TrackWriter.hpp>
 #include <corsika/modules/ParticleCut.hpp>
 
@@ -38,9 +37,6 @@
   executable. If you include the header below multiple times and
   link this togehter, it will fail.
  */
-#include <corsika/modules/sibyll/Random.hpp>
-#include <corsika/modules/urqmd/Random.hpp>
-
 #include <iostream>
 #include <limits>
 #include <typeinfo>
@@ -83,7 +79,6 @@ private:
 int main() {
 
   logging::set_level(logging::level::info);
-  corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
 
   CORSIKA_LOG_INFO("boundary_example");
 
@@ -120,10 +115,6 @@ int main() {
   // setup processes, decays and interactions
   setup::Tracking tracking;
 
-  RNGManager::getInstance().registerRandomStream("sibyll");
-  corsika::sibyll::Interaction sibyll;
-  corsika::sibyll::Decay decay;
-
   ParticleCut cut(50_GeV, true, true);
 
   TrackWriter trackWriter;
@@ -132,7 +123,7 @@ int main() {
   MyBoundaryCrossingProcess<true> boundaryCrossing("crossings.dat");
 
   // assemble all processes into an ordered process list
-  auto sequence = make_sequence(sibyll, decay, cut, boundaryCrossing, trackWriter);
+  auto sequence = make_sequence(cut, boundaryCrossing, trackWriter);
 
   // setup particle stack, and add primary particles
   setup::Stack stack;
@@ -173,7 +164,9 @@ int main() {
   // define air shower object, run simulation
   Cascade EAS(env, tracking, sequence, output, stack);
 
+  output.startOfShower();
   EAS.run();
+  output.endOfShower();
 
   CORSIKA_LOG_INFO("Result: E0={}GeV", E0 / 1_GeV);
   cut.showResults();
diff --git a/examples/cascade_example.cpp b/examples/cascade_example.cpp
index 2852428d2c074b0afaa22fcdc58ec5eeeaa0825c..623f443120ec69c52e5bacbecf6dfc02039e87e9 100644
--- a/examples/cascade_example.cpp
+++ b/examples/cascade_example.cpp
@@ -153,7 +153,9 @@ int main() {
   // define air shower object, run simulation
   Cascade EAS(env, tracking, sequence, output, stack);
 
+  output.startOfShower();
   EAS.run();
+  output.endOfShower();
 
   eLoss.printProfile(); // print longitudinal profile
 
diff --git a/examples/cascade_proton_example.cpp b/examples/cascade_proton_example.cpp
index 71be6f4f9edc102cbc5d5d190955abe757d07191..5491258ac68006e5e7d483878543b85a3a1d2de8 100644
--- a/examples/cascade_proton_example.cpp
+++ b/examples/cascade_proton_example.cpp
@@ -59,7 +59,6 @@ using namespace std;
 int main() {
 
   logging::set_level(logging::level::info);
-  corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
 
   std::cout << "cascade_proton_example" << std::endl;
 
@@ -147,7 +146,9 @@ int main() {
 
   // define air shower object, run simulation
   Cascade EAS(env, tracking, sequence, output, stack);
+  output.startOfShower();
   EAS.run();
+  output.endOfShower();
 
   cout << "Result: E0=" << E0 / 1_GeV << endl;
   cut.showResults();
diff --git a/examples/corsika.cpp b/examples/corsika.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..a67236dd5cd712cb2c4413f120fd5d0c8297f9fc
--- /dev/null
+++ b/examples/corsika.cpp
@@ -0,0 +1,341 @@
+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * 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.
+ */
+
+#define TRACE
+
+/* clang-format off */
+// InteractionCounter used boost/histogram, which
+// fails if boost/type_traits have been included before. Thus, we have
+// to include it first...
+#include <corsika/framework/process/InteractionCounter.hpp>
+/* clang-format on */
+#include <corsika/framework/geometry/Plane.hpp>
+#include <corsika/framework/geometry/Sphere.hpp>
+#include <corsika/framework/core/Logging.hpp>
+#include <corsika/framework/utility/SaveBoostHistogram.hpp>
+#include <corsika/framework/process/ProcessSequence.hpp>
+#include <corsika/framework/process/SwitchProcessSequence.hpp>
+#include <corsika/framework/process/InteractionCounter.hpp>
+#include <corsika/framework/random/RNGManager.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/utility/CorsikaFenv.hpp>
+#include <corsika/framework/core/Cascade.hpp>
+#include <corsika/framework/geometry/PhysicalGeometry.hpp>
+
+#include <corsika/output/OutputManager.hpp>
+
+#include <corsika/media/Environment.hpp>
+#include <corsika/media/FlatExponential.hpp>
+#include <corsika/media/HomogeneousMedium.hpp>
+#include <corsika/media/IMagneticFieldModel.hpp>
+#include <corsika/media/LayeredSphericalAtmosphereBuilder.hpp>
+#include <corsika/media/NuclearComposition.hpp>
+#include <corsika/media/MediumPropertyModel.hpp>
+#include <corsika/media/UniformMagneticField.hpp>
+#include <corsika/media/ShowerAxis.hpp>
+#include <corsika/media/SlidingPlanarExponential.hpp>
+
+#include <corsika/modules/BetheBlochPDG.hpp>
+#include <corsika/modules/LongitudinalProfile.hpp>
+#include <corsika/modules/ObservationPlane.hpp>
+#include <corsika/modules/OnShellCheck.hpp>
+#include <corsika/modules/StackInspector.hpp>
+#include <corsika/modules/TrackWriter.hpp>
+#include <corsika/modules/ParticleCut.hpp>
+#include <corsika/modules/Pythia8.hpp>
+#include <corsika/modules/Sibyll.hpp>
+#include <corsika/modules/UrQMD.hpp>
+#include <corsika/modules/PROPOSAL.hpp>
+#include <corsika/modules/QGSJetII.hpp>
+
+#include <corsika/setup/SetupStack.hpp>
+#include <corsika/setup/SetupTrajectory.hpp>
+
+#include <iomanip>
+#include <iostream>
+#include <limits>
+#include <string>
+
+/*
+  NOTE, WARNING, ATTENTION
+
+  The .../Random.hpppp implement the hooks of external modules to the C8 random
+  number generator. It has to occur excatly ONCE per linked
+  executable. If you include the header below multiple times and
+  link this togehter, it will fail.
+ */
+#include <corsika/modules/sibyll/Random.hpp>
+#include <corsika/modules/urqmd/Random.hpp>
+#include <corsika/modules/qgsjetII/Random.hpp>
+
+using namespace corsika;
+using namespace std;
+
+using Particle = setup::Stack::particle_type;
+
+void registerRandomStreams(const int seed) {
+  RNGManager::getInstance().registerRandomStream("cascade");
+  RNGManager::getInstance().registerRandomStream("qgsjet");
+  RNGManager::getInstance().registerRandomStream("sibyll");
+  RNGManager::getInstance().registerRandomStream("pythia");
+  RNGManager::getInstance().registerRandomStream("urqmd");
+  RNGManager::getInstance().registerRandomStream("proposal");
+
+  if (seed == 0)
+    RNGManager::getInstance().seedAll();
+  else
+    RNGManager::getInstance().seedAll(seed);
+}
+
+template <typename T>
+using MyExtraEnv = MediumPropertyModel<UniformMagneticField<T>>;
+
+// argv : 1.number of nucleons, 2.number of protons,
+//        3.total energy in GeV, 4.number of showers,
+//        5.seed (0 by default to generate random values for all)
+
+int main(int argc, char** argv) {
+
+  logging::set_level(logging::level::info);
+
+  CORSIKA_LOG_INFO("corsika");
+
+  if (argc < 5) {
+    std::cerr << "usage: corsika <A> <Z> <energy/GeV> <Nevt> [seed] \n"
+                 "       if A=0, Z is interpreted as PDG code \n"
+                 "       if no seed is given, a random seed is chosen \n"
+              << std::endl;
+    return 1;
+  }
+  feenableexcept(FE_INVALID);
+
+  int seed = 0;
+  int number_showers = std::stoi(std::string(argv[4]));
+
+  if (argc > 5) { seed = std::stoi(std::string(argv[5])); }
+
+  // initialize random number sequence(s)
+  registerRandomStreams(seed);
+
+  // setup environment, geometry
+  using EnvType = setup::Environment;
+  EnvType env;
+  CoordinateSystemPtr const& rootCS = env.getCoordinateSystem();
+  Point const center{rootCS, 0_m, 0_m, 0_m};
+  auto builder = make_layered_spherical_atmosphere_builder<
+      setup::EnvironmentInterface, MyExtraEnv>::create(center,
+                                                       constants::EarthRadius::Mean,
+                                                       Medium::AirDry1Atm,
+                                                       MagneticFieldVector{rootCS, 0_T,
+                                                                           50_uT, 0_T});
+  builder.setNuclearComposition(
+      {{Code::Nitrogen, Code::Oxygen},
+       {0.7847f, 1.f - 0.7847f}}); // values taken from AIRES manual, Ar removed for now
+
+  builder.addExponentialLayer(1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 2_km);
+  builder.addExponentialLayer(1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 4_km);
+  builder.addExponentialLayer(1144.9069_g / (1_cm * 1_cm), 878153.55_cm, 10_km);
+  builder.addExponentialLayer(1305.5948_g / (1_cm * 1_cm), 636143.04_cm, 40_km);
+  builder.addExponentialLayer(540.1778_g / (1_cm * 1_cm), 772170.16_cm, 100_km);
+  builder.addLinearLayer(1e9_cm, 112.8_km + constants::EarthRadius::Mean);
+  builder.assemble(env);
+
+  // pre-setup particle stack
+  unsigned short const A = std::stoi(std::string(argv[1]));
+  Code beamCode;
+  HEPEnergyType mass;
+  unsigned short Z = 0;
+  if (A > 0) {
+    beamCode = Code::Nucleus;
+    Z = std::stoi(std::string(argv[2]));
+    mass = get_nucleus_mass(A, Z);
+  } else {
+    int pdg = std::stoi(std::string(argv[2]));
+    beamCode = convert_from_PDG(PDGCode(pdg));
+    mass = get_mass(beamCode);
+  }
+  HEPEnergyType const E0 = 1_GeV * std::stof(std::string(argv[3]));
+  double theta = 20.;
+  double phi = 180.;
+  auto const thetaRad = theta / 180. * M_PI;
+  auto const phiRad = phi / 180. * M_PI;
+
+  auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
+    return sqrt((Elab - m) * (Elab + m));
+  };
+  HEPMomentumType P0 = elab2plab(E0, mass);
+  auto momentumComponents = [](double theta, double phi, HEPMomentumType ptot) {
+    return std::make_tuple(ptot * sin(theta) * cos(phi), ptot * sin(theta) * sin(phi),
+                           -ptot * cos(theta));
+  };
+
+  auto const [px, py, pz] = momentumComponents(thetaRad, phiRad, P0);
+  auto plab = MomentumVector(rootCS, {px, py, pz});
+  cout << "input particle: " << beamCode << endl;
+  cout << "input angles: theta=" << theta << ",phi=" << phi << endl;
+  cout << "input momentum: " << plab.getComponents() / 1_GeV
+       << ", norm = " << plab.getNorm() << endl;
+
+  auto const observationHeight = 0_km + builder.getEarthRadius();
+  auto const injectionHeight = 111.75_km + builder.getEarthRadius();
+  auto const t = -observationHeight * cos(thetaRad) +
+                 sqrt(-static_pow<2>(sin(thetaRad) * observationHeight) +
+                      static_pow<2>(injectionHeight));
+  Point const showerCore{rootCS, 0_m, 0_m, observationHeight};
+  Point const injectionPos =
+      showerCore + DirectionVector{rootCS,
+                                   {-sin(thetaRad) * cos(phiRad),
+                                    -sin(thetaRad) * sin(phiRad), cos(thetaRad)}} *
+                       t;
+
+  std::cout << "point of injection: " << injectionPos.getCoordinates() << std::endl;
+
+  // we make the axis much longer than the inj-core distance since the
+  // profile will go beyond the core, depending on zenith angle
+  std::cout << "shower axis length: " << (showerCore - injectionPos).getNorm() * 1.2
+            << std::endl;
+
+  ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.2, env};
+
+  // create the output manager that we then register outputs with
+  OutputManager output("corsika_outputs");
+
+  // setup processes, decays and interactions
+
+  // corsika::qgsjetII::Interaction qgsjet;
+  corsika::sibyll::Interaction sibyll;
+  InteractionCounter sibyllCounted(sibyll);
+
+  corsika::sibyll::NuclearInteraction sibyllNuc(sibyll, env);
+  InteractionCounter sibyllNucCounted(sibyllNuc);
+
+  corsika::pythia8::Decay decayPythia;
+
+  // use sibyll decay routine for decays of particles unknown to pythia
+  corsika::sibyll::Decay decaySibyll{{
+      Code::N1440Plus,
+      Code::N1440MinusBar,
+      Code::N1440_0,
+      Code::N1440_0Bar,
+      Code::N1710Plus,
+      Code::N1710MinusBar,
+      Code::N1710_0,
+      Code::N1710_0Bar,
+
+      Code::Pi1300Plus,
+      Code::Pi1300Minus,
+      Code::Pi1300_0,
+
+      Code::KStar0_1430_0,
+      Code::KStar0_1430_0Bar,
+      Code::KStar0_1430_Plus,
+      Code::KStar0_1430_MinusBar,
+  }};
+
+  decaySibyll.printDecayConfig();
+
+  ParticleCut cut{50_GeV, 50_GeV, 50_GeV, 50_GeV, false};
+  corsika::proposal::Interaction emCascade(env);
+  corsika::proposal::ContinuousProcess emContinuous(env);
+  InteractionCounter emCascadeCounted(emCascade);
+
+  LongitudinalProfile longprof{showerAxis};
+
+  corsika::urqmd::UrQMD urqmd;
+  InteractionCounter urqmdCounted{urqmd};
+  StackInspector<setup::Stack> stackInspect(5000, false, E0);
+
+  // assemble all processes into an ordered process list
+  struct EnergySwitch {
+    HEPEnergyType cutE_;
+    EnergySwitch(HEPEnergyType cutE)
+        : cutE_(cutE) {}
+    bool operator()(const Particle& p) { return (p.getKineticEnergy() < cutE_); }
+  };
+  auto hadronSequence = make_select(EnergySwitch(80_GeV), urqmdCounted,
+                                    make_sequence(sibyllNucCounted, sibyllCounted));
+  auto decaySequence = make_sequence(decayPythia, decaySibyll);
+
+  TrackWriter trackWriter;
+  output.add("tracks", trackWriter); // register TrackWriter
+
+  Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.}));
+  ObservationPlane observationLevel(obsPlane, DirectionVector(rootCS, {1., 0., 0.}));
+  // register the observation plane with the output
+  output.add("particles", observationLevel);
+
+  auto sequence =
+      make_sequence(stackInspect, hadronSequence, decaySequence, emCascadeCounted,
+                    emContinuous, cut, trackWriter, observationLevel, longprof);
+
+  // define air shower object, run simulation
+  setup::Tracking tracking;
+  setup::Stack stack;
+  Cascade EAS(env, tracking, sequence, output, stack);
+
+  output.startOfLibrary();
+
+  for (int i_shower = 1; i_shower < number_showers + 1; i_shower++) {
+
+    output.startOfShower();
+
+    // directory for outputs
+    string const labHist_file = "inthist_lab_verticalEAS_" + to_string(i_shower) + ".npz";
+    string const cMSHist_file = "inthist_cms_verticalEAS_" + to_string(i_shower) + ".npz";
+    string const longprof_file = "longprof_verticalEAS_" + to_string(i_shower) + ".txt";
+
+    CORSIKA_LOG_INFO("Shower {} / {} ", i_shower, number_showers);
+
+    // setup particle stack, and add primary particle
+    stack.clear();
+
+    if (A > 1) {
+      stack.addParticle(std::make_tuple(beamCode, plab, injectionPos, 0_ns, A, Z));
+
+    } else {
+      if (A == 1) {
+        if (Z == 1) {
+          stack.addParticle(std::make_tuple(Code::Proton, plab, injectionPos, 0_ns));
+        } else if (Z == 0) {
+          stack.addParticle(std::make_tuple(Code::Neutron, plab, injectionPos, 0_ns));
+        } else {
+          std::cerr << "illegal parameters" << std::endl;
+          return EXIT_FAILURE;
+        }
+      } else {
+        stack.addParticle(std::make_tuple(beamCode, plab, injectionPos, 0_ns));
+      }
+    }
+
+    // to fix the point of first interaction, uncomment the following two lines:
+    //  EAS.forceInteraction();
+
+    EAS.run();
+
+    cut.showResults();
+    // emContinuous.showResults();
+    observationLevel.showResults();
+    const HEPEnergyType Efinal = cut.getCutEnergy() + cut.getInvEnergy() +
+                                 cut.getEmEnergy() + // emContinuous.getEnergyLost() +
+                                 observationLevel.getEnergyGround();
+    cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl
+         << "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
+    observationLevel.reset();
+    cut.reset();
+    // emContinuous.reset();
+
+    auto const hists = sibyllCounted.getHistogram() + sibyllNucCounted.getHistogram() +
+                       urqmdCounted.getHistogram();
+
+    save_hist(hists.labHist(), labHist_file, true);
+    save_hist(hists.CMSHist(), cMSHist_file, true);
+    longprof.save(longprof_file);
+    output.endOfShower();
+  }
+  output.endOfLibrary();
+}
diff --git a/examples/em_shower.cpp b/examples/em_shower.cpp
index 12400723bf9c29550f8c424c9b55d117c8b7a82c..9e8d2a39e6740bf8e6403f75be20b9ff80b9741b 100644
--- a/examples/em_shower.cpp
+++ b/examples/em_shower.cpp
@@ -139,7 +139,8 @@ int main(int argc, char** argv) {
             << std::endl;
 
   OutputManager output("em_shower_outputs");
-  ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.02, env};
+  ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.02, env,
+                              false, 1000};
 
   // setup processes, decays and interactions
 
@@ -169,7 +170,9 @@ int main(int argc, char** argv) {
   //  EAS.setNodes();
   //  EAS.forceInteraction();
 
+  output.startOfShower();
   EAS.run();
+  output.endOfShower();
 
   cut.showResults();
   emContinuous.showResults();
diff --git a/examples/hybrid_MC.cpp b/examples/hybrid_MC.cpp
index 25603ad8af4c3f797d5cde7ae47d0190cf960a29..dc4f7e8ee784fa1e5188785a38fa3b6dbc3d4be1 100644
--- a/examples/hybrid_MC.cpp
+++ b/examples/hybrid_MC.cpp
@@ -173,7 +173,8 @@ int main(int argc, char** argv) {
             << std::endl;
 
   OutputManager output("hybrid_MC_outputs");
-  ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.02, env};
+  ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.02, env, true,
+                              1000};
 
   // setup processes, decays and interactions
 
@@ -249,7 +250,9 @@ int main(int argc, char** argv) {
   //  EAS.SetNodes();
   //  EAS.forceInteraction();
 
+  output.startOfShower();
   EAS.run();
+  output.endOfShower();
 
   cut.showResults();
   eLoss.showResults();
diff --git a/examples/stack_example.cpp b/examples/stack_example.cpp
index 41beda4f9c5c6f220411f675fe8e3d8c76e7eda3..da97201afad3001f6a4b9e5cec6f2701df227091 100644
--- a/examples/stack_example.cpp
+++ b/examples/stack_example.cpp
@@ -44,11 +44,11 @@ void read(VectorStack& s) {
 int main() {
 
   logging::set_level(logging::level::info);
-  corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
 
-  std::cout << "stack_example" << std::endl;
+  CORSIKA_LOG_INFO("stack_example");
   VectorStack s;
   fill(s);
   read(s);
+  CORSIKA_LOG_INFO("done");
   return 0;
 }
diff --git a/examples/stopping_power.cpp b/examples/stopping_power.cpp
index 4c9e43efc1883e933433aa83240896442491bb9b..e971a279d27edd940098b228d335efdca9d27594 100644
--- a/examples/stopping_power.cpp
+++ b/examples/stopping_power.cpp
@@ -30,7 +30,6 @@ using namespace std;
 int main() {
 
   logging::set_level(logging::level::info);
-  corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
 
   CORSIKA_LOG_INFO("stopping_power");
 
@@ -48,15 +47,16 @@ int main() {
       rootCS, 0_m, 0_m,
       112.8_km); // this is the CORSIKA 7 start of atmosphere/universe
 
-  ShowerAxis showerAxis{injectionPos, Vector<length_d>{rootCS, 0_m, 0_m, 1_m}, env};
+  ShowerAxis showerAxis{injectionPos, Vector<length_d>{rootCS, 0_m, 0_m, 1_m}, env, false,
+                        100};
   BetheBlochPDG eLoss{showerAxis};
 
   setup::Stack stack;
 
   std::ofstream file("dEdX.dat");
-  file << "# beta*gamma, dE/dX / eV/(g/cm²)" << std::endl;
+  file << "# beta*gamma, dE/dX / MeV/(g/cm²)" << std::endl;
 
-  for (HEPEnergyType E0 = 300_MeV; E0 < 1_PeV; E0 *= 1.05) {
+  for (HEPEnergyType E0 = 200_MeV; E0 < 1_PeV; E0 *= 1.05) {
     stack.clear();
     const Code beamCode = Code::MuPlus;
     const HEPMassType mass = get_mass(beamCode);
@@ -72,16 +72,14 @@ int main() {
                              -ptot * cos(theta));
     };
     auto const [px, py, pz] =
-        momentumComponents(theta / 180. * M_PI, phi / 180. * M_PI, P0);
+        momentumComponents(theta / 180. * constants::pi, phi / 180. * constants::pi, P0);
     auto plab = MomentumVector(rootCS, {px, py, pz});
-    cout << "input particle: " << beamCode << endl;
-    cout << "input angles: theta=" << theta << " phi=" << phi << endl;
-    cout << "input momentum: " << plab.getComponents() / 1_GeV << endl;
 
     stack.addParticle(std::make_tuple(beamCode, plab, injectionPos, 0_ns));
 
     auto const p = stack.getNextParticle();
     HEPEnergyType dE = eLoss.getTotalEnergyLoss(p, 1_g / square(1_cm));
-    file << P0 / mass << "\t" << -dE / 1_eV << std::endl;
+    file << P0 / mass << "\t" << -dE / 1_MeV << std::endl;
   }
+  CORSIKA_LOG_INFO("finished writing dEdX.dat");
 }
diff --git a/examples/vertical_EAS.cpp b/examples/vertical_EAS.cpp
index 21354adcfc31d8fec2ec28381a7aeb4ceb39e39e..b8e16ab432e45620a2007888c20a4a633848e95f 100644
--- a/examples/vertical_EAS.cpp
+++ b/examples/vertical_EAS.cpp
@@ -43,15 +43,12 @@
 #include <corsika/modules/BetheBlochPDG.hpp>
 #include <corsika/modules/LongitudinalProfile.hpp>
 #include <corsika/modules/ObservationPlane.hpp>
-#include <corsika/modules/OnShellCheck.hpp>
 #include <corsika/modules/StackInspector.hpp>
 #include <corsika/modules/TrackWriter.hpp>
 #include <corsika/modules/ParticleCut.hpp>
 #include <corsika/modules/Pythia8.hpp>
 #include <corsika/modules/Sibyll.hpp>
 #include <corsika/modules/UrQMD.hpp>
-#include <corsika/modules/PROPOSAL.hpp>
-#include <corsika/modules/QGSJetII.hpp>
 
 #include <corsika/setup/SetupStack.hpp>
 #include <corsika/setup/SetupTrajectory.hpp>
@@ -71,7 +68,6 @@
  */
 #include <corsika/modules/sibyll/Random.hpp>
 #include <corsika/modules/urqmd/Random.hpp>
-#include <corsika/modules/qgsjetII/Random.hpp>
 
 using namespace corsika;
 using namespace std;
@@ -80,7 +76,6 @@ using Particle = setup::Stack::particle_type;
 
 void registerRandomStreams(const int seed) {
   RNGManager::getInstance().registerRandomStream("cascade");
-  RNGManager::getInstance().registerRandomStream("qgsjet");
   RNGManager::getInstance().registerRandomStream("sibyll");
   RNGManager::getInstance().registerRandomStream("pythia");
   RNGManager::getInstance().registerRandomStream("urqmd");
@@ -105,8 +100,8 @@ int main(int argc, char** argv) {
 
   CORSIKA_LOG_INFO("vertical_EAS");
 
-  if (argc < 5) {
-    std::cerr << "usage: vertical_EAS <A> <Z> <energy/GeV> <Nevt> [seed] \n"
+  if (argc < 4) {
+    std::cerr << "usage: vertical_EAS <A> <Z> <energy/GeV> [seed] \n"
                  "       if A=0, Z is interpreted as PDG code \n"
                  "       if no seed is given, a random seed is chosen \n"
               << std::endl;
@@ -115,9 +110,8 @@ int main(int argc, char** argv) {
   feenableexcept(FE_INVALID);
 
   int seed = 0;
-  int number_showers = std::stoi(std::string(argv[4]));
 
-  if (argc > 5) { seed = std::stoi(std::string(argv[5])); }
+  if (argc > 4) { seed = std::stoi(std::string(argv[4])); }
 
   // initialize random number sequence(s)
   registerRandomStreams(seed);
@@ -162,8 +156,8 @@ int main(int argc, char** argv) {
   HEPEnergyType const E0 = 1_GeV * std::stof(std::string(argv[3]));
   double theta = 0.;
   double phi = 180.;
-  auto const thetaRad = theta / 180. * M_PI;
-  auto const phiRad = phi / 180. * M_PI;
+  auto const thetaRad = theta / 180. * constants::pi;
+  auto const phiRad = phi / 180. * constants::pi;
 
   auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
     return sqrt((Elab - m) * (Elab + m));
@@ -200,14 +194,14 @@ int main(int argc, char** argv) {
   std::cout << "shower axis length: " << (showerCore - injectionPos).getNorm() * 1.5
             << std::endl;
 
-  ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.5, env};
+  ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.5, env, false,
+                              1000};
 
   // create the output manager that we then register outputs with
   OutputManager output("vertical_EAS_outputs");
 
   // setup processes, decays and interactions
 
-  // corsika::qgsjetII::Interaction qgsjet;
   corsika::sibyll::Interaction sibyll;
   InteractionCounter sibyllCounted(sibyll);
 
@@ -240,15 +234,6 @@ int main(int argc, char** argv) {
   decaySibyll.printDecayConfig();
 
   ParticleCut cut{60_GeV, 60_GeV, 60_GeV, 60_GeV, true};
-  corsika::proposal::Interaction emCascade(env);
-  corsika::proposal::ContinuousProcess emContinuous(env);
-  InteractionCounter emCascadeCounted(emCascade);
-
-  OnShellCheck reset_particle_mass(1.e-3, 1.e-1, false);
-
-  LongitudinalProfile longprof{showerAxis};
-
-  Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.}));
 
   corsika::urqmd::UrQMD urqmd;
   InteractionCounter urqmdCounted{urqmd};
@@ -265,170 +250,73 @@ int main(int argc, char** argv) {
                                     make_sequence(sibyllNucCounted, sibyllCounted));
   auto decaySequence = make_sequence(decayPythia, decaySibyll);
 
-  for (int i_shower = 1; i_shower < number_showers + 1; i_shower++) {
-
-    // directory for outputs
-    string const labHist_file = "inthist_lab_verticalEAS_" + to_string(i_shower) + ".npz";
-    string const cMSHist_file = "inthist_cms_verticalEAS_" + to_string(i_shower) + ".npz";
-    string const longprof_file = "longprof_verticalEAS_" + to_string(i_shower) + ".txt";
-
-    std::cout << std::endl;
-    std::cout << "Shower " << i_shower << "/" << number_showers << std::endl;
-
-    // setup particle stack, and add primary particle
-    setup::Stack stack;
-    stack.clear();
-    unsigned short const A = std::stoi(std::string(argv[1]));
-    Code beamCode;
-    HEPEnergyType mass;
-    unsigned short Z = 0;
-    if (A > 0) {
-      beamCode = Code::Nucleus;
-      Z = std::stoi(std::string(argv[2]));
-      mass = get_nucleus_mass(A, Z);
-    } else {
-      int pdg = std::stoi(std::string(argv[2]));
-      beamCode = convert_from_PDG(PDGCode(pdg));
-      mass = get_mass(beamCode);
-    }
-    HEPEnergyType const E0 = 1_GeV * std::stof(std::string(argv[3]));
-    double theta = 0.;
-    auto const thetaRad = theta / 180. * M_PI;
-
-    auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
-      return sqrt((Elab - m) * (Elab + m));
-    };
-    HEPMomentumType P0 = elab2plab(E0, mass);
-    auto momentumComponents = [](double thetaRad, HEPMomentumType ptot) {
-      return std::make_tuple(ptot * sin(thetaRad), 0_eV, -ptot * cos(thetaRad));
-    };
-
-    auto const [px, py, pz] = momentumComponents(thetaRad, P0);
-    auto plab = MomentumVector(rootCS, {px, py, pz});
-    cout << "input particle: " << beamCode << endl;
-    cout << "input angles: theta=" << theta << endl;
-    cout << "input momentum: " << plab.getComponents() / 1_GeV
-         << ", norm = " << plab.getNorm() << endl;
-
-    auto const observationHeight = 0_km + builder.getEarthRadius();
-    auto const injectionHeight = 111.75_km + builder.getEarthRadius();
-    auto const t = (injectionHeight - observationHeight) / cos(thetaRad);
-    Point const showerCore{rootCS, 0_m, 0_m, observationHeight};
-    Point const injectionPos =
-        showerCore + DirectionVector{rootCS, {-sin(thetaRad), 0, cos(thetaRad)}} * t;
-
-    std::cout << "point of injection: " << injectionPos.getCoordinates() << std::endl;
-
-    if (A > 1) {
-      stack.addParticle(std::make_tuple(beamCode, plab, injectionPos, 0_ns, A, Z));
+  // directory for outputs
+  string const labHist_file = "inthist_lab_verticalEAS.npz";
+  string const cMSHist_file = "inthist_cms_verticalEAS.npz";
+  string const longprof_file = "longprof_verticalEAS.txt";
 
-    } else {
-      if (A == 1) {
-        if (Z == 1) {
-          stack.addParticle(std::make_tuple(Code::Proton, plab, injectionPos, 0_ns));
-        } else if (Z == 0) {
-          stack.addParticle(std::make_tuple(Code::Neutron, plab, injectionPos, 0_ns));
-        } else {
-          std::cerr << "illegal parameters" << std::endl;
-          return EXIT_FAILURE;
-        }
+  // setup particle stack, and add primary particle
+  setup::Stack stack;
+  stack.clear();
+
+  if (A > 1) {
+    stack.addParticle(std::make_tuple(beamCode, plab, injectionPos, 0_ns, A, Z));
+
+  } else {
+    if (A == 1) {
+      if (Z == 1) {
+        stack.addParticle(std::make_tuple(Code::Proton, plab, injectionPos, 0_ns));
+      } else if (Z == 0) {
+        stack.addParticle(std::make_tuple(Code::Neutron, plab, injectionPos, 0_ns));
       } else {
-        stack.addParticle(std::make_tuple(beamCode, plab, injectionPos, 0_ns));
+        std::cerr << "illegal parameters" << std::endl;
+        return EXIT_FAILURE;
       }
+    } else {
+      stack.addParticle(std::make_tuple(beamCode, plab, injectionPos, 0_ns));
     }
+  }
 
-    // we make the axis much longer than the inj-core distance since the
-    // profile will go beyond the core, depending on zenith angle
-    std::cout << "shower axis length: " << (showerCore - injectionPos).getNorm() * 1.5
-              << std::endl;
+  BetheBlochPDG emContinuous(showerAxis);
 
-    ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.5, env,
-                                false};
-
-    // setup processes, decays and interactions
-
-    // corsika::qgsjetII::Interaction qgsjet;
-    corsika::sibyll::Interaction sibyll;
-    InteractionCounter sibyllCounted(sibyll);
-
-    corsika::sibyll::NuclearInteraction sibyllNuc(sibyll, env);
-    InteractionCounter sibyllNucCounted(sibyllNuc);
-
-    corsika::pythia8::Decay decayPythia;
-
-    // use sibyll decay routine for decays of particles unknown to pythia
-    corsika::sibyll::Decay decaySibyll{{
-        Code::N1440Plus,
-        Code::N1440MinusBar,
-        Code::N1440_0,
-        Code::N1440_0Bar,
-        Code::N1710Plus,
-        Code::N1710MinusBar,
-        Code::N1710_0,
-        Code::N1710_0Bar,
-
-        Code::Pi1300Plus,
-        Code::Pi1300Minus,
-        Code::Pi1300_0,
-
-        Code::KStar0_1430_0,
-        Code::KStar0_1430_0Bar,
-        Code::KStar0_1430_Plus,
-        Code::KStar0_1430_MinusBar,
-    }};
-
-    decaySibyll.printDecayConfig();
-
-    ParticleCut cut{60_GeV, true, true};
-    // corsika::proposal::Interaction emCascade(env);
-    // corsika::proposal::ContinuousProcess emContinuous(env);
-    // InteractionCounter emCascadeCounted(emCascade);
-    BetheBlochPDG emContinuous(showerAxis);
-
-    OnShellCheck reset_particle_mass(1.e-3, 1.e-1, false);
-    TrackWriter trackWriter;
-    output.add("tracks", trackWriter); // register TrackWriter
-
-    LongitudinalProfile longprof{showerAxis};
-
-    Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.}));
-    ObservationPlane observationLevel(obsPlane, DirectionVector(rootCS, {1., 0., 0.}));
-    // register the observation plane with the output
-    output.add("particles", observationLevel);
-
-    auto sequence = make_sequence( // emCascadeCounted,
-        stackInspect, hadronSequence, reset_particle_mass, decaySequence,
-        // emContinuous,
-        BetheBlochPDG(showerAxis), cut, trackWriter, observationLevel, longprof);
-
-    // define air shower object, run simulation
-    setup::Tracking tracking;
-    Cascade EAS(env, tracking, sequence, output, stack);
-
-    // to fix the point of first interaction, uncomment the following two lines:
-    //  EAS.forceInteraction();
-
-    EAS.run();
-
-    cut.showResults();
-    // emContinuous.showResults();
-    observationLevel.showResults();
-    const HEPEnergyType Efinal = cut.getCutEnergy() + cut.getInvEnergy() +
-                                 cut.getEmEnergy() + // emContinuous.getEnergyLost() +
-                                 observationLevel.getEnergyGround();
-    cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl
-         << "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
-    observationLevel.reset();
-    cut.reset();
-    // emContinuous.reset();
-
-    auto const hists = sibyllCounted.getHistogram() + sibyllNucCounted.getHistogram() +
-                       urqmdCounted.getHistogram();
-
-    save_hist(hists.labHist(), labHist_file, true);
-    save_hist(hists.CMSHist(), cMSHist_file, true);
-    longprof.save(longprof_file);
-
-    output.endOfLibrary();
-  }
+  TrackWriter trackWriter;
+  output.add("tracks", trackWriter); // register TrackWriter
+
+  LongitudinalProfile longprof{showerAxis};
+
+  Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.}));
+  ObservationPlane observationLevel(obsPlane, DirectionVector(rootCS, {1., 0., 0.}));
+  // register the observation plane with the output
+  output.add("particles", observationLevel);
+
+  auto sequence = make_sequence(stackInspect, hadronSequence, decaySequence, emContinuous,
+                                cut, trackWriter, observationLevel, longprof);
+
+  // define air shower object, run simulation
+  setup::Tracking tracking;
+  Cascade EAS(env, tracking, sequence, output, stack);
+  output.startOfShower();
+  EAS.run();
+  output.endOfShower();
+
+  cut.showResults();
+  // emContinuous.showResults();
+  observationLevel.showResults();
+  const HEPEnergyType Efinal = cut.getCutEnergy() + cut.getInvEnergy() +
+                               cut.getEmEnergy() + // emContinuous.getEnergyLost() +
+                               observationLevel.getEnergyGround();
+  cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl
+       << "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
+  observationLevel.reset();
+  cut.reset();
+  // emContinuous.reset();
+
+  auto const hists = sibyllCounted.getHistogram() + sibyllNucCounted.getHistogram() +
+                     urqmdCounted.getHistogram();
+
+  save_hist(hists.labHist(), labHist_file, true);
+  save_hist(hists.CMSHist(), cMSHist_file, true);
+  longprof.save(longprof_file);
+
+  output.endOfLibrary();
 }
diff --git a/tests/common/testTestStack.hpp b/tests/common/testTestStack.hpp
index 70f4d5f147dd7f8d2499f89a1b4809057f77edaa..c3b4687ddb28f49f0e5ec25f2de5123b1bfca246 100644
--- a/tests/common/testTestStack.hpp
+++ b/tests/common/testTestStack.hpp
@@ -71,7 +71,7 @@ public:
   */
 
   // default version for particle-creation from input data
-  void setParticleData(const std::tuple<double> v) { setData(std::get<0>(v)); }
+  void setParticleData(std::tuple<double> v) { setData(std::get<0>(v)); }
   void setParticleData(TestParticleInterface<StackIteratorInterface>& /*parent*/,
                        std::tuple<double> v) {
     setData(std::get<0>(v));
diff --git a/tests/framework/testCascade.hpp b/tests/framework/testCascade.hpp
index 8e6b06bd5d5b1f4fabc6b6f6f0f11b40187fdd83..638424f893c9fa1fef7dd957b1d7fc0010285cc2 100644
--- a/tests/framework/testCascade.hpp
+++ b/tests/framework/testCascade.hpp
@@ -30,7 +30,7 @@ using StackWithGeometryInterface =
                                        SetupGeometryDataInterface, StackIter>;
 
 using TestCascadeStack = corsika::CombinedStack<
-    typename corsika::nuclear_stack::ParticleDataStack::stack_implementation_type,
+    typename corsika::nuclear_stack::ParticleDataStack::stack_data_type,
     corsika::node::GeometryData<TestEnvironmentType>, StackWithGeometryInterface>;
 
 /*
@@ -38,7 +38,7 @@ using TestCascadeStack = corsika::CombinedStack<
 */
 #if defined(__clang__)
 using TestCascadeStackView =
-    corsika::SecondaryView<typename TestCascadeStack::stack_implementation_type,
+    corsika::SecondaryView<typename TestCascadeStack::stack_data_type,
                            StackWithGeometryInterface>;
 #elif defined(__GNUC__) || defined(__GNUG__)
 using TestCascadeStackView = corsika::MakeView<TestCascadeStack>::type;
diff --git a/tests/framework/testCombinedStack.cpp b/tests/framework/testCombinedStack.cpp
index de6defc6ad7a0f94fcbcc289a012de38761112dd..7e86bd88b2beb540328e192969b4a6a878fed57d 100644
--- a/tests/framework/testCombinedStack.cpp
+++ b/tests/framework/testCombinedStack.cpp
@@ -291,8 +291,8 @@ using CombinedTestInterfaceType2 =
     corsika::CombinedParticleInterface<StackTest::pi_type, TestParticleInterface3,
                                        TStackIter>;
 
-using StackTest2 = CombinedStack<typename StackTest::stack_implementation_type,
-                                 TestStackData3, CombinedTestInterfaceType2>;
+using StackTest2 = CombinedStack<typename StackTest::stack_data_type, TestStackData3,
+                                 CombinedTestInterfaceType2>;
 
 TEST_CASE("Combined Stack - multi", "[stack]") {
 
@@ -379,12 +379,12 @@ using CombinedTestInterfaceType2 =
     corsika::CombinedParticleInterface<StackTest::pi_type, TestParticleInterface3,
                                        TStackIter>;
 
-using StackTest2 = CombinedStack<typename StackTest::stack_implementation_type,
-                                 TestStackData3, CombinedTestInterfaceType2>;
+using StackTest2 = CombinedStack<typename StackTest::stack_data_type, TestStackData3,
+                                 CombinedTestInterfaceType2>;
 
 #if defined(__clang__)
-using StackTestView = SecondaryView<typename StackTest2::stack_implementation_type,
-                                    CombinedTestInterfaceType2>;
+using StackTestView =
+    SecondaryView<typename StackTest2::stack_data_type, CombinedTestInterfaceType2>;
 #elif defined(__GNUC__) || defined(__GNUG__)
 using StackTestView = corsika::MakeView<StackTest2>::type;
 #endif
diff --git a/tests/modules/testPythia8.cpp b/tests/modules/testPythia8.cpp
index 846c734877badaaa8cbb508f30da1716ed24b78d..09a7c067263bc4d30e6d477cd7f68e1231203ec1 100644
--- a/tests/modules/testPythia8.cpp
+++ b/tests/modules/testPythia8.cpp
@@ -96,8 +96,10 @@ TEST_CASE("Pythia8Interface", "modules") {
   logging::set_level(logging::level::info);
   auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Proton);
   auto const& cs = *csPtr;
-  [[maybe_unused]] auto const& env_dummy = env;
-  [[maybe_unused]] auto const& node_dummy = nodePtr;
+  {
+    [[maybe_unused]] auto const& env_dummy = env;
+    [[maybe_unused]] auto const& node_dummy = nodePtr;
+  }
 
   SECTION("pythia decay") {
     HEPEnergyType const P0 = 10_GeV;
@@ -228,15 +230,18 @@ TEST_CASE("Pythia8Interface", "modules") {
 
     // incompatible target
     auto [env_Fe, csPtr_Fe, nodePtr_Fe] = setup::testing::setup_environment(Code::Iron);
-    [[maybe_unused]] auto const& cs_Fe = *csPtr_Fe;
-    [[maybe_unused]] auto const& env_dummy_Fe = env_Fe;
-    [[maybe_unused]] auto const& node_dummy_Fe = nodePtr_Fe;
+    {
+      [[maybe_unused]] auto const& cs_Fe = *csPtr_Fe;
+      [[maybe_unused]] auto const& env_dummy_Fe = env_Fe;
+      [[maybe_unused]] auto const& node_dummy_Fe = nodePtr_Fe;
+    }
 
     // resonable projectile, but tool low energy
     auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
         Code::Proton, 0, 0, 1_GeV, (setup::Environment::BaseNodeType* const)nodePtr_Fe,
         *csPtr_Fe);
     auto& view = *secViewPtr;
+    { [[maybe_unused]] auto const& dummy_StackPtr = stackPtr; }
 
     corsika::pythia8::Interaction collision;
 
diff --git a/tests/modules/testQGSJetII.cpp b/tests/modules/testQGSJetII.cpp
index 3d68eb0cd4d14cd07c09bfd9bbc8afc0a2116350..1dfe2f57379d46206ea2d36e7268205b3bcec9c2 100644
--- a/tests/modules/testQGSJetII.cpp
+++ b/tests/modules/testQGSJetII.cpp
@@ -52,15 +52,16 @@ TEST_CASE("CORSIKA_DATA", "[processes]") {
 
   SECTION("check CORSIKA_DATA") {
 
-    const char* data = std::getenv("CORSIKA_DATA");
+    const char* CORSIKA_DATA = std::getenv("CORSIKA_DATA");
     // these CHECKS are needed:
-    CHECK(data != 0);
-    CHECK(boost::filesystem::is_directory(boost::filesystem::path(data) / "QGSJetII"));
+    CHECK(CORSIKA_DATA != 0);
+    CHECK(boost::filesystem::is_directory(boost::filesystem::path(CORSIKA_DATA) /
+                                          "QGSJetII"));
     CORSIKA_LOG_INFO(
         "data: {}"
         " isDir: {}"
         "/QGSJetII",
-        data, boost::filesystem::is_directory(data));
+        CORSIKA_DATA, boost::filesystem::is_directory(CORSIKA_DATA));
   }
 }
 
diff --git a/tests/stack/testGeometryNodeStackExtension.cpp b/tests/stack/testGeometryNodeStackExtension.cpp
index 8256927d92ed44f57d84a0ddfa370244630c9483..e0673f3c33da61acedec35ddb0a1cae365d40044 100644
--- a/tests/stack/testGeometryNodeStackExtension.cpp
+++ b/tests/stack/testGeometryNodeStackExtension.cpp
@@ -34,9 +34,8 @@ using StackWithGeometryInterface =
     CombinedParticleInterface<dummy_stack::DummyStack::pi_type,
                               DummyGeometryDataInterface, TStackIter>;
 
-using TestStack =
-    CombinedStack<typename dummy_stack::DummyStack::stack_implementation_type,
-                  node::GeometryData<DummyEnv>, StackWithGeometryInterface>;
+using TestStack = CombinedStack<typename dummy_stack::DummyStack::stack_data_type,
+                                node::GeometryData<DummyEnv>, StackWithGeometryInterface>;
 
 TEST_CASE("GeometryNodeStackExtension", "[stack]") {
 
diff --git a/tests/stack/testHistoryStack.cpp b/tests/stack/testHistoryStack.cpp
index b393291bfe3598e3a65f848864d62f7610ba7a71..485b794a531a25a0b0a55c7ec637a40432745872 100644
--- a/tests/stack/testHistoryStack.cpp
+++ b/tests/stack/testHistoryStack.cpp
@@ -41,7 +41,7 @@ using StackWithHistoryInterface =
                               TStackIter>;
 
 using TestStack =
-    CombinedStack<typename dummy_stack::DummyStack::stack_implementation_type,
+    CombinedStack<typename dummy_stack::DummyStack::stack_data_type,
                   history::HistoryData<DummyEvent>, StackWithHistoryInterface>;
 
 using EvtPtr = std::shared_ptr<DummyEvent>;
diff --git a/tests/stack/testHistoryView.cpp b/tests/stack/testHistoryView.cpp
index 56e11820979f3b7b0bf9ee2dd72a473118cf7d16..1633b8212c5c59903e99617fd4b2851f8ab133f7 100644
--- a/tests/stack/testHistoryView.cpp
+++ b/tests/stack/testHistoryView.cpp
@@ -32,7 +32,7 @@ using StackWithHistoryInterface =
                               history::HistoryEventDataInterface, TStackIter>;
 
 using TestStack =
-    CombinedStack<typename nuclear_stack::ParticleDataStack::stack_implementation_type,
+    CombinedStack<typename nuclear_stack::ParticleDataStack::stack_data_type,
                   history::HistoryEventData, StackWithHistoryInterface>;
 
 /*
@@ -48,8 +48,8 @@ using TestStack =
   */
 #if defined(__clang__)
 using TheTestStackView =
-    SecondaryView<typename TestStack::stack_implementation_type,
-                  StackWithHistoryInterface, history::HistorySecondaryProducer>;
+    SecondaryView<typename TestStack::stack_data_type, StackWithHistoryInterface,
+                  history::HistorySecondaryProducer>;
 #elif defined(__GNUC__) || defined(__GNUG__)
 using TheTestStackView = MakeView<TestStack, history::HistorySecondaryProducer>::type;
 #endif
diff --git a/tests/stack/testWeightStackExtension.cpp b/tests/stack/testWeightStackExtension.cpp
index 2754c367913e336325491945f6570e4bd014899b..d984e1e5cc36588fac2fd36623df996cd924ab44 100644
--- a/tests/stack/testWeightStackExtension.cpp
+++ b/tests/stack/testWeightStackExtension.cpp
@@ -28,9 +28,8 @@ using StackWithGeometryInterface =
     CombinedParticleInterface<dummy_stack::DummyStack::pi_type, DummyWeightDataInterface,
                               TStackIter>;
 
-using TestStack =
-    CombinedStack<typename dummy_stack::DummyStack::stack_implementation_type,
-                  weights::WeightData, StackWithGeometryInterface>;
+using TestStack = CombinedStack<typename dummy_stack::DummyStack::stack_data_type,
+                                weights::WeightData, StackWithGeometryInterface>;
 
 TEST_CASE("WeightStackExtension", "[stack]") {