added new global root-coordinate-system machinery

CMakeLists.txt

@@ -21,6 +21,22 @@ set (CMAKE_CXX_STANDARD 17)
 enable_testing ()
 set (CTEST_OUTPUT_ON_FAILURE 1)
 
+# 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)
 ##set(COVERAGE_LCOV_EXCLUDES 'Documentation/*')

Documentation/Examples/geometry_example.cc

@@ -9,7 +9,7 @@
  * the license.
  */
 
-#include <corsika/geometry/CoordinateSystem.h>
+#include <corsika/geometry/RootCoordinateSystem.h>
 #include <corsika/geometry/Point.h>
 #include <corsika/geometry/Sphere.h>
 #include <corsika/geometry/Vector.h>
@@ -19,12 +19,13 @@
 #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,7 +9,7 @@
  * the license.
  */
 
-#include <corsika/geometry/CoordinateSystem.h>
+#include <corsika/geometry/RootCoordinateSystem.h>
 #include <corsika/geometry/Helix.h>
 #include <corsika/geometry/Point.h>
 #include <corsika/geometry/Vector.h>
@@ -18,11 +18,12 @@
 #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

@@ -89,12 +89,7 @@ void modular() {
 
   DummyData p;
   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));
+  Trajectory t;
 
   const int n = 100000000;
   for (int i = 0; i < n; ++i) { sequence.DoContinuous(p, t, s); }

Framework/Cascade/testCascade.cc

@@ -36,12 +36,10 @@ public:
 
   template <typename Particle>
   void MinStepLength(Particle&, Trajectory& ) const {
-    //return 0;
   }
 
   template <typename Particle, typename Stack>
   EProcessReturn DoContinuous(Particle&, Trajectory&, Stack&) const {
-    // corsika::utls::ignore(p);
     return EProcessReturn::eOk;
   }
 

Framework/Geometry/BaseVector.h

@@ -31,6 +31,8 @@ namespace corsika::geometry {
     BaseVector(CoordinateSystem const& pCS, QuantityVector<dim> pQVector)
         : qVector(pQVector)
         , cs(&pCS) {}
+
+
   };
 
 } // namespace corsika::geometry

Framework/Geometry/CMakeLists.txt

