Merge branch 'sibyll' of gitlab.ikp.kit.edu:AirShowerPhysics/corsika into sibyll

CMakeLists.txt

@@ -7,6 +7,9 @@ project (
   LANGUAGES CXX
   )
 
+# as long as there still are modules using it:
+enable_language (Fortran)
+
 # ignore many irrelevant Up-to-date messages during install
 set (CMAKE_INSTALL_MESSAGE LAZY)
 
@@ -21,7 +24,21 @@ set (CMAKE_CXX_STANDARD 17)
 enable_testing ()
 set (CTEST_OUTPUT_ON_FAILURE 1)
 
-ENABLE_LANGUAGE (Fortran)
+# Set a default build type if none was specified
+set(default_build_type "Release")
+if(EXISTS "${CMAKE_SOURCE_DIR}/.git")
+  set(default_build_type "Debug")
+endif()
+ 
+if(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES)
+  message(STATUS "Setting build type to '${default_bluild_type}' as no other was specified.")
+  set(CMAKE_BUILD_TYPE "${default_build_type}" CACHE
+      STRING "Choose the type of build." FORCE)
+  # Set the possible values of build type for cmake-gui
+  set_property(CACHE CMAKE_BUILD_TYPE PROPERTY STRINGS
+    "Debug" "Release" "MinSizeRel" "RelWithDebInfo")
+endif()
+
 
 # unit testing coverage, does not work yet
 #include (CodeCoverage)

Documentation/Examples/CMakeLists.txt

@@ -31,6 +31,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/cascade_example.cc

@@ -12,14 +12,15 @@
 #include <corsika/cascade/Cascade.h>
 #include <corsika/process/ProcessSequence.h>
 #include <corsika/process/stack_inspector/StackInspector.h>
+#include <corsika/process/tracking_line/TrackingLine.h>
 
 #include <corsika/setup/SetupStack.h>
 #include <corsika/setup/SetupTrajectory.h>
 
 #include <corsika/random/RNGManager.h>
 
-#include <corsika/cascade/sibyll2.3c.h>
 #include <corsika/cascade/SibStack.h>
+#include <corsika/cascade/sibyll2.3c.h>
 #include <corsika/process/sibyll/ParticleConversion.h>
 
 #include <corsika/units/PhysicalUnits.h>
@@ -29,9 +30,10 @@ using namespace corsika::process;
 using namespace corsika::units;
 using namespace corsika::particles;
 using namespace corsika::random;
+using namespace corsika::setup;
 
-#include <typeinfo>
 #include <iostream>
+#include <typeinfo>
 using namespace std;
 
 static int fCount = 0;
@@ -41,23 +43,25 @@ public:
   ProcessSplit() {}
 
   template <typename Particle>
