this mostly resolves it

Documentation/Examples/cascade_example.cc

@@ -36,13 +36,12 @@ using namespace std;
 
 static int fCount = 0;
 
-class ProcessSplit : public corsika::process::BaseProcess<ProcessSplit> {
+class ProcessSplit : public corsika::process::DiscreteProcess<ProcessSplit> {
 public:
   ProcessSplit() {}
 
   template <typename Particle, typename Track>
-  double MinStepLength(Particle& p, Track&) const {
-
+  double GetInteractionLength(Particle& p, Track&) const {
     // beam particles for sibyll : 1, 2, 3 for p, pi, k
     // read from cross section code table
     int kBeam = 1;
@@ -81,17 +80,20 @@ public:
     std::cout << "ProcessSplit: "
               << "interaction length (g/cm2): " << int_length << std::endl;
     // add exponential sampling
-    int a = 0;
-    const double next_step = -int_length * log(s_rndm_(a));
     /*
       what are the units of the output? slant depth or 3space length?
 
     */
-    std::cout << "ProcessSplit: "
-              << "next step (g/cm2): " << next_step << std::endl;
-    return next_step;
+    return int_length;
+    //
+    //int a = 0;
+    //const double next_step = -int_length * log(s_rndm_(a));
+    //std::cout << "ProcessSplit: "
+    //        << "next step (g/cm2): " << next_step << std::endl;
+    //return next_step;
   }
 
+
   template <typename Particle, typename Track, typename Stack>
   EProcessReturn DoContinuous(Particle&, Track&, Stack&) const {
     // corsika::utls::ignore(p);

Framework/Cascade/Cascade.h

@@ -14,11 +14,12 @@
 
 #include <corsika/process/ProcessReturn.h>
 #include <corsika/units/PhysicalUnits.h>
+#include <corsika/setup/SetupTrajectory.h>
+#include <corsika/random/RNGManager.h>
 
 #include <type_traits>
 
-#include <corsika/setup/SetupTrajectory.h>
-
+using namespace corsika;
 using namespace corsika::units::si;
 
 namespace corsika::cascade {
@@ -50,24 +51,53 @@ namespace corsika::cascade {
         }
         // do cascade equations, which can put new particles on Stack,
         // thus, the double loop
-        // DoCascadeEquations(); //
+        // DoCascadeEquations();
       }
     }
 
-    void Step(Particle& particle) {
+    void Step(Particle& particle) {      
+      
       corsika::setup::Trajectory step = fTracking.GetTrack(particle);
-      fProcesseList.MinStepLength(particle, step);
-
+      const double total_inv_lambda = fProcesseList.GetTotalInverseInteractionLength(particle, step);
+      
+      // sample random exponential step length
+      // ....
+      static corsika::random::RNG& rmng =
+	corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");
+      const double sample_step = rmng() / (double)rmng.max();
+      const double next_step = -log(sample_step) / total_inv_lambda;
+
+      const double min_step = fProcesseList.MinStepLength(particle, step);
+      // convert next_step from grammage to length
+      // Environment::GetDistance(step, next_step);
+      // ....
+      
+      // update step with actual min(min_step, next_step)
+      // ....
+      
       /// here the particle is actually moved along the trajectory to new position:
       // std::visit(corsika::setup::ParticleUpdate<Particle>{particle}, step);
       particle.SetPosition(step.GetPosition(1));
-
+      
       corsika::process::EProcessReturn status =
           fProcesseList.DoContinuous(particle, step, fStack);
+      
       if (status == corsika::process::EProcessReturn::eParticleAbsorbed) {
-        fStack.Delete(particle); // TODO: check if this is really needed
+        // fStack.Delete(particle); // TODO: check if this is really needed
       } else {
-        fProcesseList.DoDiscrete(particle, fStack);
+
+	// check if min_step < random_step  ->  skip discrete if rndm>min_step/random_step
+	if (min_step<next_step) {
+	  const double p_next = min_step/next_step;
+	  const double sample_skip = rmng() / (double)rmng.max();
+	  if (sample_skip>p_next) {
+	    return;// corsika::process::EProcessReturn::eOk;
+	  }
+	}
+	const double sample_process = rmng() / (double)rmng.max();
+	double inv_lambda_count = 0;
+	fProcesseList.SelectDiscrete(particle, fStack, total_inv_lambda,
+				     sample_process, inv_lambda_count);
       }
     }
 

Framework/Cascade/testCascade.cc

@@ -17,6 +17,8 @@
 
 #include <corsika/stack/super_stupid/SuperStupidStack.h>
 
+#include <corsika/particles/ParticleProperties.h>
+
 #include <corsika/geometry/Point.h>
 #include <corsika/geometry/RootCoordinateSystem.h>
 #include <corsika/geometry/Vector.h>
@@ -42,15 +44,18 @@ static int fCount = 0;
 class ProcessSplit : public corsika::process::BaseProcess<ProcessSplit> {
 public:
   ProcessSplit() {}
-
-  template <typename Particle>
-  void MinStepLength(Particle&, setup::Trajectory&) const {}
-
-  template <typename Particle, typename Stack>
-  EProcessReturn DoContinuous(Particle&, setup::Trajectory&, Stack&) const {
+  
+  template <typename Particle, typename T>
+  double MinStepLength(Particle&, T&) const { return 1; }
+  
+  template <typename Particle, typename T>
+  double GetInteractionLength(Particle&, T&) const { return 1; }
+  
+  template <typename Particle, typename T, typename Stack>
+  EProcessReturn DoContinuous(Particle&, T&, Stack&) const {
     return EProcessReturn::eOk;
   }
-
+  
   template <typename Particle, typename Stack>
   void DoDiscrete(Particle& p, Stack& s) const {
     EnergyType E = p.GetEnergy();
@@ -60,7 +65,9 @@ public:
     } else {
       p.SetEnergy(E / 2);
       auto pnew = s.NewParticle();
-      // s.Copy(p, pnew);
+      // s.Copy(p, pnew); fix that .... todo
+      pnew.SetPID(p.GetPID());
+      pnew.SetTime(p.GetTime());
       pnew.SetEnergy(E / 2);
       pnew.SetPosition(p.GetPosition());
       pnew.SetMomentum(p.GetMomentum());
@@ -76,27 +83,34 @@ private:
 
 TEST_CASE("Cascade", "[Cascade]") {
 
-  tracking_line::TrackingLine<setup::Stack> tracking;
+  corsika::random::RNGManager& rmng = corsika::random::RNGManager::GetInstance();
+  ;
+  const std::string str_name = "s_rndm";
+  rmng.RegisterRandomStream(str_name);
 
+  tracking_line::TrackingLine<setup::Stack> tracking;
+  
   stack_inspector::StackInspector<setup::Stack> p0(true);
   ProcessSplit p1;
   const auto sequence = p0 + p1;
   setup::Stack stack;
-
+  
   corsika::cascade::Cascade EAS(tracking, sequence, stack);
-
   CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
-
+  
   stack.Clear();
   auto particle = stack.NewParticle();
   EnergyType E0 = 100_GeV;
+  particle.SetPID(particles::Code::Electron);
   particle.SetEnergy(E0);
   particle.SetPosition(Point(rootCS, {0_m, 0_m, 10_km}));
   particle.SetMomentum(corsika::stack::super_stupid::MomentumVector(
       rootCS, {0 * newton * second, 0 * newton * second, -1 * newton * second}));
+  particle.SetTime(0_ns);
   EAS.Init();
   EAS.Run();
-
+  
+  /*
   SECTION("sectionTwo") {
     for (int i = 0; i < 0; ++i) {
       stack.Clear();
@@ -105,8 +119,9 @@ TEST_CASE("Cascade", "[Cascade]") {
       particle.SetEnergy(E0);
       EAS.Init();
       EAS.Run();
-
+      
       // cout << "Result: E0=" << E0 / 1_GeV << "GeV, count=" << p1.GetCount() << endl;
     }
   }
+  */
 }

Framework/ProcessSequence/BaseProcess.h

@@ -29,14 +29,27 @@ namespace corsika::process {
   struct BaseProcess {
     derived& GetRef() { return static_cast<derived&>(*this); }
     const derived& GetRef() const { return static_cast<const derived&>(*this); }
-
+    
     template <typename Particle, typename Stack>
     inline EProcessReturn DoDiscrete(Particle&, Stack&) const; // {}
-
+    
     template <typename Particle, typename Track, typename Stack>
     inline EProcessReturn DoContinuous(Particle&, Track&, Stack&) const; // {}
-  };
 
+    template <typename Particle, typename Track>
+      inline double GetInverseInteractionLength(Particle& p, Track& t) const {
+      return 1./GetRef().GetInteractionLength(p, t);
+    }
+};
+  
+  /*
+  template<template<typename, typename> class T, typename A, typename B>
+    typename BaseProcess< T<A, B> >::is_process_sequence 
+    {
+      static const bool value = true;
+    };
+  */
+  
   /*
   template <typename T>
   struct is_base {

Framework/ProcessSequence/DiscreteProcess.h

@@ -37,6 +37,12 @@ namespace corsika::process {
     // -> enforce derived to implement DoDiscrete...
     template <typename Particle, typename Stack>
     inline EProcessReturn DoDiscrete(Particle&, Stack&) const; // {}
+
+    template <typename Particle, typename Track>
+    inline double GetInverseInteractionLength(Particle& p, Track& t) const {
+      return 1./GetRef().GetInteractionLength(p, t);
+    }
+    
   };
 
 } // namespace corsika::process

Framework/ProcessSequence/ProcessReturn.h

@@ -23,7 +23,8 @@ namespace corsika::process {
   enum class EProcessReturn {
     eOk = 1,
     eParticleAbsorbed = 2,
-  };
+    eInteracted = 3,
+    };
 } // namespace corsika::process
 
 #endif

Framework/ProcessSequence/ProcessSequence.h

@@ -17,6 +17,9 @@
 #include <corsika/process/DiscreteProcess.h>
 #include <corsika/process/ProcessReturn.h>
 
+#include <limits>
+#include <cmath>
+
 //#include <corsika/setup/SetupTrajectory.h>
 // using corsika::setup::Trajectory;
 //#include <variant>
@@ -25,7 +28,7 @@
 
 namespace corsika::process {
 
-  /* namespace detail { */
+  // namespace detail {
 
   /*   /\* template<typename TT1, typename TT2, typename Type = void> *\/ */
   /*   /\*   struct CallHello { *\/ */
@@ -41,7 +44,7 @@ namespace corsika::process {
   /*   /\* 	static void Call(const TT1&, const TT2&) { *\/ */
   /*   /\* 	  std::cout << "special" << std::endl; *\/ */
   /*   /\* 	}	 *\/ */
-  /*   /\*   }; *\/ */
+  /*          }; */
 
   /*   template<typename T1, typename T2, typename Particle, typename Trajectory, typename
    * Stack> //, typename Type = void> */
@@ -121,7 +124,7 @@ namespace corsika::process {
   /*     } */
   /*   }; */
   /*   *\/ */
-  /* } // end namespace detail */
+  //} // end namespace detail
 
   /**
      \class ProcessSequence
@@ -134,8 +137,24 @@ namespace corsika::process {
      https://en.wikipedia.org/wiki/Curiously_recurring_template_pattern
    */
 
+  //template <typename T1, typename T2>
+  //class ProcessSequence : public BaseProcess<ProcessSequence<T1, T2> >;
+
+  // this is a marker to track which BaseProcess is also a ProcessSequence
+    template <typename T>
+    struct is_process_sequence {
+      static const bool value = false;
+    };
+    
+
   template <typename T1, typename T2>
   class ProcessSequence : public BaseProcess<ProcessSequence<T1, T2> > {
+    
+    /*    template <>
+      struct BaseProcess<ProcessSequence<T1, T2> >::is_process_sequence<true> {
+      static const bool value = true;
+      };*/
+    
   public:
     const T1& A;
     const T2& B;
@@ -162,20 +181,38 @@ namespace corsika::process {
     }
 
     template <typename Particle, typename Track>
-    inline void MinStepLength(Particle& p, Track& track) const {
-      A.MinStepLength(p, track);
-      B.MinStepLength(p, track);
+    inline double MinStepLength(Particle& p, Track& track) const {
+      double min_length = std::numeric_limits<double>::infinity();
+      if constexpr (!std::is_base_of<DiscreteProcess<T1>, T1>::value) {
+	  min_length = std::min(min_length, A.MinStepLength(p,track));
+	}
+      if constexpr (!std::is_base_of<DiscreteProcess<T2>, T2>::value) {
+	  min_length = std::min(min_length, B.MinStepLength(p,track));
+	}
+      return min_length;
+    }
+
+    template <typename Particle, typename Track>
+    inline double GetTotalInteractionLength(Particle& p, Track& t) const {
+      return 1. / GetInverseInteractionLength(p, t);
+    }
+
+    template <typename Particle, typename Track>
+    inline double GetTotalInverseInteractionLength(Particle& p, Track& t) const {
+      return GetInverseInteractionLength(p, t);
     }
 
-    /*
     template <typename Particle, typename Track>
-    inline Track Transport(Particle& p, double& length) const {
-      A.Transport(p, length); // todo: maybe check (?) if there is more than one Transport
-                              // process implemented??
-      return B.Transport(
-          p, length); // need to do this also to decide which Track to return!!!!
+    inline double GetInverseInteractionLength(Particle& p, Track& t) const {
+      double tot = 0;
+      if constexpr (!std::is_base_of<ContinuousProcess<T1>, T1>::value) {
+        tot += A.GetInverseInteractionLength(p, t);
+      }
+      if constexpr (!std::is_base_of<ContinuousProcess<T2>, T2>::value) {
+        tot += B.GetInverseInteractionLength(p, t);
+      }
+      return tot;
     }
-    */
 
     template <typename Particle, typename Stack>
     inline EProcessReturn DoDiscrete(Particle& p, Stack& s) const {
@@ -188,6 +225,47 @@ namespace corsika::process {
       return EProcessReturn::eOk;
     }
 
+    template <typename Particle, typename Stack>
+      inline EProcessReturn SelectDiscrete(Particle& p, Stack& s,
+					   const double lambda_inv_tot,
+					   const double rndm_select,
+					   double& lambda_inv_count) const {
+      if constexpr (is_process_sequence<T1>::value) {
+	  // if A is a process sequence --> check inside
+	  const EProcessReturn ret = A.SelectDiscrete(p, s, lambda_inv_count, rndm_select, lambda_inv_count);
+	  // if A did suceed, stop routine
+	  if (ret != EProcessReturn::eOk) {
+	    return ret;
+	  }
+	} else if constexpr (!std::is_base_of<ContinuousProcess<T1>, T1>::value) {
+	  // if this is not a ContinuousProcess --> evaluate probability
+	  lambda_inv_count += A.GetInverseInteractionLength(p, s);
+	  // check if we should execute THIS process and then EXIT
+	  if (rndm_select<lambda_inv_count/lambda_inv_tot) {
+	    A.DoDiscrete(p, s);
+	    return EProcessReturn::eInteracted;
+	  }
+	} // end branch A
+      
+      if constexpr (is_process_sequence<T2>::value) {
+	  // if A is a process sequence --> check inside
+	  const EProcessReturn ret = B.SelectDiscrete(p, s, lambda_inv_count, rndm_select, lambda_inv_count);
+	  // if A did suceed, stop routine
+	  if (ret != EProcessReturn::eOk) {
+	    return ret;
+	  }
+	} else if constexpr (!std::is_base_of<ContinuousProcess<T2>, T2>::value) {
+	  // if this is not a ContinuousProcess --> evaluate probability
+	  lambda_inv_count += B.GetInverseInteractionLength(p, s);
+	  // check if we should execute THIS process and then EXIT
+	  if (rndm_select<lambda_inv_count/lambda_inv_tot) {
+	    B.DoDiscrete(p, s);
+	    return EProcessReturn::eInteracted;
+	  }
+	} // end branch A
+      return EProcessReturn::eOk;
+    }
+
     /// TODO the const_cast is not nice, think about the constness here
     inline void Init() const {
       const_cast<T1*>(&A)->Init();
@@ -289,6 +367,12 @@ namespace corsika::process {
     }
   */
 
+        template<template<typename, typename> class T, typename A, typename B>
+      struct is_process_sequence< T<A,B> > 
+    {
+      static const bool value = true;
+    };
+
 } // namespace corsika::process
 
 #endif

Framework/ProcessSequence/testProcessSequence.cc

@@ -100,6 +100,11 @@ public:
     cout << "Process2::DoDiscrete" << endl;
     return EProcessReturn::eOk;
   }
+  template <typename Particle, typename Track>
+  inline double GetInteractionLength(Particle&, Track&) const {
+    cout << "Process2::GetInteractionLength" << endl;
+    return 3;
+  }
 };
 
 class Process3 : public DiscreteProcess<Process3> {
@@ -118,6 +123,11 @@ public:
     cout << "Process3::DoDiscrete" << endl;
     return EProcessReturn::eOk;
   }
+  template <typename Particle, typename Track>
+  inline double GetInteractionLength(Particle&, Track&) const {
+    cout << "Process3::GetInteractionLength" << endl;
+    return 1.;
+  }
 };
 
 class Process4 : public BaseProcess<Process4> {
@@ -169,6 +179,20 @@ TEST_CASE("Process Sequence", "[Process Sequence]") {
     // REQUIRE_THROWS(sequence_wrong.Init());
   }
 
+  SECTION("interaction length") {
+    ContinuousProcess1 cp1(0);
+    Process2 m2(1);
+    Process3 m3(2);
+
+    DummyStack s;
+    DummyTrajectory t;
+
+    const auto sequence2 = cp1 + m2 + m3;
+    double tot = sequence2.GetTotalInteractionLength(s, t);
+    double tot_inv = sequence2.GetTotalInverseInteractionLength(s, t);
+    cout << "lambda_tot=" << tot << " lambda_tot_inv=" << tot_inv << endl;
+  }
+  
   SECTION("sectionTwo") {
 
     ContinuousProcess1 cp1(0);

Framework/StackInterface/Stack.h

@@ -67,6 +67,9 @@ namespace corsika::stack {
       IncrementSize();
       return StackIterator(*this, GetSize() - 1);
     }
+    void Copy(StackIterator& a, StackIterator& b) {
+      Copy(a.GetIndex(), b.GetIndex());
+    }
     /// delete this particle
     void Delete(StackIterator& p) {
       if (GetSize() == 0) { /*error*/

Processes/StackInspector/StackInspector.cc

@@ -17,6 +17,7 @@
 
 #include <corsika/setup/SetupTrajectory.h>
 
+#include <limits>
 #include <iostream>
 using namespace std;
 
@@ -56,8 +57,8 @@ process::EProcessReturn StackInspector<Stack>::DoContinuous(Particle&, setup::Tr
 }
 
 template <typename Stack>
-void StackInspector<Stack>::MinStepLength(Particle&, setup::Trajectory&) const {
-  // return 0;
+double StackInspector<Stack>::MinStepLength(Particle&, setup::Trajectory&) const {
+  return std::numeric_limits<double>::infinity();
 }
 
 template <typename Stack>

Processes/StackInspector/StackInspector.h

@@ -36,7 +36,7 @@ namespace corsika::process {
       EProcessReturn DoContinuous(Particle&, corsika::setup::Trajectory&, Stack& s) const;
 
       //      template <typename Particle>
-      void MinStepLength(Particle&, corsika::setup::Trajectory&) const;
+      double MinStepLength(Particle&, corsika::setup::Trajectory&) const;
 
     private:
       bool fReport;

Stack/SuperStupidStack/SuperStupidStack.h

@@ -128,7 +128,7 @@ namespace corsika::stack {
       void Swap(const int i1, const int i2) {
         std::swap(fDataPID[i2], fDataPID[i1]);
         std::swap(fDataE[i2], fDataE[i1]);
-        std::swap(fMomentum[i2], fMomentum[i1]); // should be Momentum !!!!
+        std::swap(fMomentum[i2], fMomentum[i1]);
         std::swap(fPosition[i2], fPosition[i1]);
         std::swap(fTime[i2], fTime[i1]);
       }