first moving particles in corsika

Framework/Cascade/Cascade.h

@@ -60,20 +60,12 @@ namespace corsika::cascade {
     }
 
     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);
-      */
-      // //[[maybe_unused]] double nextStep = fProcesseList.MinStepLength(particle);
       corsika::setup::Trajectory step = fTracking.GetTrack(particle);
-      fProcesseList.MinStepLength(particle, step);     
+      fProcesseList.MinStepLength(particle, step);
+
+      /// here the particle is actually moved along the trajectory to new position:
+      std::visit(corsika::setup::ParticleUpdate<Particle>{particle}, step);
+
       corsika::process::EProcessReturn status =
 	fProcesseList.DoContinuous(particle, step, fStack);
       if (status == corsika::process::EProcessReturn::eParticleAbsorbed) {

Framework/Cascade/testCascade.cc

@@ -10,9 +10,15 @@
  */
 
 #include <corsika/cascade/Cascade.h>
+
 #include <corsika/process/ProcessSequence.h>
 #include <corsika/process/stack_inspector/StackInspector.h>
 
+#include <corsika/stack/super_stupid/SuperStupidStack.h>
+
+#include <corsika/geometry/Point.h>
+#include <corsika/geometry/Vector.h>
+
 #include <corsika/setup/SetupStack.h>
 #include <corsika/setup/SetupTrajectory.h>
 using corsika::setup::Trajectory;
@@ -24,6 +30,7 @@ using corsika::setup::Trajectory;
 using namespace corsika;
 using namespace corsika::process;
 using namespace corsika::units;
+using namespace corsika::geometry;
 
 #include <iostream>
 using namespace std;
@@ -35,11 +42,10 @@ public:
   ProcessSplit() {}
 
   template <typename Particle>
-  void MinStepLength(Particle&, Trajectory& ) const {
-  }
+  void MinStepLength(Particle&, setup::Trajectory&) const {}
 
   template <typename Particle, typename Stack>
-  EProcessReturn DoContinuous(Particle&, Trajectory&, Stack&) const {
+  EProcessReturn DoContinuous(Particle&, setup::Trajectory&, Stack&) const {
     return EProcessReturn::eOk;
   }
 
@@ -51,7 +57,11 @@ public:
       fCount++;
     } else {
       p.SetEnergy(E / 2);
-      s.NewParticle().SetEnergy(E / 2);
+      auto pnew = s.NewParticle();
+      // s.Copy(p, pnew);
+      pnew.SetEnergy(E / 2);
+      pnew.SetPosition(p.GetPosition());
+      pnew.SetMomentum(p.GetMomentum());
     }
   }
 
@@ -62,23 +72,28 @@ public:
 private:
 };
 
-template<typename Stack>
+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::Vector<SpeedType::dimension_type> v =
+        p.GetDirection(); 
     corsika::geometry::Line traj(p.GetPosition(), v);
-    return corsika::geometry::Trajectory<corsika::geometry::Line>(traj, 0_s, 100_ns);
+    {
+      CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
+      cout << v.GetComponents(rootCS) << endl;
+    }
+    return corsika::geometry::Trajectory<corsika::geometry::Line>(traj, 100_ns);
   }
 };
 
