working cascade with tracking

Documentation/Examples/CMakeLists.txt

@@ -18,6 +18,7 @@ add_executable (staticsequence_example staticsequence_example.cc)
 target_link_libraries (staticsequence_example
   CORSIKAprocesssequence
   CORSIKAunits
+  CORSIKAgeometry
   CORSIKAlogging)
 install (TARGETS staticsequence_example DESTINATION share/examples)
 

Documentation/Examples/staticsequence_example.cc

@@ -15,6 +15,11 @@
 
 #include <corsika/process/ProcessSequence.h>
 
+#include <corsika/setup/SetupTrajectory.h> // TODO: try to break this dependence later
+using corsika::setup::Trajectory;
+#include <corsika/units/PhysicalUnits.h> // dito
+using namespace corsika::units::si;
+
 using namespace std;
 using namespace corsika::process;
 
@@ -72,7 +77,6 @@ struct DummyData {
   double p[10];
 };
 struct DummyStack {};
-struct DummyTrajectory {};
 
 void modular() {
 
@@ -84,9 +88,14 @@ void modular() {
   const auto sequence = m1 + m2 + m3 + m4;
 
   DummyData p;
-  DummyTrajectory t;
   DummyStack s;
 
+  auto const root = corsika::geometry::CoordinateSystem::CreateRootCS();
+  corsika::geometry::Point pos(root, {0_m, 0_m, 0_m});
+  corsika::geometry::Vector<SpeedType::dimension_type> vec(root, {1_m/1_s,0_m/1_s,0_m/1_s});
+  corsika::geometry::Line traj(pos, vec);
+  Trajectory t(corsika::geometry::Trajectory<corsika::geometry::Line>(traj, 0_s, 100_ns));
+
   const int n = 100000000;
   for (int i = 0; i < n; ++i) { sequence.DoContinuous(p, t, s); }
 

Framework/Cascade/CMakeLists.txt

@@ -9,10 +9,12 @@ set (
 set (
   CORSIKAcascade_HEADERS
   Cascade.h
+  #TrackingStep.h
   )
 
 #set (
 #  CORSIKAcascade_SOURCES
+#  TrackingStep.cc
 #  Cascade.cc
 #  )
 

Framework/Cascade/Cascade.h

@@ -9,19 +9,21 @@
  * the license.
  */
 
-#ifndef _include_Cascade_h_
-#define _include_Cascade_h_
+#ifndef _include_corsika_cascade_Cascade_h_
+#define _include_corsika_cascade_Cascade_h_
 
 #include <corsika/process/ProcessReturn.h>
 #include <corsika/units/PhysicalUnits.h>
 
+#include <type_traits>
+
 #include <corsika/setup/SetupTrajectory.h>
 
 using namespace corsika::units::si;
 
 namespace corsika::cascade {
 
-  template <typename ProcessList, typename Stack> //, typename Trajectory>
+  template <typename Tracking, typename ProcessList, typename Stack>
   class Cascade {
 
     typedef typename Stack::ParticleType Particle;
@@ -29,32 +31,26 @@ namespace corsika::cascade {
     Cascade() = delete;
 
   public:
-    Cascade(ProcessList& pl, Stack& stack)
-        : fProcesseList(pl)
-        , fStack(stack) {}
+    Cascade(Tracking& tr, ProcessList& pl, Stack& stack)
+        : fTracking(tr)
+        , fProcesseList(pl)
+        , fStack(stack) {
+      // static_assert(std::is_member_function_pointer<decltype(&ProcessList::DoDiscrete)>::value,
+      //"ProcessList has not function DoDiscrete.");
+      // static_assert(std::is_member_function_pointer<decltype(&ProcessList::DoContinuous)>::value,
+      //	    "ProcessList has not function DoContinuous.");
+    }
 
     void Init() {
-      fStack.Init();
+      fTracking.Init();
       fProcesseList.Init();
+      fStack.Init();
     }
 
     void Run() {
       while (!fStack.IsEmpty()) {
         while (!fStack.IsEmpty()) {
           Particle& pNext = *fStack.GetNextParticle();
-          /*
-            EnergyType Emin;
-            typename Stack::StackIterator pNext(fStack, 0);
-            bool first = true;
-            for (typename Stack::StackIterator ip = fStack.begin(); ip != fStack.end();
-            ++ip) {
-            if (first || ip.GetEnergy() < Emin) {
-            first = false;
-            pNext = ip;
-            Emin = pMin.GetEnergy();
-            }
-            }
-          */
           Step(pNext);
         }
         // do cascade equations, which can put new particles on Stack,
@@ -62,27 +58,33 @@ namespace corsika::cascade {
         // DoCascadeEquations(); //
       }
     }
-    
+
     void Step(Particle& particle) {
+      /*
       [[maybe_unused]] double nextStep = fProcesseList.MinStepLength(particle);
       // corsika::utls::ignore(nextStep);
       auto const root = corsika::geometry::CoordinateSystem::CreateRootCS();
       corsika::geometry::Trajectory<corsika::geometry::Line>
           trajectory( // trajectory is not yet used. this is a dummy.
-		     corsika::geometry::Line(corsika::geometry::Point(root, {0_m, 0_m, 0_m}),
-					     corsika::geometry::Vector<corsika::units::si::SpeedType::dimension_type>(
-														      root, 0 * 1_m / second, 0 * 1_m / second, 1 * 1_m / second)),
-		     0_s, 1_s);
+                     corsika::geometry::Line(corsika::geometry::Point(root, {0_m, 0_m,
+      0_m}), corsika::geometry::Vector<corsika::units::si::SpeedType::dimension_type>(
+                                                                                                                      root, 0 * 1_m / second, 0 * 1_m / second, 1 * 1_m / second)),
+                     0_s, 1_s);
+      */
+      // //[[maybe_unused]] double nextStep = fProcesseList.MinStepLength(particle);
+      corsika::setup::Trajectory step = fTracking.GetTrack(particle);
+      fProcesseList.MinStepLength(particle, step);     
       corsika::process::EProcessReturn status =
-          fProcesseList.DoContinuous(particle, trajectory, fStack);
+	fProcesseList.DoContinuous(particle, step, fStack);
       if (status == corsika::process::EProcessReturn::eParticleAbsorbed) {
-        fStack.Delete(particle);
+        fStack.Delete(particle); // TODO: check if this is really needed
       } else {
         fProcesseList.DoDiscrete(particle, fStack);
       }
     }
-
+    
   private:
+    Tracking& fTracking;
     ProcessList& fProcesseList;
     Stack& fStack;
   };

Framework/Cascade/testCascade.cc

@@ -15,6 +15,7 @@
 
 #include <corsika/setup/SetupStack.h>
 #include <corsika/setup/SetupTrajectory.h>
+using corsika::setup::Trajectory;
 
 #define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one
                           // cpp file
@@ -34,11 +35,11 @@ public:
   ProcessSplit() {}
 
   template <typename Particle>
-  double MinStepLength(Particle&) const {
-    return 0;
+  void MinStepLength(Particle&, Trajectory& ) const {
+    //return 0;
   }
 
-  template <typename Particle, typename Trajectory, typename Stack>
+  template <typename Particle, typename Stack>
   EProcessReturn DoContinuous(Particle&, Trajectory&, Stack&) const {
     // corsika::utls::ignore(p);
     return EProcessReturn::eOk;
@@ -58,19 +59,34 @@ public:
 
   void Init() { fCount = 0; }
 
-  int GetCount() { return fCount; }
+  int GetCount() const { return fCount; }
 
 private:
 };
 
+template<typename Stack>
+class NewtonTracking { // naja.. not yet
+  typedef typename Stack::ParticleType Particle;
+public:
+  void Init() {}
+  corsika::setup::Trajectory GetTrack(Particle& p) {
+    corsika::geometry::Vector<SpeedType::dimension_type> v = p.GetDirection();// * 1_m/1_s;
+    corsika::geometry::Line traj(p.GetPosition(), v);
+    return corsika::geometry::Trajectory<corsika::geometry::Line>(traj, 0_s, 100_ns);
+  }
+};
+
+
 TEST_CASE("Cascade", "[Cascade]") {
 
-  stack_inspector::StackInspector<setup::Stack, setup::Trajectory> p0(true);
+  NewtonTracking<setup::Stack> tracking;
+  
+  stack_inspector::StackInspector<setup::Stack> p0(true);
   ProcessSplit p1;
   const auto sequence = p0 + p1;
   setup::Stack stack;
 
-  corsika::cascade::Cascade EAS(sequence, stack);
+  corsika::cascade::Cascade EAS(tracking, sequence, stack);
 
   stack.Clear();
   auto particle = stack.NewParticle();

Framework/Geometry/CMakeLists.txt

@@ -15,7 +15,7 @@ set (
   BaseVector.h
   QuantityVector.h
   Trajectory.h
-  BaseTrajectory.h
+  # BaseTrajectory.h
   )
 
 set (

Framework/Geometry/Helix.h

@@ -31,7 +31,7 @@ namespace corsika::geometry {
    */
 
   class Helix {
-    
+
     using VelocityVec = Vector<corsika::units::si::SpeedType::dimension_type>;
 
     Point const r0;
@@ -58,8 +58,8 @@ namespace corsika::geometry {
 
     auto GetRadius() const { return radius; }
 
-    LengthType DistanceBetween(corsika::units::si::TimeType t1,
-                               corsika::units::si::TimeType t2) const {
+    LengthType GetDistanceBetween(corsika::units::si::TimeType t1,
+                                  corsika::units::si::TimeType t2) const {
       return (vPar + vPerp).norm() * (t2 - t1);
     }
   };

Framework/Geometry/Line.h

@@ -34,9 +34,9 @@ namespace corsika::geometry {
       return r0 + v0 * t;
     }
     
-    LengthType DistanceBetween(corsika::units::si::TimeType t1,
-                               corsika::units::si::TimeType t2) const {
-      assert(t2 >= t1);
+    LengthType GetDistanceBetween(corsika::units::si::TimeType t1,
+				  corsika::units::si::TimeType t2) const {
+      // assert(t2 >= t1);
       return v0.norm() * (t2 - t1);
     }
   };

Framework/Geometry/Trajectory.h

@@ -12,37 +12,43 @@
 #ifndef _include_TRAJECTORY_H
 #define _include_TRAJECTORY_H
 
-#include <corsika/geometry/BaseTrajectory.h>
+//#include <corsika/geometry/BaseTrajectory.h>
 #include <corsika/units/PhysicalUnits.h>
 
 namespace corsika::geometry {
 
   template <typename T>
-  class Trajectory : public BaseTrajectory {
+  class Trajectory : public T { // BaseTrajectory {
 
-    T fTraj;
+    // T fTraj;
+    corsika::units::si::TimeType fTStart, fTEnd;
 
   public:
+    using T::GetPosition;
+    using T::GetDistanceBetween;
+
     Trajectory(T const& theT, corsika::units::si::TimeType pTStart,
                corsika::units::si::TimeType pTEnd)
-        //: T(theT), fTStart(pTStart), fTEnd(pTEnd) {}
-        : BaseTrajectory(pTStart, pTEnd)
-        , fTraj(theT) {}
+        : T(theT)
+        , fTStart(pTStart)
+        , fTEnd(pTEnd) {}
+    //: BaseTrajectory(pTStart, pTEnd)
+    //  , fTraj(theT) {}
 
-    Point GetPosition(corsika::units::si::TimeType t) const {
+    /*Point GetPosition(corsika::units::si::TimeType t) const {
       return fTraj.GetPosition(t + fTStart);
-    }
+      }*/
 
     Point GetPosition(double u) const {
-      return GetPosition(fTEnd * u + fTStart * (1 - u));
+      return T::GetPosition(fTEnd * u + fTStart * (1 - u));
     }
-    
+
+    /*
     LengthType GetDistance(corsika::units::si::TimeType t1,
-			   corsika::units::si::TimeType t2) const {
+                           corsika::units::si::TimeType t2) const {
       return fTraj.DistanceBetween(t1, t2);
     }
-     
-
+    */
   };
 
 } // namespace corsika::geometry

Framework/Geometry/testGeometry.cc

@@ -154,7 +154,7 @@ TEST_CASE("Trajectories") {
     CHECK(line.GetPosition(2_s).GetCoordinates() ==
           base.GetPosition(2_s).GetCoordinates());
 
-    CHECK(base.GetDistance(1_s, 2_s) / 1_m == Approx(1));
+    CHECK(base.GetDistanceBetween(1_s, 2_s) / 1_m == Approx(1));
   }
 
   SECTION("Helix") {
@@ -182,6 +182,6 @@ TEST_CASE("Trajectories") {
     CHECK(helix.GetPosition(1234_s).GetCoordinates() ==
           base.GetPosition(1234_s).GetCoordinates());
 
-    CHECK(base.GetDistance(1_s, 2_s) / 1_m == Approx(5));
+    CHECK(base.GetDistanceBetween(1_s, 2_s) / 1_m == Approx(5));
   }
 }

Framework/ProcessSequence/CMakeLists.txt

@@ -42,6 +42,7 @@ add_executable (
 
 target_link_libraries (
   testProcessSequence
+  CORSIKAgeometry
   CORSIKAprocesssequence
   CORSIKAthirdparty # for catch2
   )

Framework/ProcessSequence/ProcessSequence.h

@@ -17,8 +17,12 @@
 #include <corsika/process/DiscreteProcess.h>
 #include <corsika/process/ProcessReturn.h>
 
+#include <corsika/setup/SetupTrajectory.h>
+
 //#include <type_traits> // still needed ?
 
+using corsika::setup::Trajectory;
+
 namespace corsika::process {
 
   /* namespace detail { */
@@ -143,7 +147,7 @@ namespace corsika::process {
     // example for a trait-based call:
     // void Hello() const  { detail::CallHello<T1,T2>::Call(A, B); }
 
-    template <typename Particle, typename Trajectory, typename Stack>
+    template <typename Particle, typename Stack>
     inline EProcessReturn DoContinuous(Particle& p, Trajectory& t, Stack& s) const {
       EProcessReturn ret = EProcessReturn::eOk;
       if constexpr (!std::is_base_of<DiscreteProcess<T1>, T1>::value) {
@@ -155,9 +159,10 @@ namespace corsika::process {
       return ret;
     }
 
-    template <typename D>
-    inline double MinStepLength(D& d) const {
-      return std::min(A.MinStepLength(d), B.MinStepLength(d));
+    template <typename Particle>
+      inline void MinStepLength(Particle& p, Trajectory& step) const {
+      A.MinStepLength(p, step);
+      B.MinStepLength(p, step);
     }
 
     /*

Framework/ProcessSequence/testProcessSequence.cc

@@ -19,17 +19,26 @@
 
 #include <corsika/process/ProcessSequence.h>
 
+#include <corsika/setup/SetupTrajectory.h> // TODO: maybe try to break this dependency later!
+using corsika::setup::Trajectory;
+#include <corsika/units/PhysicalUnits.h>
+using namespace corsika::units::si;
+
 using namespace std;
 using namespace corsika::process;
 
+
+static const int nData = 10;
+
+
 class ContinuousProcess1 : public ContinuousProcess<ContinuousProcess1> {
 public:
   ContinuousProcess1() {}
   void Init() { cout << "ContinuousProcess1::Init" << endl; }
-  template <typename D, typename T, typename S>
-  inline EProcessReturn DoContinuous(D& d, T&, S&) const {
+  template <typename D, typename S>
+  inline EProcessReturn DoContinuous(D& d, Trajectory&, S&) const {
     cout << "ContinuousProcess1::DoContinuous" << endl;
-    for (int i = 0; i < 10; ++i) d.p[i] += 0.933;
+    for (int i = 0; i < nData; ++i) d.p[i] += 0.933;
     return EProcessReturn::eOk;
   }
 
@@ -44,10 +53,10 @@ class ContinuousProcess2 : public ContinuousProcess<ContinuousProcess2> {
 public:
   ContinuousProcess2() {}
   void Init() { cout << "ContinuousProcess2::Init" << endl; }
-  template <typename D, typename T, typename S>
-  inline EProcessReturn DoContinuous(D& d, T&, S&) const {
+  template <typename D, typename S>
+  inline EProcessReturn DoContinuous(D& d, Trajectory&, S&) const {
     cout << "ContinuousProcess2::DoContinuous" << endl;
-    for (int i = 0; i < 20; ++i) d.p[i] += 0.933;
+    for (int i = 0; i < nData; ++i) d.p[i] += 0.933;
     return EProcessReturn::eOk;
   }
 
@@ -62,9 +71,9 @@ class Process1 : public DiscreteProcess<Process1> {
 public:
   Process1() {}
   void Init() { cout << "Process1::Init" << endl; }
-  template <typename D, typename T, typename S>
-  inline EProcessReturn DoContinuous(D& d, T&, S&) const {
-    for (int i = 0; i < 10; ++i) d.p[i] += 1 + i;
+  template <typename D, typename S>
+  inline EProcessReturn DoContinuous(D& d, Trajectory&, S&) const {
+    for (int i = 0; i < nData; ++i) d.p[i] += 1 + i;
     return EProcessReturn::eOk;
   }
   template <typename Particle, typename Stack>
@@ -78,9 +87,9 @@ class Process2 : public DiscreteProcess<Process2> {
 public:
   Process2() {}
   void Init() { cout << "Process2::Init" << endl; }
-  template <typename D, typename T, typename S>
-  inline EProcessReturn DoContinuous(D& d, T&, S&) const {
-    for (int i = 0; i < 10; ++i) d.p[i] *= 0.7;
+  template <typename D, typename S>
+  inline EProcessReturn DoContinuous(D& d, Trajectory&, S&) const {
+    for (int i = 0; i < nData; ++i) d.p[i] *= 0.7;
     return EProcessReturn::eOk;
   }
   template <typename Particle, typename Stack>
@@ -94,9 +103,9 @@ class Process3 : public DiscreteProcess<Process3> {
 public:
   Process3() {}
   void Init() { cout << "Process3::Init" << endl; }
-  template <typename D, typename T, typename S>
-  inline EProcessReturn DoContinuous(D& d, T&, S&) const {
-    for (int i = 0; i < 10; ++i) d.p[i] += 0.933;
+  template <typename D, typename S>
+  inline EProcessReturn DoContinuous(D& d, Trajectory&, S&) const {
+    for (int i = 0; i < nData; ++i) d.p[i] += 0.933;
     return EProcessReturn::eOk;
   }
   template <typename Particle, typename Stack>
@@ -110,9 +119,9 @@ class Process4 : public BaseProcess<Process4> {
 public:
   Process4() {}
   void Init() { cout << "Process4::Init" << endl; }
-  template <typename D, typename T, typename S>
-  inline EProcessReturn DoContinuous(D& d, T&, S&) const {
-    for (int i = 0; i < 10; ++i) d.p[i] /= 1.2;
+  template <typename D, typename S>
+  inline EProcessReturn DoContinuous(D& d, Trajectory&, S&) const {
+    for (int i = 0; i < nData; ++i) d.p[i] /= 1.2;
     return EProcessReturn::eOk;
   }
   // inline double MinStepLength(D& d) {
@@ -120,10 +129,9 @@ public:
 };
 
 struct DummyData {
-  double p[10] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
+  double p[nData] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
 };
 struct DummyStack {};
-struct DummyTrajectory {};
 
 TEST_CASE("Cascade", "[Cascade]") {
 
@@ -142,12 +150,18 @@ TEST_CASE("Cascade", "[Cascade]") {
     const auto sequence2 = cp1 + m2 + m3 + cp2;
 
     DummyData p;
-    DummyTrajectory t;
     DummyStack s;
-
+    
     cout << "-->init" << endl;
     sequence2.Init();
     cout << "-->docont" << endl;
+
+    //auto const root = corsika::geometry::CoordinateSystem::CreateRootCS();
+    //corsika::geometry::Point pos(root, {0_m, 0_m, 0_m});
+    //corsika::geometry::Vector<SpeedType::dimension_type> vec(root, {1_m/1_s,0_m/1_s,0_m/1_s});
+    //corsika::geometry::Line traj(pos, vec);
+    Trajectory t;//(corsika::geometry::Trajectory<corsika::geometry::Line>(traj, 0_s, 100_ns));
+
     sequence2.DoContinuous(p, t, s);
     cout << "-->dodisc" << endl;
     sequence2.DoDiscrete(p, s);
@@ -155,11 +169,11 @@ TEST_CASE("Cascade", "[Cascade]") {
 
     sequence.Init();
 
-    const int n = 100;
-    cout << "Running loop with n=" << n << endl;
-    for (int i = 0; i < n; ++i) { sequence.DoContinuous(p, t, s); }
-
-    for (int i = 0; i < 10; i++) { cout << "data[" << i << "]=" << p.p[i] << endl; }
+    const int nLoop = 5;
+    cout << "Running loop with n=" << nLoop << endl;
+    for (int i = 0; i < nLoop; ++i) { sequence.DoContinuous(p, t, s); }
+    for (int i = 0; i < nData; i++) { cout << "data[" << i << "]=" << p.p[i] << endl; }
+    cout << "done" << endl;
   }
 
   SECTION("sectionThree") {}

Processes/StackInspector/StackInspector.cc

@@ -14,6 +14,8 @@
 
 #include <corsika/logging/Logger.h>
 
+#include <corsika/setup/SetupTrajectory.h>
+
 #include <iostream>
 using namespace std;
 
@@ -21,18 +23,16 @@ using namespace corsika;
 using namespace corsika::units::si;
 using namespace corsika::process::stack_inspector;
 
-template <typename Stack, typename Trajectory>
-StackInspector<Stack, Trajectory>::StackInspector(const bool aReport)
+template <typename Stack>
+StackInspector<Stack>::StackInspector(const bool aReport)
     : fReport(aReport) {}
 
-template <typename Stack, typename Trajectory>
-StackInspector<Stack, Trajectory>::~StackInspector() {}
-
-template <typename Stack, typename Trajectory>
-process::EProcessReturn StackInspector<Stack, Trajectory>::DoContinuous(Particle&,
-                                                                        Trajectory&,
-                                                                        Stack& s) const {
+template <typename Stack>
+StackInspector<Stack>::~StackInspector() {}
 
+template <typename Stack>
+process::EProcessReturn StackInspector<Stack>::DoContinuous(
+    Particle&, corsika::setup::Trajectory&, Stack& s) const {
   static int countStep = 0;
   if (!fReport) return EProcessReturn::eOk;
   [[maybe_unused]] int i = 0;
@@ -40,22 +40,26 @@ process::EProcessReturn StackInspector<Stack, Trajectory>::DoContinuous(Particle
   for (auto& iterP : s) {
     EnergyType E = iterP.GetEnergy();
     Etot += E;
+    cout << "i=" << setw(5) << fixed << (i++)
+	 << ", id=" << setw(30) << iterP.GetPID()
+         << " E=" << setw(15) << scientific << (E / 1_GeV) << " GeV, "
+         //<< " pos=" << iterP.GetPosition()
+         << endl;
   }
   countStep++;
   cout << countStep << " " << s.GetSize() << " " << Etot / 1_GeV << " " << endl;
   return EProcessReturn::eOk;
 }
 
-template <typename Stack, typename Trajectory>
-double StackInspector<Stack, Trajectory>::MinStepLength(Particle&) const {
-  return 0;
+template <typename Stack>
+void StackInspector<Stack>::MinStepLength(Particle&,
+                                          corsika::setup::Trajectory&) const {
+  // return 0;
 }
 
-template <typename Stack, typename Trajectory>
-void StackInspector<Stack, Trajectory>::Init() {}
+template <typename Stack>
+void StackInspector<Stack>::Init() {}
 
 #include <corsika/setup/SetupStack.h>
-#include <corsika/setup/SetupTrajectory.h>
 
-template class corsika::process::stack_inspector::StackInspector<setup::Stack,
-                                                                 setup::Trajectory>;
+template class corsika::process::stack_inspector::StackInspector<setup::Stack>;

Processes/StackInspector/StackInspector.h

@@ -14,13 +14,19 @@
 
 #include <corsika/process/ContinuousProcess.h>
 
+#include <corsika/setup/SetupTrajectory.h>
+
+//namespace corsika::setup {
+//class BaseTrajectory;
+//}
+
 namespace corsika::process {
 
   namespace stack_inspector {
 
-    template <typename Stack, typename Trajectory>
+    template <typename Stack>
     class StackInspector
-        : public corsika::process::ContinuousProcess<StackInspector<Stack, Trajectory>> {
+        : public corsika::process::ContinuousProcess<StackInspector<Stack>> {
 
       typedef typename Stack::ParticleType Particle;
 
@@ -31,10 +37,10 @@ namespace corsika::process {
       void Init();
 
       // template <typename Particle, typename Trajectory, typename Stack>
-      EProcessReturn DoContinuous(Particle&, Trajectory&, Stack& s) const;
+      EProcessReturn DoContinuous(Particle&, corsika::setup::Trajectory&, Stack& s) const;
 
-      // template <typename Particle>
-      double MinStepLength(Particle&) const;
+      //      template <typename Particle>
+      void MinStepLength(Particle&, corsika::setup::Trajectory&) const;
 
     private:
       bool fReport;

Setup/SetupTrajectory.h

@@ -12,12 +12,22 @@
 #ifndef _corsika_setup_setuptrajectory_h_
 #define _corsika_setup_setuptrajectory_h_
 
+#include <corsika/geometry/Helix.h>
 #include <corsika/geometry/Line.h>
 #include <corsika/geometry/Trajectory.h>
 
+#include <variant>
+
 namespace corsika::setup {
 
-  typedef corsika::geometry::Trajectory<corsika::geometry::Line> Trajectory;
-}
+  using corsika::geometry::Helix;
+  using corsika::geometry::Line;
+
+  typedef std::variant<std::monostate,
+                       corsika::geometry::Trajectory<Line>,
+                       corsika::geometry::Trajectory<Helix>>
+      Trajectory;
+
+} // namespace corsika::setup
 
 #endif

Stack/SuperStupidStack/SuperStupidStack.h

@@ -30,6 +30,7 @@ namespace corsika::stack {
     using corsika::particles::Code;
     using corsika::units::si::EnergyType;
     using corsika::units::si::TimeType;
+    using corsika::units::si::SpeedType;
     using corsika::units::si::second;
     using corsika::units::si::meter;
     using corsika::units::si::joule;
@@ -62,6 +63,12 @@ namespace corsika::stack {
       MomentumVector GetMomentum() const { return GetStackData().GetMomentum(GetIndex()); }
       Point GetPosition() const { return GetStackData().GetPosition(GetIndex()); }
       TimeType GetTime() const { return GetStackData().GetTime(GetIndex()); }
+      
+#warning this does not really work, nor make sense:
+      Vector<SpeedType::dimension_type> GetDirection() const {
+	auto P = GetMomentum();
+	return P/P.norm() * (corsika::units::si::meter/corsika::units::si::second); }
+    
     };
 
     /**