@@ -11,6 +11,7 @@ set (
   Line.h
   Sphere.h
   CoordinateSystem.h
+  RootCoordinateSystem.h
   Helix.h
   BaseVector.h
   QuantityVector.h

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 {
@@ -35,8 +37,11 @@ namespace corsika::geometry {
         : // 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);

Framework/Geometry/Point.h

@@ -34,8 +34,10 @@ 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
+    
+  };
+  
+}
+
+#endif

Framework/Geometry/testGeometry.cc

@@ -14,6 +14,7 @@
 #include <catch2/catch.hpp>
 
 #include <corsika/geometry/CoordinateSystem.h>
+#include <corsika/geometry/RootCoordinateSystem.h>
 #include <corsika/geometry/Helix.h>
 #include <corsika/geometry/Line.h>
 #include <corsika/geometry/Point.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") {

Framework/ProcessSequence/CMakeLists.txt

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

Processes/StackInspector/StackInspector.cc

@@ -11,6 +11,7 @@
 
 #include <corsika/process/stack_inspector/StackInspector.h>
 #include <corsika/units/PhysicalUnits.h>
+#include <corsika/geometry/RootCoordinateSystem.h>
 
 #include <corsika/logging/Logger.h>
 
@@ -32,7 +33,7 @@ StackInspector<Stack>::~StackInspector() {}
 
 template <typename Stack>
 process::EProcessReturn StackInspector<Stack>::DoContinuous(
-    Particle&, corsika::setup::Trajectory&, Stack& s) const {
+    Particle&, setup::Trajectory&, Stack& s) const {
   static int countStep = 0;
   if (!fReport) return EProcessReturn::eOk;
   [[maybe_unused]] int i = 0;
@@ -40,20 +41,22 @@ process::EProcessReturn StackInspector<Stack>::DoContinuous(
   for (auto& iterP : s) {
     EnergyType E = iterP.GetEnergy();
     Etot += E;
-    cout << "i=" << setw(5) << fixed << (i++)
+    geometry::CoordinateSystem& rootCS = geometry::RootCoordinateSystem::GetInstance().GetRootCS(); // for printout
+    auto pos = iterP.GetPosition().GetCoordinates(rootCS);
+    cout << "StackInspector: i=" << setw(5) << fixed << (i++)
 	 << ", id=" << setw(30) << iterP.GetPID()
          << " E=" << setw(15) << scientific << (E / 1_GeV) << " GeV, "
-         //<< " pos=" << iterP.GetPosition()
+      //<< " pos=" << pos
          << endl;
   }
   countStep++;
-  cout << countStep << " " << s.GetSize() << " " << Etot / 1_GeV << " " << endl;
+  cout << "StackInspector: nStep=" << countStep << " stackSize=" << s.GetSize() << " Estack=" << Etot / 1_GeV << " GeV" << endl;
   return EProcessReturn::eOk;
 }
 
 template <typename Stack>
 void StackInspector<Stack>::MinStepLength(Particle&,
-                                          corsika::setup::Trajectory&) const {
+                                          setup::Trajectory&) const {
   // return 0;
 }
 
@@ -62,4 +65,4 @@ void StackInspector<Stack>::Init() {}
 
 #include <corsika/setup/SetupStack.h>
 
-template class corsika::process::stack_inspector::StackInspector<setup::Stack>;
+template class process::stack_inspector::StackInspector<setup::Stack>;

Stack/SuperStupidStack/SuperStupidStack.h

@@ -16,27 +16,29 @@
 #include <corsika/stack/Stack.h>
 #include <corsika/units/PhysicalUnits.h>
 
-#include <corsika/geometry/CoordinateSystem.h> // remove 
+#include <corsika/geometry/RootCoordinateSystem.h> // remove 
 #include <corsika/geometry/Point.h>
 #include <corsika/geometry/Vector.h>
 
 #include <vector>
 #include <algorithm>
 
+using namespace corsika;
+
 namespace corsika::stack {
 
   namespace super_stupid {
 
-    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;
-    using corsika::units::si::energy_d;
-    using corsika::geometry::Point;
-    using corsika::geometry::Vector;
+    using particles::Code;
+    using units::si::EnergyType;
+    using units::si::TimeType;
+    using units::si::SpeedType;
+    using units::si::second;
+    using units::si::meter;
+    using units::si::joule;
+    using units::si::energy_d;
+    using geometry::Point;
+    using geometry::Vector;
 
 #warning replace this with a proper momentum vector:
     typedef Vector<energy_d> MomentumVector; // should be momentum_d !!!
@@ -67,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() * (corsika::units::si::meter/corsika::units::si::second); }
+	return P/P.norm() * (units::si::meter/units::si::second); }
     
     };
 
@@ -128,7 +130,7 @@ namespace corsika::stack {
         fDataPID.push_back(Code::Unknown);
         fDataE.push_back(0 * joule);
 #warning this here makes no sense: see issue #48
-	auto const dummyCS = corsika::geometry::CoordinateSystem::CreateRootCS();
+	geometry::CoordinateSystem& dummyCS = geometry::RootCoordinateSystem::GetInstance().GetRootCS(); 
 	fMomentum.push_back(MomentumVector(dummyCS,
 					   {0 * joule, 0 * joule, 0 * joule}));	
 	fPosition.push_back(Point(dummyCS,
@@ -150,7 +152,7 @@ namespace corsika::stack {
 
       std::vector<Code> fDataPID;
       std::vector<EnergyType> fDataE;
-      std::vector<Vector<corsika::units::si::energy_d>> fMomentum; // should be Momentum !!!!
+      std::vector<Vector<units::si::energy_d>> fMomentum; // should be Momentum !!!!
       std::vector<Point> fPosition;
       std::vector<TimeType> fTime;
 

Stack/SuperStupidStack/testSuperStupidStack.cc

@@ -11,6 +11,7 @@
 
 #include <corsika/stack/super_stupid/SuperStupidStack.h>
 #include <corsika/units/PhysicalUnits.h>
+#include <corsika/geometry/RootCoordinateSystem.h> 
 
 using namespace corsika::geometry;
 using namespace corsika::units::si;
@@ -33,9 +34,9 @@ TEST_CASE("SuperStupidStack", "[stack]") {
 
     SuperStupidStack s;
     auto p = s.NewParticle();
-    p.SetPID(corsika::particles::Code::Electron);
+    p.SetPID(particles::Code::Electron);
     p.SetEnergy(1.5_GeV);
-    auto const dummyCS = corsika::geometry::CoordinateSystem::CreateRootCS();
+    geometry::CoordinateSystem& dummyCS = geometry::RootCoordinateSystem::GetInstance().GetRootCS(); 
     p.SetMomentum(MomentumVector(dummyCS, {1 * joule, 1 * joule, 1 * joule}));	
     p.SetPosition(Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}));
     p.SetTime(100_s);
@@ -43,7 +44,7 @@ TEST_CASE("SuperStupidStack", "[stack]") {
     // read
     REQUIRE(s.GetSize() == 1);
     auto pout = s.GetNextParticle();
-    REQUIRE(pout.GetPID() == corsika::particles::Code::Electron);
+    REQUIRE(pout.GetPID() == particles::Code::Electron);
     REQUIRE(pout.GetEnergy() == 1.5_GeV);
 #warning Fix the next two lines:
     //REQUIRE(pout.GetMomentum() == MomentumVector(dummyCS, {1 * joule, 1 * joule, 1 * joule}));