-
 TEST_CASE("Cascade", "[Cascade]") {
 
   NewtonTracking<setup::Stack> tracking;
-  
+
   stack_inspector::StackInspector<setup::Stack> p0(true);
   ProcessSplit p1;
   const auto sequence = p0 + p1;
@@ -86,10 +101,15 @@ TEST_CASE("Cascade", "[Cascade]") {
 
   corsika::cascade::Cascade EAS(tracking, sequence, stack);
 
+  CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
+
   stack.Clear();
   auto particle = stack.NewParticle();
   EnergyType E0 = 100_GeV;
   particle.SetEnergy(E0);
+  particle.SetPosition(Point(rootCS, {0_m, 0_m, 10_km}));
+  particle.SetMomentum(
+      corsika::stack::super_stupid::MomentumVector(rootCS, {0_GeV, 0_GeV, -1_MeV}));
   EAS.Init();
   EAS.Run();
 

Framework/Geometry/Trajectory.h

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

Framework/Geometry/testGeometry.cc

@@ -152,7 +152,7 @@ TEST_CASE("Trajectories") {
               .norm()
               .magnitude() == Approx(0).margin(absMargin));
 
-    Trajectory<Line> base(line, 0_s, 1_s);
+    Trajectory<Line> base(line, 1_s);
     CHECK(line.GetPosition(2_s).GetCoordinates() ==
           base.GetPosition(2_s).GetCoordinates());
 
@@ -180,7 +180,7 @@ TEST_CASE("Trajectories") {
               .norm()
               .magnitude() == Approx(0).margin(absMargin));
 
-    Trajectory<Helix> const base(helix, 0_s, 1_s);
+    Trajectory<Helix> const base(helix, 1_s);
     CHECK(helix.GetPosition(1234_s).GetCoordinates() ==
           base.GetPosition(1234_s).GetCoordinates());
 

Framework/ProcessSequence/ProcessSequence.h

@@ -19,6 +19,8 @@
 
 #include <corsika/setup/SetupTrajectory.h>
 
+#include <variant>
+
 //#include <type_traits> // still needed ?
 
 using corsika::setup::Trajectory;

Processes/StackInspector/StackInspector.cc

@@ -46,7 +46,7 @@ process::EProcessReturn StackInspector<Stack>::DoContinuous(
     cout << "StackInspector: i=" << setw(5) << fixed << (i++)
 	 << ", id=" << setw(30) << iterP.GetPID()
          << " E=" << setw(15) << scientific << (E / 1_GeV) << " GeV, "
-      //<< " pos=" << pos
+	 << " pos=" << pos
          << endl;
   }
   countStep++;

Setup/SetupTrajectory.h

@@ -16,6 +16,8 @@
 #include <corsika/geometry/Line.h>
 #include <corsika/geometry/Trajectory.h>
 
+#include <corsika/units/PhysicalUnits.h>
+
 #include <variant>
 
 namespace corsika::setup {
@@ -23,11 +25,38 @@ namespace corsika::setup {
   using corsika::geometry::Helix;
   using corsika::geometry::Line;
 
-  typedef std::variant<std::monostate,
-                       corsika::geometry::Trajectory<Line>,
+  /// definition of Trajectory base class, to be used in tracking and cascades
+  typedef std::variant<std::monostate, corsika::geometry::Trajectory<Line>,
                        corsika::geometry::Trajectory<Helix>>
       Trajectory;
 
+  /// helper visitor to modify Particle by moving along Trajectory
+  template <typename Particle>
+  class ParticleUpdate {
+    
+    Particle& fP;
+
+  public:
+    ParticleUpdate(Particle& p)
+        : fP(p) {}
+    void operator()(std::monostate const&) {}
+
+    template <typename T>
+    void operator()(T const& trajectory) {
+      fP.SetPosition(trajectory.GetPosition(1));
+    }
+  };
+
+  /// helper visitor to modify Particle by moving along Trajectory
+  class GetDuration {    
+  public:
+    corsika::units::si::TimeType operator()(std::monostate const&) { return 0*corsika::units::si::second; }
+    template <typename T>
+      corsika::units::si::TimeType operator()(T const &trajectory) {
+      return trajectory.GetDuration();
+    }
+  };
+
 } // namespace corsika::setup
 
 #endif

Stack/SuperStupidStack/SuperStupidStack.h

@@ -69,7 +69,7 @@ namespace corsika::stack {
 #warning this does not really work, nor make sense:
       Vector<SpeedType::dimension_type> GetDirection() const {
 	auto P = GetMomentum();
-	return P/P.norm() * (units::si::meter/units::si::second); }
+	return P/P.norm() * 1e10 * (units::si::meter/units::si::second); }
     
     };