-  double MinStepLength(Particle& p) const {
+  double MinStepLength(Particle& p, setup::Trajectory&) const {
 
     const Code corsikaBeamId = p.GetPID();
     
     // beam particles for sibyll : 1, 2, 3 for p, pi, k
     // read from cross section code table
+
     int kBeam   =  process::sibyll::GetSibyllXSCode( corsikaBeamId );   
 
     bool kInteraction = process::sibyll::CanInteract( corsikaBeamId );
     
     /* 
        the target should be defined by the Environment,
-       ideally as full particle object so that the four momenta 
+       ideally as full particle object so that the four momenta
        and the boosts can be defined..
      */
     // target nuclei: A < 18
     // FOR NOW: assume target is oxygen
+
     int kTarget = 16;
     double beamEnergy =  p.GetEnergy() / 1_GeV;
 #warning boost to cm. still missing, sibyll cross section input is cm. energy!
@@ -98,25 +102,27 @@ public:
     
     /*
       what are the units of the output? slant depth or 3space length?
-      
+
     */
-    std::cout << "ProcessSplit: " << "next step (g/cm2): " << next_step << std::endl;
+    std::cout << "ProcessSplit: "
+              << "next step (g/cm2): " << next_step << std::endl;
     return next_step;
   }
 
-  template <typename Particle, typename Trajectory, typename Stack>
-  EProcessReturn DoContinuous(Particle&, Trajectory&, Stack&) const {
+  template <typename Particle, typename Stack>
+  EProcessReturn DoContinuous(Particle&, setup::Trajectory&, Stack&) const {
     // corsika::utls::ignore(p);
     return EProcessReturn::eOk;
   }
 
   template <typename Particle, typename Stack>
   void DoDiscrete(Particle& p, Stack& s) const {
+
     cout << "DoDiscrete: " << p.GetPID() << " interaction? " << process::sibyll::CanInteract( p.GetPID() )  << endl; 
     if( process::sibyll::CanInteract( p.GetPID() ) ){
       cout << "defining coordinates" << endl;
       // coordinate system, get global frame of reference
-      CoordinateSystem rootCS = CoordinateSystem::CreateRootCS();
+      CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
 
       QuantityVector<length_d> const coordinates{0_m, 0_m, 0_m};
       Point pOrig(rootCS, coordinates);
@@ -253,84 +259,87 @@ public:
       //cout << "tot. momentum final (GeV/c): " << Ptot_final.GetComponents() / 1_GeV * si::constants::c << endl;
     }
     }else
+
       p.Delete();
   }
-  
-  
-  void Init() 
-  {
+
+  void Init() {
     fCount = 0;
 
     // define reference frame? --> defines boosts between corsika stack and model stack.
-    
+
     // initialize random numbers for sibyll
     // FOR NOW USE SIBYLL INTERNAL !!!
     //    rnd_ini_();
-    
-    corsika::random::RNGManager & rmng = corsika::random::RNGManager::GetInstance();;
+
+    corsika::random::RNGManager& rmng = corsika::random::RNGManager::GetInstance();
+    ;
     const std::string str_name = "s_rndm";
     rmng.RegisterRandomStream(str_name);
-    
+
     // //    corsika::random::RNG srng;
     // auto srng = rmng.GetRandomStream("s_rndm");
 
     // test random number generator
-    std::cout << "ProcessSplit: " << " test sequence of random numbers."  << std::endl;
+    std::cout << "ProcessSplit: "
+              << " test sequence of random numbers." << std::endl;
     int a = 0;
-    for(int i=0; i<8; ++i)
-      std::cout << i << " " << s_rndm_(a) << std::endl;
-    
-    //initialize Sibyll
+    for (int i = 0; i < 8; ++i) std::cout << i << " " << s_rndm_(a) << std::endl;
+
+    // initialize Sibyll
     sibyll_ini_();
 
     // set particles stable / unstable
     // use stack to loop over particles
     setup::Stack ds;
-    ds.NewParticle().SetPID( Code::Proton  );
-    ds.NewParticle().SetPID( Code::Neutron );
-    ds.NewParticle().SetPID( Code::PiPlus  );
-    ds.NewParticle().SetPID( Code::PiMinus );
-    ds.NewParticle().SetPID( Code::KPlus   );
-    ds.NewParticle().SetPID( Code::KMinus  );
-    ds.NewParticle().SetPID( Code::K0Long  );
-    ds.NewParticle().SetPID( Code::K0Short );
-
-    for( auto &p: ds){
-      int s_id = process::sibyll::ConvertToSibyllRaw( p.GetPID() );
+    ds.NewParticle().SetPID(Code::Proton);
+    ds.NewParticle().SetPID(Code::Neutron);
+    ds.NewParticle().SetPID(Code::PiPlus);
+    ds.NewParticle().SetPID(Code::PiMinus);
+    ds.NewParticle().SetPID(Code::KPlus);
+    ds.NewParticle().SetPID(Code::KMinus);
+    ds.NewParticle().SetPID(Code::K0Long);
+    ds.NewParticle().SetPID(Code::K0Short);
+
+    for (auto& p : ds) {
+      int s_id = process::sibyll::ConvertToSibyllRaw(p.GetPID());
       // set particle stable by setting table value negative
-      cout << "ProcessSplit: Init: setting " << p.GetPID() << "(" << s_id << ")" << " stable in Sibyll .." << endl;
-      s_csydec_.idb[ s_id ] = -s_csydec_.idb[ s_id-1 ];
+      cout << "ProcessSplit: Init: setting " << p.GetPID() << "(" << s_id << ")"
+           << " stable in Sibyll .." << endl;
+      s_csydec_.idb[s_id] = -s_csydec_.idb[s_id - 1];
       p.Delete();
     }
   }
-  
+
   int GetCount() { return fCount; }
 
 private:
 };
 
-double s_rndm_(int &)
-{
-  static corsika::random::RNG& rmng = corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");;
-  return rmng()/(double)rmng.max();
+double s_rndm_(int&) {
+  static corsika::random::RNG& rmng =
+      corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");
+  ;
+  return rmng() / (double)rmng.max();
 }
 
-
-int main(){
+int main() {
 
   // coordinate system, get global frame of reference
-  CoordinateSystem rootCS = CoordinateSystem::CreateRootCS();
+  CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
   
   QuantityVector<length_d> const coordinates{0_m, 0_m, 0_m};
   Point pOrig(rootCS, coordinates);    
   
-  stack_inspector::StackInspector<setup::Stack, setup::Trajectory> p0(true);
+  //  stack_inspector::StackInspector<setup::Stack, setup::Trajectory> p0(true);
+  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(sequence, stack);
+  corsika::cascade::Cascade EAS(tracking, sequence, stack);
 
   stack.Clear();
   auto particle = stack.NewParticle();
@@ -341,10 +350,10 @@ int main(){
 					   0. * 1_GeV / si::constants::c,
 					   P0);
   particle.SetEnergy(E0);
+
   particle.SetMomentum(plab);
   particle.SetPID( Code::Proton );
   EAS.Init();
   EAS.Run();
   cout << "Result: E0=" << E0 / 1_GeV << "GeV, count=" << p1.GetCount() << endl;
-  
 }

Documentation/Examples/geometry_example.cc

@@ -9,8 +9,8 @@
  * the license.
  */
 
-#include <corsika/geometry/CoordinateSystem.h>
 #include <corsika/geometry/Point.h>
+#include <corsika/geometry/RootCoordinateSystem.h>
 #include <corsika/geometry/Sphere.h>
 #include <corsika/geometry/Vector.h>
 #include <corsika/units/PhysicalUnits.h>
@@ -19,12 +19,14 @@
 #include <iostream>
 #include <typeinfo>
 
+using namespace corsika;
 using namespace corsika::geometry;
 using namespace corsika::units::si;
 
 int main() {
   // define the root coordinate system
-  auto const root = CoordinateSystem::CreateRootCS();
+  geometry::CoordinateSystem& root =
+      geometry::RootCoordinateSystem::GetInstance().GetRootCS();
 
   // another CS defined by a translation relative to the root CS
   CoordinateSystem cs2 = root.translate({0_m, 0_m, 1_m});

Documentation/Examples/helix_example.cc

@@ -9,20 +9,22 @@
  * the license.
  */
 
-#include <corsika/geometry/CoordinateSystem.h>
 #include <corsika/geometry/Helix.h>
 #include <corsika/geometry/Point.h>
+#include <corsika/geometry/RootCoordinateSystem.h>
 #include <corsika/geometry/Vector.h>
 #include <corsika/units/PhysicalUnits.h>
 #include <array>
 #include <cstdlib>
 #include <iostream>
 
+using namespace corsika;
 using namespace corsika::geometry;
 using namespace corsika::units::si;
 
 int main() {
-  auto const root = CoordinateSystem::CreateRootCS();
+  geometry::CoordinateSystem& root =
+      geometry::RootCoordinateSystem::GetInstance().GetRootCS();
 
   Point const r0(root, {0_m, 0_m, 0_m});
   auto const omegaC = 2 * M_PI * 1_Hz;

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,8 +88,8 @@ void modular() {
   const auto sequence = m1 + m2 + m3 + m4;
 
   DummyData p;
-  DummyTrajectory t;
   DummyStack s;
+  Trajectory t;
 
   const int n = 100000000;
   for (int i = 0; i < n; ++i) { sequence.DoContinuous(p, t, s); }

Framework/Cascade/CMakeLists.txt

@@ -15,6 +15,7 @@ set (
 
 #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,25 @@ 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::setup::Trajectory step = fTracking.GetTrack(particle);
+      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, 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/SibStack.h

 #ifndef _include_sibstack_h_
 #define _include_sibstack_h_
 
-#include <vector>
 #include <string>
+#include <vector>
 
-#include <corsika/stack/Stack.h>
 #include <corsika/cascade/sibyll2.3c.h>
 #include <corsika/process/sibyll/ParticleConversion.h>
+#include <corsika/stack/Stack.h>
 
 using namespace std;
 using namespace corsika::stack;
@@ -14,10 +14,10 @@ using namespace corsika::units;
 using namespace corsika::geometry;
 
 class SibStackData {
-  
- public:
+
+public:
   void Init();
-  
+
   void Clear() { s_plist_.np = 0; }
   
   int GetSize() const { return s_plist_.np;  }
@@ -40,7 +40,7 @@ class SibStackData {
   
   super_stupid::MomentumVector GetMomentum(const int i) const
   {
-    CoordinateSystem rootCS = CoordinateSystem::CreateRootCS();
+    CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
     corsika::geometry::QuantityVector<momentum_d> components{ s_plist_.p[0][i] * 1_GeV / si::constants::c , s_plist_.p[1][i] * 1_GeV / si::constants::c, s_plist_.p[2][i] * 1_GeV / si::constants::c};
     super_stupid::MomentumVector v1(rootCS,components);
     return v1;
@@ -50,18 +50,20 @@ class SibStackData {
     s_plist_.llist[i1] = s_plist_.llist[i2];
     s_plist_.p[3][i1] = s_plist_.p[3][i2];
   }
-  
- protected:
+
+protected:
   void IncrementSize() { s_plist_.np++; }
-  void DecrementSize() { if ( s_plist_.np>0) { s_plist_.np--; } }
+  void DecrementSize() {
+    if (s_plist_.np > 0) { s_plist_.np--; }
+  }
 };
 
-
-template<typename StackIteratorInterface>
+template <typename StackIteratorInterface>
 class ParticleInterface : public ParticleBase<StackIteratorInterface> {
   using ParticleBase<StackIteratorInterface>::GetStackData;
   using ParticleBase<StackIteratorInterface>::GetIndex;
- public:
+
+public:
   void SetEnergy(const double v) { GetStackData().SetEnergy(GetIndex(), v); }
   EnergyType GetEnergy() const { return GetStackData().GetEnergy(GetIndex()); }
   void SetPID(const int v) { GetStackData().SetId(GetIndex(), v); }
@@ -70,7 +72,6 @@ class ParticleInterface : public ParticleBase<StackIteratorInterface> {
   
 };
 
-
 typedef Stack<SibStackData, ParticleInterface> SibStack;
 
 #endif

Framework/Cascade/sibyll2.3c.h

@@ -3,94 +3,88 @@
 //----------------------------------------------
 //  C++ interface for the SIBYLL event generator
 //----------------------------------------------
-//wrapper
+// wrapper
 
+extern "C" {
 
+typedef char s_name[6];
 
+// SIBYLL particle stack (FORTRAN COMMON)
+// variables are: np : numer of particles on stack
+//                 p : 4momentum + mass of particles on stack
+//             llist : id of particles on stack
+extern struct {
+  double p[5][8000];
+  int llist[8000];
+  int np;
+} s_plist_;
 
-extern"C"{
+extern struct {
+  double cbr[223 + 16 + 12 + 8];
+  int kdec[1338 + 6 * (16 + 12 + 8)];
+  int lbarp[99];
+  int idb[99];
+} s_csydec_;
 
-  typedef char s_name[6];
+// additional particle stack for the mother particles of unstable particles
+// stable particles have entry zero
+extern struct { int llist1[8000]; } s_plist1_;
 
-  // SIBYLL particle stack (FORTRAN COMMON)
-  // variables are: np : numer of particles on stack
-  //                 p : 4momentum + mass of particles on stack
-  //             llist : id of particles on stack
-  extern struct { 
-    double p[5][8000]; 
-    int llist[8000]; 
-    int np; 
-  }s_plist_;
+// tables with particle properties
+// charge, strangeness and baryon number
+extern struct {
+  int ichp[99];
+  int istr[99];
+  int ibar[99];
+} s_chp_;
 
+// tables with particle properties
+// mass and mass squared
+extern struct {
+  double am[99];
+  double am2[99];
+} s_mass1_;
 
-  extern struct {
-    double  cbr[223+16+12+8];
-    int    kdec[1338+6*(16+12+8)];
-    int    lbarp[99];
-    int     idb[99];
-  }s_csydec_;
+// table with particle names
+extern struct { char namp[6][99]; } s_cnam_;
 
-  
-  // additional particle stack for the mother particles of unstable particles
-  // stable particles have entry zero
-  extern struct { 
-    int llist1[8000]; 
-  }s_plist1_;
+// debug info
+extern struct {
+  int ncall;
+  int ndebug;
+  int lun;
+} s_debug_;
 
+// lund random generator setup
+// extern struct {int mrlu[6]; float rrlu[100]; }ludatr_;
 
-  // tables with particle properties
-  // charge, strangeness and baryon number
-  extern struct {
-    int ichp[99]; 
-    int istr[99]; 
-    int ibar[99]; 
-  }s_chp_;
+// sibyll main subroutine
+void sibyll_(int&, int&, double&);
 
-  // tables with particle properties
-  // mass and mass squared
-  extern struct {
-    double am[99]; 
-    double am2[99]; 
-  }s_mass1_;
+// subroutine to initiate sibyll
+void sibyll_ini_();
 
-  // table with particle names
-  extern struct {char namp[6][99];}s_cnam_;
+// subroutine to SET DECAYS
+void dec_ini_();
 
-  // debug info
-  extern struct {int ncall; int ndebug; int lun; }s_debug_;
+// subroutine to initiate random number generator
+// void rnd_ini_();
 
-  // lund random generator setup
-  //extern struct {int mrlu[6]; float rrlu[100]; }ludatr_;
+// print event
+void sib_list_(int&);
 
+// decay routine
+void decsib_();
 
-  // sibyll main subroutine
-  void sibyll_(int&,int&,double&);
-    
-  // subroutine to initiate sibyll
-  void sibyll_ini_();
+// interaction length
+// double fpni_(double&, int&);
 
-  // subroutine to SET DECAYS
-  void dec_ini_();
+void sib_sigma_hnuc_(int&, int&, double&, double&, double&);
+void sib_sigma_hp_(int&, double&, double&, double&, double&, double*, double&, double&);
 
-  // subroutine to initiate random number generator
-  //void rnd_ini_();
-  
-  // print event
-  void sib_list_(int&);
+double s_rndm_(int&);
 
-  // decay routine
-  void decsib_();
-
-  // interaction length
-  //double fpni_(double&, int&);
-
-  void sib_sigma_hnuc_(int&,int&,double&,double&,double&);
-  void sib_sigma_hp_(int&,double&,double&,double&,double&,double*,double&,double&);
-  
-  double s_rndm_(int&);
-
-  // phojet random generator setup
-  void pho_rndin_(int&,int&,int&,int&);
+// phojet random generator setup
+void pho_rndin_(int&, int&, int&, int&);
 }
 #endif
-

Framework/Cascade/testCascade.cc

@@ -10,11 +10,19 @@
  */
 
 #include <corsika/cascade/Cascade.h>
+
 #include <corsika/process/ProcessSequence.h>
 #include <corsika/process/stack_inspector/StackInspector.h>
+#include <corsika/process/tracking_line/TrackingLine.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;
 
 #define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one
                           // cpp file
@@ -23,6 +31,7 @@
 using namespace corsika;
 using namespace corsika::process;
 using namespace corsika::units;
+using namespace corsika::geometry;
 
 #include <iostream>
 using namespace std;
@@ -34,13 +43,10 @@ public:
   ProcessSplit() {}
 
   template <typename Particle>
-  double MinStepLength(Particle&) const {
-    return 0;
-  }
+  void MinStepLength(Particle&, setup::Trajectory&) const {}
 
-  template <typename Particle, typename Trajectory, typename Stack>
-  EProcessReturn DoContinuous(Particle&, Trajectory&, Stack&) const {
-    // corsika::utls::ignore(p);
+  template <typename Particle, typename Stack>
+  EProcessReturn DoContinuous(Particle&, setup::Trajectory&, Stack&) const {
     return EProcessReturn::eOk;
   }
 
@@ -52,30 +58,42 @@ 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());
     }
   }
 
   void Init() { fCount = 0; }
 
-  int GetCount() { return fCount; }
+  int GetCount() const { return fCount; }
 
 private:
 };
 
+
 TEST_CASE("Cascade", "[Cascade]") {
 
-  stack_inspector::StackInspector<setup::Stack, setup::Trajectory> p0(true);
+  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(sequence, stack);
+  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*newton*second, 0*newton*second, -1*newton*second}));
   EAS.Init();
   EAS.Run();
 

Framework/Geometry/CMakeLists.txt

@@ -11,11 +11,12 @@ set (
   Line.h
   Sphere.h
   CoordinateSystem.h
+  RootCoordinateSystem.h
   Helix.h
   BaseVector.h
   QuantityVector.h
   Trajectory.h
-  BaseTrajectory.h
+  # BaseTrajectory.h
   )
 
 set (

Framework/Geometry/CoordinateSystem.h

@@ -21,6 +21,8 @@ typedef Eigen::Translation<double, 3> EigenTranslation;
 
 namespace corsika::geometry {
 
+  class RootCoordinateSystem;
+
   using corsika::units::si::length_d;
 
   class CoordinateSystem {
@@ -30,14 +32,17 @@ namespace corsika::geometry {
     CoordinateSystem(CoordinateSystem const& reference, EigenTransform const& transf)
         : reference(&reference)
         , transf(transf) {}
-        
+
     CoordinateSystem()
         : // for creating the root CS
         transf(EigenTransform::Identity()) {}
 
+  protected:
+    static auto CreateCS() { return CoordinateSystem(); }
+    friend corsika::geometry::RootCoordinateSystem; /// this is the only class that can
+                                                    /// creat ONE unique root CS
+
   public:
-    static auto CreateRootCS() { return CoordinateSystem(); }
-  
     static EigenTransform GetTransformation(CoordinateSystem const& c1,
                                             CoordinateSystem const& c2);
 
@@ -55,9 +60,9 @@ namespace corsika::geometry {
 
     auto rotate(QuantityVector<phys::units::length_d> axis, double angle) const {
       if (axis.eVector.isZero()) {
-          throw std::string("null-vector given as axis parameter");
+        throw std::string("null-vector given as axis parameter");
       }
-      
+
       EigenTransform const rotation{Eigen::AngleAxisd(angle, axis.eVector.normalized())};
 
       return CoordinateSystem(*this, rotation);
@@ -66,9 +71,9 @@ namespace corsika::geometry {
     auto translateAndRotate(QuantityVector<phys::units::length_d> translation,
                             QuantityVector<phys::units::length_d> axis, double angle) {
       if (axis.eVector.isZero()) {
-          throw std::string("null-vector given as axis parameter");
+        throw std::string("null-vector given as axis parameter");
       }
-                                
+
       EigenTransform const transf{Eigen::AngleAxisd(angle, axis.eVector.normalized()) *
                                   EigenTranslation(translation.eVector)};
 

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

@@ -19,7 +19,7 @@
 namespace corsika::geometry {
 
   class Line {
-    
+
     using VelocityVec = Vector<corsika::units::si::SpeedType::dimension_type>;
 
     Point const r0;
@@ -30,13 +30,11 @@ namespace corsika::geometry {
         : r0(pR0)
         , v0(pV0) {}
 
-    Point GetPosition(corsika::units::si::TimeType t) const {
-      return r0 + v0 * t;
-    }
-    
-    LengthType DistanceBetween(corsika::units::si::TimeType t1,
-                               corsika::units::si::TimeType t2) const {
-      assert(t2 >= t1);
+    Point GetPosition(corsika::units::si::TimeType t) const { return r0 + v0 * t; }
+
+    LengthType GetDistanceBetween(corsika::units::si::TimeType t1,
+                                  corsika::units::si::TimeType t2) const {
+      // assert(t2 >= t1);
       return v0.norm() * (t2 - t1);
     }
   };

Framework/Geometry/Point.h

@@ -34,8 +34,11 @@ namespace corsika::geometry {
     Point(CoordinateSystem const& cs, LengthType x, LengthType y, LengthType z)
         : BaseVector<phys::units::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; }
 
+    /// 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;

Framework/Geometry/RootCoordinateSystem.h

+#ifndef _include_corsika_geometry_rootcoordinatesystem_h_
+#define _include_corsika_geometry_rootcoordinatesystem_h_
+
+#include <corsika/utl/Singleton.h>
+
+#include <corsika/geometry/CoordinateSystem.h>
+
+/*!
+ * This is the only way to get a root-coordinate system, and it is a
+ * singleton. All other CoordinateSystems must be relative to the
+ * RootCoordinateSystem
+ */
+
+namespace corsika::geometry {
+
+  class RootCoordinateSystem : public corsika::utl::Singleton<RootCoordinateSystem> {
+
+    friend class corsika::utl::Singleton<RootCoordinateSystem>;
+
+  protected:
+    RootCoordinateSystem() {}
+
+  public:
+    corsika::geometry::CoordinateSystem& GetRootCS() { return fRootCS; }
+    const corsika::geometry::CoordinateSystem& GetRootCS() const { return fRootCS; }
+
+  private:
+    corsika::geometry::CoordinateSystem fRootCS; // THIS IS IT
+  };
+
+} // namespace corsika::geometry
+
+#endif

Framework/Geometry/Trajectory.h

@@ -12,37 +12,39 @@
 #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 BaseTrajectory {
+  class Trajectory : public T {
 
-    T fTraj;
+    corsika::units::si::TimeType fTimeLength;
 
   public:
-    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) {}
+    using T::GetDistanceBetween;
+    using T::GetPosition;
+
+    Trajectory(T const& theT, corsika::units::si::TimeType timeLength)
+        : T(theT)
+        , fTimeLength(timeLength) {}
 
-    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));
-    }
-    
-    LengthType GetDistance(corsika::units::si::TimeType t1,
-			   corsika::units::si::TimeType t2) const {
-      return fTraj.DistanceBetween(t1, t2);
-    }
-     
+    Point GetPosition(const double u) const { return T::GetPosition(fTimeLength * u); }
+
+    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

@@ -17,6 +17,7 @@
 #include <corsika/geometry/Helix.h>
 #include <corsika/geometry/Line.h>
 #include <corsika/geometry/Point.h>
+#include <corsika/geometry/RootCoordinateSystem.h>
 #include <corsika/geometry/Sphere.h>
 #include <corsika/geometry/Trajectory.h>
 #include <corsika/units/PhysicalUnits.h>
@@ -28,7 +29,7 @@ using namespace corsika::units::si;
 double constexpr absMargin = 1.0e-8;
 
 TEST_CASE("transformations between CoordinateSystems") {
-  CoordinateSystem rootCS = CoordinateSystem::CreateRootCS();
+  CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
 
   REQUIRE(CoordinateSystem::GetTransformation(rootCS, rootCS)
               .isApprox(EigenTransform::Identity()));
@@ -44,10 +45,11 @@ TEST_CASE("transformations between CoordinateSystems") {
   REQUIRE((p1.GetCoordinates(rootCS) - coordinates).norm().magnitude() ==
           Approx(0).margin(absMargin));
 
+  /*
   SECTION("unconnected CoordinateSystems") {
     CoordinateSystem rootCS2 = CoordinateSystem::CreateRootCS();
     REQUIRE_THROWS(CoordinateSystem::GetTransformation(rootCS, rootCS2));
-  }
+    }*/
 
   SECTION("translations") {
     QuantityVector<length_d> const translationVector{0_m, 4_m, 0_m};
@@ -126,7 +128,7 @@ TEST_CASE("transformations between CoordinateSystems") {
 }
 
 TEST_CASE("Sphere") {
-  CoordinateSystem rootCS = CoordinateSystem::CreateRootCS();
+  CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
   Point center(rootCS, {0_m, 3_m, 4_m});
   Sphere sphere(center, 5_m);
 
@@ -137,7 +139,7 @@ TEST_CASE("Sphere") {
 }
 
 TEST_CASE("Trajectories") {
-  CoordinateSystem rootCS = CoordinateSystem::CreateRootCS();
+  CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
   Point r0(rootCS, {0_m, 0_m, 0_m});
 
   SECTION("Line") {
@@ -145,16 +147,16 @@ TEST_CASE("Trajectories") {
                                          {1_m / second, 0_m / second, 0_m / second});
 
     Line const line(r0, v0);
-    CHECK((line.GetPosition(2_s).GetCoordinates() -
-           QuantityVector<length_d>(2_m, 0_m, 0_m))
-              .norm()
-              .magnitude() == Approx(0).margin(absMargin));
+    CHECK(
+        (line.GetPosition(2_s).GetCoordinates() - QuantityVector<length_d>(2_m, 0_m, 0_m))
+            .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());
 
-    CHECK(base.GetDistance(1_s, 2_s) / 1_m == Approx(1));
+    CHECK(base.GetDistanceBetween(1_s, 2_s) / 1_m == Approx(1));
   }
 
   SECTION("Helix") {
@@ -178,10 +180,10 @@ 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());
 
-    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,8 @@ add_executable (
 
 target_link_libraries (
   testProcessSequence
+  CORSIKAsetup
+  CORSIKAgeometry
   CORSIKAprocesssequence
   CORSIKAthirdparty # for catch2
   )