74 trajectories

Documentation/Examples/cascade_example.cc

@@ -40,8 +40,8 @@ class ProcessSplit : public corsika::process::BaseProcess<ProcessSplit> {
 public:
   ProcessSplit() {}
 
-  template <typename Particle>
-  double MinStepLength(Particle& p, setup::Trajectory&) const {
+  template <typename Particle, typename Track>
+  double MinStepLength(Particle& p, Track&) const {
 
     // beam particles for sibyll : 1, 2, 3 for p, pi, k
     // read from cross section code table
@@ -92,8 +92,8 @@ public:
     return next_step;
   }
 
-  template <typename Particle, typename Stack>
-  EProcessReturn DoContinuous(Particle&, Trajectory&, Stack&) const {
+  template <typename Particle, typename Track, typename Stack>
+  EProcessReturn DoContinuous(Particle&, Track&, Stack&) const {
     // corsika::utls::ignore(p);
     return EProcessReturn::eOk;
   }

Documentation/Examples/geometry_example.cc

@@ -50,12 +50,12 @@ int main() {
   std::cout << "p2-p1 components in cs3: " << diff.GetComponents(cs3)
             << std::endl; // but not under rotations
   std::cout << "p2-p1 norm^2: " << norm << std::endl;
-  assert(norm == 1 * meter*meter);
-  
+  assert(norm == 1 * meter * meter);
+
   Sphere s(p1, 10_m); // define a sphere around a point with a radius
   std::cout << "p1 inside s:  " << s.Contains(p2) << std::endl;
   assert(s.Contains(p2) == 1);
-  
+
   Sphere s2(p1, 3_um); // another sphere
   std::cout << "p1 inside s2: " << s2.Contains(p2) << std::endl;
   assert(s2.Contains(p2) == 0);

Documentation/Examples/stack_example.cc

@@ -11,9 +11,9 @@
 
 #include <corsika/particles/ParticleProperties.h>
 #include <corsika/stack/super_stupid/SuperStupidStack.h>
+#include <cassert>
 #include <iomanip>
 #include <iostream>
-#include <cassert>
 
 using namespace corsika::units::si;
 using namespace corsika::stack;
@@ -28,7 +28,7 @@ void fill(corsika::stack::super_stupid::SuperStupidStack& s) {
 }
 
 void read(corsika::stack::super_stupid::SuperStupidStack& s) {
-  assert(s.GetSize() == 11);  // stack has 11 particles
+  assert(s.GetSize() == 11); // stack has 11 particles
 
   EnergyType total_energy;
   int i = 0;
@@ -38,7 +38,7 @@ void read(corsika::stack::super_stupid::SuperStupidStack& s) {
     assert(p.GetPID() == corsika::particles::Code::Electron);
     assert(p.GetEnergy() == 1.5_GeV * (i++));
   }
-  //assert(total_energy == 82.5_GeV);
+  // assert(total_energy == 82.5_GeV);
 }
 
 int main() {

Documentation/Examples/staticsequence_example.cc

@@ -15,20 +15,23 @@
 
 #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;
+#include <corsika/geometry/Point.h>
+#include <corsika/geometry/RootCoordinateSystem.h>
+#include <corsika/geometry/Vector.h>
 
-using namespace std;
+using namespace corsika;
+using namespace corsika::units::si;
 using namespace corsika::process;
+using namespace std;
+
+const int nData = 10;
 
 class Process1 : public BaseProcess<Process1> {
 public:
   Process1() {}
   template <typename D, typename T, typename S>
   EProcessReturn DoContinuous(D& d, T&, S&) const {
-    for (int i = 0; i < 10; ++i) d.p[i] += 1;
+    for (int i = 0; i < nData; ++i) d.p[i] += 1;
     return EProcessReturn::eOk;
   }
 };
@@ -39,66 +42,64 @@ public:
 
   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.1*i;
+    for (int i = 0; i < nData; ++i) d.p[i] -= 0.1 * i;
     return EProcessReturn::eOk;
   }
 };
 
 class Process3 : public BaseProcess<Process3> {
 public:
-  // Process3(const int v) :fV(v) {}
   Process3() {}
 
   template <typename D, typename T, typename S>
-  inline EProcessReturn DoContinuous(D& /*d*/, T& /*t*/, S& /*s*/) const {
-    // for (int i=0; i<10; ++i) d.p[i] += fV;
+  inline EProcessReturn DoContinuous(D&, T&, S&) const {
     return EProcessReturn::eOk;
   }
-
-private:
-  // int fV;
 };
 
 class Process4 : public BaseProcess<Process4> {
 public:
-  // Process4(const int v) : fV(v) {}
-  Process4() {}
+  Process4(const double v)
+      : fV(v) {}
   template <typename D, typename T, typename S>
-  inline EProcessReturn DoContinuous(D& /*d*/, T& /*t*/, S& /*s*/) const {
-    // for (int i=0; i<10; ++i) d.p[i] /= fV;
+  inline EProcessReturn DoContinuous(D& d, T&, S&) const {
+    for (int i = 0; i < nData; ++i) d.p[i] *= fV;
     return EProcessReturn::eOk;
   }
 
 private:
-  // int fV;
+  double fV;
 };
 
 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 {
+  void clear() {}
 };
-struct DummyStack {};
+struct DummyTrajectory {};
 
 void modular() {
 
   Process1 m1;
   Process2 m2;
   Process3 m3;
-  Process4 m4;
+  Process4 m4(0.9);
 
   const auto sequence = m1 + m2 + m3 + m4;
-  
+
   DummyData p;
   DummyStack s;
-  Trajectory t;
+  DummyTrajectory t;
 
   const int n = 1000;
   for (int i = 0; i < n; ++i) { sequence.DoContinuous(p, t, s); }
-  
-  for(int i=0; i<10; ++i) {
-    //cout << p.p[i] << endl;
-    //assert(p.p[i] == n-i*100);
+
+  for (int i = 0; i < nData; ++i) {
+    // cout << p.p[i] << endl;
+    // assert(p.p[i] == n-i*100);
   }
-  
+
   cout << " done (nothing...) " << endl;
 }
 

Framework/Cascade/Cascade.h

@@ -34,12 +34,7 @@ namespace corsika::cascade {
     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.");
-    }
+        , fStack(stack) {}
 
     void Init() {
       fTracking.Init();
@@ -64,7 +59,8 @@ namespace corsika::cascade {
       fProcesseList.MinStepLength(particle, step);
 
       /// here the particle is actually moved along the trajectory to new position:
-      std::visit(corsika::setup::ParticleUpdate<Particle>{particle}, step);
+      // std::visit(corsika::setup::ParticleUpdate<Particle>{particle}, step);
+      particle.SetPosition(step.GetPosition(1));
 
       corsika::process::EProcessReturn status =
           fProcesseList.DoContinuous(particle, step, fStack);

Framework/Cascade/testCascade.cc

@@ -18,6 +18,7 @@
 #include <corsika/stack/super_stupid/SuperStupidStack.h>
 
 #include <corsika/geometry/Point.h>
+#include <corsika/geometry/RootCoordinateSystem.h>
 #include <corsika/geometry/Vector.h>
 
 #include <corsika/setup/SetupStack.h>

Framework/Geometry/Point.h

@@ -26,27 +26,26 @@ namespace corsika::geometry {
    * A Point represents a point in position space. It is defined by its
    * coordinates with respect to some CoordinateSystem.
    */
-  class Point : public BaseVector<phys::units::length_d> {
+  class Point : public BaseVector<length_d> {
   public:
-    Point(CoordinateSystem const& pCS, QuantityVector<phys::units::length_d> pQVector)
-        : BaseVector<phys::units::length_d>(pCS, pQVector) {}
+    Point(CoordinateSystem const& pCS, QuantityVector<length_d> pQVector)
+        : BaseVector<length_d>(pCS, pQVector) {}
 
     Point(CoordinateSystem const& cs, LengthType x, LengthType y, LengthType z)
-        : BaseVector<phys::units::length_d>(cs, {x, y, z}) {}
+        : BaseVector<length_d>(cs, {x, y, z}) {}
 
     // TODO: this should be private or protected, we don NOT want to expose numbers
     // without reference to outside:
-    auto GetCoordinates() const { return BaseVector<phys::units::length_d>::qVector; }
+    auto GetCoordinates() const { return BaseVector<length_d>::qVector; }
 
     /// this always returns a QuantityVector as triple
     auto GetCoordinates(CoordinateSystem const& pCS) const {
-      if (&pCS == BaseVector<phys::units::length_d>::cs) {
-        return BaseVector<phys::units::length_d>::qVector;
+      if (&pCS == BaseVector<length_d>::cs) {
+        return BaseVector<length_d>::qVector;
       } else {
-        return QuantityVector<phys::units::length_d>(
-            CoordinateSystem::GetTransformation(*BaseVector<phys::units::length_d>::cs,
-                                                pCS) *
-            BaseVector<phys::units::length_d>::qVector.eVector);
+        return QuantityVector<length_d>(
+            CoordinateSystem::GetTransformation(*BaseVector<length_d>::cs, pCS) *
+            BaseVector<length_d>::qVector.eVector);
       }
     }
 
@@ -55,22 +54,21 @@ namespace corsika::geometry {
      * coordinates interally
      */
     void rebase(CoordinateSystem const& pCS) {
-      BaseVector<phys::units::length_d>::qVector = GetCoordinates(pCS);
-      BaseVector<phys::units::length_d>::cs = &pCS;
+      BaseVector<length_d>::qVector = GetCoordinates(pCS);
+      BaseVector<length_d>::cs = &pCS;
     }
 
-    Point operator+(Vector<phys::units::length_d> const& pVec) const {
-      return Point(
-          *BaseVector<phys::units::length_d>::cs,
-          GetCoordinates() + pVec.GetComponents(*BaseVector<phys::units::length_d>::cs));
+    Point operator+(Vector<length_d> const& pVec) const {
+      return Point(*BaseVector<length_d>::cs,
+                   GetCoordinates() + pVec.GetComponents(*BaseVector<length_d>::cs));
     }
 
     /*!
      * returns the distance Vector between two points
      */
-    Vector<phys::units::length_d> operator-(Point const& pB) const {
-      auto& cs = *BaseVector<phys::units::length_d>::cs;
-      return Vector<phys::units::length_d>(cs, GetCoordinates() - pB.GetCoordinates(cs));
+    Vector<length_d> operator-(Point const& pB) const {
+      auto& cs = *BaseVector<length_d>::cs;
+      return Vector<length_d>(cs, GetCoordinates() - pB.GetCoordinates(cs));
     }
   };
 

Framework/ProcessSequence/BaseProcess.h

@@ -29,8 +29,15 @@ 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 T>
   struct is_base {
     static const bool value = false;
@@ -39,6 +46,7 @@ namespace corsika::process {
   struct is_base<BaseProcess<T>> {
     static const bool value = true;
   };
+  */
 
 } // namespace corsika::process
 

Framework/ProcessSequence/ContinuousProcess.h

@@ -13,6 +13,7 @@
 #define _include_corsika_continuousprocess_h_
 
 #include <corsika/process/ProcessReturn.h> // for convenience
+//#include <corsika/setup/SetupTrajectory.h>
 
 namespace corsika::process {
 
@@ -29,6 +30,11 @@ namespace corsika::process {
   struct ContinuousProcess {
     derived& GetRef() { return static_cast<derived&>(*this); }
     const derived& GetRef() const { return static_cast<const derived&>(*this); }
+
+    // here starts the interface part
+    // -> enforce derived to implement DoContinuous...
+    template <typename D, typename T, typename S>
+    inline EProcessReturn DoContinuous(D&, T&, S&) const;
   };
 
 } // namespace corsika::process

Framework/ProcessSequence/DiscreteProcess.h

@@ -13,7 +13,7 @@
 #define _include_corsika_discreteprocess_h_
 
 #include <corsika/process/ProcessReturn.h> // for convenience
-
+#include <corsika/setup/SetupTrajectory.h>
 #include <iostream> // debug
 
 namespace corsika::process {
@@ -30,25 +30,13 @@ namespace corsika::process {
   template <typename derived>
   struct DiscreteProcess {
 
-    // DiscreteProcess() {
-    // static_assert(mustProvide<derived>::mustProvide, "");
-    //}
-
     derived& GetRef() { return static_cast<derived&>(*this); }
     const derived& GetRef() const { return static_cast<const derived&>(*this); }
 
-    // here starts the interface part
+    /// here starts the interface-definition part
     // -> enforce derived to implement DoDiscrete...
     template <typename Particle, typename Stack>
     inline EProcessReturn DoDiscrete(Particle&, Stack&) const; // {}
-
-    //  private:
-    template <typename D, typename T, typename S>
-    inline EProcessReturn DoContinuous(D& d, T&, S&) const {
-      std::cout << "yeah" << std::endl;
-      return EProcessReturn::eOk;
-    } // find out how to make this FINAL
-    // void DoContinuous;
   };
 
 } // namespace corsika::process

Framework/ProcessSequence/ProcessSequence.h

@@ -17,14 +17,12 @@
 #include <corsika/process/DiscreteProcess.h>
 #include <corsika/process/ProcessReturn.h>
 
-#include <corsika/setup/SetupTrajectory.h>
-
-#include <variant>
+//#include <corsika/setup/SetupTrajectory.h>
+// using corsika::setup::Trajectory;
+//#include <variant>
 
 //#include <type_traits> // still needed ?
 
-using corsika::setup::Trajectory;
-
 namespace corsika::process {
 
   /* namespace detail { */
@@ -149,31 +147,33 @@ namespace corsika::process {
     // example for a trait-based call:
     // void Hello() const  { detail::CallHello<T1,T2>::Call(A, B); }
 
-    template <typename Particle, typename Stack>
-    inline EProcessReturn DoContinuous(Particle& p, Trajectory& t, Stack& s) const {
+    template <typename Particle, typename Track, typename Stack>
+    inline EProcessReturn DoContinuous(Particle& p, Track& t, Stack& s) const {
       EProcessReturn ret = EProcessReturn::eOk;
       if constexpr (!std::is_base_of<DiscreteProcess<T1>, T1>::value) {
+        // A.DoContinuous(std::forward<Particle>(p), t, std::forward<Stack>(s));
         A.DoContinuous(p, t, s);
       }
       if constexpr (!std::is_base_of<DiscreteProcess<T2>, T2>::value) {
+        // B.DoContinuous(std::forward<Particle>(p), t, std::forward<Stack>(s));
         B.DoContinuous(p, t, s);
       }
       return ret;
     }
 
-    template <typename Particle>
-    inline void MinStepLength(Particle& p, Trajectory& step) const {
-      A.MinStepLength(p, step);
-      B.MinStepLength(p, step);
+    template <typename Particle, typename Track>
+    inline void MinStepLength(Particle& p, Track& track) const {
+      A.MinStepLength(p, track);
+      B.MinStepLength(p, track);
     }
 
     /*
-    template <typename Particle, typename Trajectory>
-    inline Trajectory Transport(Particle& p, double& length) const {
+    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 Trajectory to return!!!!
+          p, length); // need to do this also to decide which Track to return!!!!
     }
     */
 

Framework/ProcessSequence/testProcessSequence.cc

@@ -19,76 +19,81 @@
 
 #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;
 using namespace corsika::units::si;
-
-using namespace std;
 using namespace corsika::process;
+using namespace std;
 
 static const int nData = 10;
 
+int globalCount = 0;
+
 class ContinuousProcess1 : public ContinuousProcess<ContinuousProcess1> {
+  int fV = 0;
+
 public:
-  ContinuousProcess1() {}
-  void Init() { cout << "ContinuousProcess1::Init" << endl; }
-  template <typename D, typename S>
-  inline EProcessReturn DoContinuous(D& d, Trajectory&, S&) const {
+  ContinuousProcess1(const int v)
+      : fV(v) {}
+  void Init() {
+    cout << "ContinuousProcess1::Init" << endl;
+    assert(globalCount == fV);
+    globalCount++;
+  }
+  template <typename D, typename T, typename S>
+  inline EProcessReturn DoContinuous(D& d, T&, S&) const {
     cout << "ContinuousProcess1::DoContinuous" << endl;
     for (int i = 0; i < nData; ++i) d.p[i] += 0.933;
     return EProcessReturn::eOk;
   }
-
-  template <typename Particle, typename Stack>
-  inline EProcessReturn DoDiscrete(Particle&, Stack&) const {
-    cout << "ContinuousProcess1::DoDiscrete" << endl;
-    return EProcessReturn::eOk;
-  }
 };
 
 class ContinuousProcess2 : public ContinuousProcess<ContinuousProcess2> {
+  int fV = 0;
+
 public:
-  ContinuousProcess2() {}
-  void Init() { cout << "ContinuousProcess2::Init" << endl; }
-  template <typename D, typename S>
-  inline EProcessReturn DoContinuous(D& d, Trajectory&, S&) const {
+  ContinuousProcess2(const int v)
+      : fV(v) {}
+  void Init() {
+    cout << "ContinuousProcess2::Init" << endl;
+    assert(globalCount == fV);
+    globalCount++;
+  }
+  template <typename D, typename T, typename S>
+  inline EProcessReturn DoContinuous(D& d, T&, S&) const {
     cout << "ContinuousProcess2::DoContinuous" << endl;
     for (int i = 0; i < nData; ++i) d.p[i] += 0.933;
     return EProcessReturn::eOk;
   }
-
-  template <typename Particle, typename Stack>
-  inline EProcessReturn DoDiscrete(Particle&, Stack&) const {
-    cout << "ContinuousProcess2::DoDiscrete" << endl;
-    return EProcessReturn::eOk;
-  }
 };
 
 class Process1 : public DiscreteProcess<Process1> {
 public:
-  Process1() {}
-  void Init() { cout << "Process1::Init" << endl; }
+  Process1(const int v)
+      : fV(v) {}
+  void Init() {
+    cout << "Process1::Init" << endl;
+    assert(globalCount == fV);
+    globalCount++;
+  }
   template <typename D, typename S>
-  inline EProcessReturn DoContinuous(D& d, Trajectory&, S&) const {
+  inline EProcessReturn DoDiscrete(D& d, S&) const {
     for (int i = 0; i < nData; ++i) d.p[i] += 1 + i;
     return EProcessReturn::eOk;
   }
-  template <typename Particle, typename Stack>
-  inline EProcessReturn DoDiscrete(Particle&, Stack&) const {
-    cout << "Process1::DoDiscrete" << endl;
-    return EProcessReturn::eOk;
-  }
+  // private:
+  int fV;
 };
 
 class Process2 : public DiscreteProcess<Process2> {
+  int fV = 0;
+
 public:
-  Process2() {}
-  void Init() { cout << "Process2::Init" << endl; }
-  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;
+  Process2(const int v)
+      : fV(v) {}
+  void Init() {
+    cout << "Process2::Init" << endl;
+    assert(globalCount == fV);
+    globalCount++;
   }
   template <typename Particle, typename Stack>
   inline EProcessReturn DoDiscrete(Particle&, Stack&) const {
@@ -98,13 +103,15 @@ public:
 };
 
 class Process3 : public DiscreteProcess<Process3> {
+  int fV = 0;
+
 public:
-  Process3() {}
-  void Init() { cout << "Process3::Init" << endl; }
-  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;
+  Process3(const int v)
+      : fV(v) {}
+  void Init() {
+    cout << "Process3::Init" << endl;
+    assert(globalCount == fV);
+    globalCount++;
   }
   template <typename Particle, typename Stack>
   inline EProcessReturn DoDiscrete(Particle&, Stack&) const {
@@ -114,66 +121,84 @@ public:
 };
 
 class Process4 : public BaseProcess<Process4> {
+  int fV = 0;
+
 public:
-  Process4() {}
-  void Init() { cout << "Process4::Init" << endl; }
-  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;
+  Process4(const int v)
+      : fV(v) {}
+  void Init() {
+    cout << "Process4::Init" << endl;
+    assert(globalCount == fV);
+    globalCount++;
+  }
+  template <typename D, typename T, typename S>
+  inline EProcessReturn DoContinuous(D& d, T&, S&) const {
+    for (int i = 0; i < nData; ++i) { d.p[i] /= 1.2; }
     return EProcessReturn::eOk;
   }
   // inline double MinStepLength(D& d) {
-  // void DoDiscrete(Particle& p, Stack& s) const {
+  template <typename Particle, typename Stack>
+  EProcessReturn DoDiscrete(Particle&, Stack&) const {
+    return EProcessReturn::eOk;
+  }
 };
 
 struct DummyData {
   double p[nData] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
 };
 struct DummyStack {};
+struct DummyTrajectory {};
 
-TEST_CASE("Cascade", "[Cascade]") {
+TEST_CASE("Process Sequence", "[Process Sequence]") {
 
-  SECTION("sectionTwo") {
-
-    Process1 m1;
-    Process2 m2;
-    Process3 m3;
-    Process4 m4;
+  SECTION("Check init order") {
+    Process1 m1(0);
+    Process2 m2(1);
+    Process3 m3(2);
+    Process4 m4(3);
 
     const auto sequence = m1 + m2 + m3 + m4;
 
-    ContinuousProcess1 cp1;
-    ContinuousProcess2 cp2;
+    globalCount = 0;
+    sequence.Init();
+    // REQUIRE_NOTHROW( (sequence.Init()) );
+
+    // const auto sequence_wrong = m3 + m2 + m1 + m4;
+    // globalCount = 0;
+    // sequence_wrong.Init();
+    // REQUIRE_THROWS(sequence_wrong.Init());
+  }
+
+  SECTION("sectionTwo") {
+
+    ContinuousProcess1 cp1(0);
+    ContinuousProcess2 cp2(3);
+    Process2 m2(1);
+    Process3 m3(2);
 
     const auto sequence2 = cp1 + m2 + m3 + cp2;
 
     DummyData p;
     DummyStack s;
+    DummyTrajectory t;
 
-    cout << "-->init" << endl;
+    cout << "-->init sequence2" << endl;
+    globalCount = 0;
     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);
     cout << "-->done" << endl;
 
-    sequence.Init();
-
     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 < nLoop; ++i) {
+      sequence2.DoContinuous(p, t, s);
+      sequence2.DoDiscrete(p, s);
+    }
     for (int i = 0; i < nData; i++) { cout << "data[" << i << "]=" << p.p[i] << endl; }
     cout << "done" << endl;
   }
-
-  SECTION("sectionThree") {}
 }

Setup/SetupTrajectory.h

@@ -26,9 +26,12 @@ namespace corsika::setup {
   using corsika::geometry::Line;
 
   /// definition of Trajectory base class, to be used in tracking and cascades
+  typedef corsika::geometry::Trajectory<Line> Trajectory;
+
+  /*
   typedef std::variant<std::monostate, corsika::geometry::Trajectory<Line>,
                        corsika::geometry::Trajectory<Helix>>
-      Trajectory;
+                       Trajectory;
 
   /// helper visitor to modify Particle by moving along Trajectory
   template <typename Particle>
@@ -58,7 +61,7 @@ namespace corsika::setup {
       return trajectory.GetDuration();
     }
   };
-
+  */
 } // namespace corsika::setup
 
 #endif

Stack/SuperStupidStack/SuperStupidStack.h

@@ -89,8 +89,11 @@ namespace corsika::stack {
       void Init() {}
 
       void Clear() {
-        fDataE.clear();
         fDataPID.clear();
+        fDataE.clear();
+        fMomentum.clear();
+        fPosition.clear();
+        fTime.clear();
       }
 
       int GetSize() const { return fDataPID.size(); }