diff --git a/corsika/detail/framework/core/Cascade.inl b/corsika/detail/framework/core/Cascade.inl
index bb741627f9cdebe233694a7c55ea5b979a47c46b..2a8b32489820ca21f123626b468a90f3a328afc1 100644
--- a/corsika/detail/framework/core/Cascade.inl
+++ b/corsika/detail/framework/core/Cascade.inl
@@ -97,6 +97,19 @@ namespace corsika {
   template <typename TTracking, typename TProcessList, typename TOutput, typename TStack>
   inline void Cascade<TTracking, TProcessList, TOutput, TStack>::forceInteraction() {
     forceInteraction_ = true;
+    if (forceDecay_) {
+      CORSIKA_LOG_ERROR("Cannot set forceInteraction and forceDecay");
+      throw std::runtime_error("Cannot set forceInteraction and forceDecay");
+    }
+  }
+
+  template <typename TTracking, typename TProcessList, typename TOutput, typename TStack>
+  inline void Cascade<TTracking, TProcessList, TOutput, TStack>::forceDecay() {
+    forceDecay_ = true;
+    if (forceInteraction_) {
+      CORSIKA_LOG_ERROR("Cannot set forceDecay and forceInteraction");
+      throw std::runtime_error("Cannot set forceDecay and forceInteraction");
+    }
   }
 
   template <typename TTracking, typename TProcessList, typename TOutput, typename TStack>
@@ -140,6 +153,22 @@ namespace corsika {
       return;
     }
 
+    if (forceDecay_) {
+      CORSIKA_LOG_TRACE("forced decay!");
+      forceDecay_ = false; // just one (first) interaction
+      stack_view_type secondaries(particle);
+      decay(secondaries, sequence_.getInverseLifetime(particle));
+      if (secondaries.getSize() == 1 && secondaries.getProjectile().getPID() ==
+                                            secondaries.getNextParticle().getPID()) {
+        throw std::runtime_error(
+            fmt::format("Particle {} decays into itself!",
+                        get_name(secondaries.getProjectile().getPID())));
+      }
+      sequence_.doSecondaries(secondaries);
+      particle.erase(); // primary particle is done
+      return;
+    }
+
     // calculate interaction length in medium
     GrammageType const total_lambda =
         (composition.getAverageMassNumber() * constants::u) / total_cx_pre;
diff --git a/corsika/framework/core/Cascade.hpp b/corsika/framework/core/Cascade.hpp
index d4b7d9774c7f951dedb33db7df5c36e4daaf9762..87bea31427df998270957159e59d87e4b4a5ccf6 100644
--- a/corsika/framework/core/Cascade.hpp
+++ b/corsika/framework/core/Cascade.hpp
@@ -92,9 +92,18 @@ namespace corsika {
      * Force an interaction of the top particle of the stack at its current position.
      * Note that setNodes() or an equivalent procedure needs to be called first if you
      * want to call forceInteraction() for the primary interaction.
+     * Incompatible with forceDecay()
      */
     void forceInteraction();
 
+    /**
+     * Force an decay of the top particle of the stack at its current position.
+     * Note that setNodes() or an equivalent procedure needs to be called first if you
+     * want to call forceDecay() for the primary interaction.
+     * Incompatible with forceInteraction()
+     */
+    void forceDecay();
+
   private:
     /**
      * The Step function is executed for each particle from the
@@ -122,6 +131,7 @@ namespace corsika {
     TStack& stack_;
     default_prng_type& rng_ = RNGManager<>::getInstance().getRandomStream("cascade");
     bool forceInteraction_;
+    bool forceDecay_;
     unsigned int count_ = 0;
 
     // but this here temporarily. Should go into dedicated file later: