From a0f0fcfeb7c4590e7bbfd064dbf5a0e7c5c6c063 Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Sun, 2 Dec 2018 12:24:21 +0100
Subject: [PATCH] added new global root-coordinate-system machinery
---
CMakeLists.txt | 16 +++++++++
Documentation/Examples/geometry_example.cc | 5 +--
Documentation/Examples/helix_example.cc | 5 +--
.../Examples/staticsequence_example.cc | 7 +---
Framework/Cascade/testCascade.cc | 2 --
Framework/Geometry/BaseVector.h | 2 ++
Framework/Geometry/CMakeLists.txt | 1 +
Framework/Geometry/CoordinateSystem.h | 7 +++-
Framework/Geometry/Point.h | 2 ++
Framework/Geometry/RootCoordinateSystem.h | 36 +++++++++++++++++++
Framework/Geometry/testGeometry.cc | 10 +++---
Framework/ProcessSequence/CMakeLists.txt | 1 +
Processes/StackInspector/StackInspector.cc | 15 ++++----
Stack/SuperStupidStack/SuperStupidStack.h | 30 ++++++++--------
.../SuperStupidStack/testSuperStupidStack.cc | 7 ++--
15 files changed, 106 insertions(+), 40 deletions(-)
create mode 100644 Framework/Geometry/RootCoordinateSystem.h
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 67149c9f..3e8c642c 100644
--- a/CMakeLists.txt
+++ b/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/*')
diff --git a/Documentation/Examples/geometry_example.cc b/Documentation/Examples/geometry_example.cc
index 048347ae..ebf824c0 100644
--- a/Documentation/Examples/geometry_example.cc
+++ b/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});
diff --git a/Documentation/Examples/helix_example.cc b/Documentation/Examples/helix_example.cc
index 2797a30a..2c223257 100644
--- a/Documentation/Examples/helix_example.cc
+++ b/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;
diff --git a/Documentation/Examples/staticsequence_example.cc b/Documentation/Examples/staticsequence_example.cc
index 9a29d6f7..07cc7561 100644
--- a/Documentation/Examples/staticsequence_example.cc
+++ b/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); }
diff --git a/Framework/Cascade/testCascade.cc b/Framework/Cascade/testCascade.cc
index 3ade3948..6d6258d2 100644
--- a/Framework/Cascade/testCascade.cc
+++ b/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;
}
diff --git a/Framework/Geometry/BaseVector.h b/Framework/Geometry/BaseVector.h
index 2cd1c3e4..50cf1b07 100644
--- a/Framework/Geometry/BaseVector.h
+++ b/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
diff --git a/Framework/Geometry/CMakeLists.txt b/Framework/Geometry/CMakeLists.txt
index 78f0dee4..2b0c8d3a 100644
--- a/Framework/Geometry/CMakeLists.txt
+++ b/Framework/Geometry/CMakeLists.txt
@@ -11,6 +11,7 @@ set (
Line.h
Sphere.h
CoordinateSystem.h
+ RootCoordinateSystem.h
Helix.h
BaseVector.h
QuantityVector.h
diff --git a/Framework/Geometry/CoordinateSystem.h b/Framework/Geometry/CoordinateSystem.h
index 089a3f80..a75ff439 100644
--- a/Framework/Geometry/CoordinateSystem.h
+++ b/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);
diff --git a/Framework/Geometry/Point.h b/Framework/Geometry/Point.h
index 26de8510..9570c880 100644
--- a/Framework/Geometry/Point.h
+++ b/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;
diff --git a/Framework/Geometry/RootCoordinateSystem.h b/Framework/Geometry/RootCoordinateSystem.h
new file mode 100644
index 00000000..58cc20c6
--- /dev/null
+++ b/Framework/Geometry/RootCoordinateSystem.h
@@ -0,0 +1,36 @@
+#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
diff --git a/Framework/Geometry/testGeometry.cc b/Framework/Geometry/testGeometry.cc
index 274ec7d5..874c2ce0 100644
--- a/Framework/Geometry/testGeometry.cc
+++ b/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") {
diff --git a/Framework/ProcessSequence/CMakeLists.txt b/Framework/ProcessSequence/CMakeLists.txt
index 9e95f0c8..c1f5a458 100644
--- a/Framework/ProcessSequence/CMakeLists.txt
+++ b/Framework/ProcessSequence/CMakeLists.txt
@@ -42,6 +42,7 @@ add_executable (
target_link_libraries (
testProcessSequence
+ CORSIKAsetup
CORSIKAgeometry
CORSIKAprocesssequence
CORSIKAthirdparty # for catch2
diff --git a/Processes/StackInspector/StackInspector.cc b/Processes/StackInspector/StackInspector.cc
index 7480d12d..94abaff7 100644
--- a/Processes/StackInspector/StackInspector.cc
+++ b/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>;
diff --git a/Stack/SuperStupidStack/SuperStupidStack.h b/Stack/SuperStupidStack/SuperStupidStack.h
index fa6bf4e2..27e97f4a 100644
--- a/Stack/SuperStupidStack/SuperStupidStack.h
+++ b/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;
diff --git a/Stack/SuperStupidStack/testSuperStupidStack.cc b/Stack/SuperStupidStack/testSuperStupidStack.cc
index 93f0a095..112b6844 100644
--- a/Stack/SuperStupidStack/testSuperStupidStack.cc
+++ b/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}));
--
GitLab