From 7f8b3f7ddf040a341b267eec675f37353bddea8e Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Sat, 1 Dec 2018 23:21:47 +0100
Subject: [PATCH] working cascade with tracking
---
Documentation/Examples/CMakeLists.txt | 1 +
.../Examples/staticsequence_example.cc | 13 +++-
Framework/Cascade/CMakeLists.txt | 2 +
Framework/Cascade/Cascade.h | 58 ++++++++--------
Framework/Cascade/testCascade.cc | 28 ++++++--
Framework/Geometry/CMakeLists.txt | 2 +-
Framework/Geometry/Helix.h | 6 +-
Framework/Geometry/Line.h | 6 +-
Framework/Geometry/Trajectory.h | 32 +++++----
Framework/Geometry/testGeometry.cc | 4 +-
Framework/ProcessSequence/CMakeLists.txt | 1 +
Framework/ProcessSequence/ProcessSequence.h | 13 ++--
.../ProcessSequence/testProcessSequence.cc | 68 +++++++++++--------
Processes/StackInspector/StackInspector.cc | 38 ++++++-----
Processes/StackInspector/StackInspector.h | 16 +++--
Setup/SetupTrajectory.h | 14 +++-
Stack/SuperStupidStack/SuperStupidStack.h | 7 ++
17 files changed, 196 insertions(+), 113 deletions(-)
diff --git a/Documentation/Examples/CMakeLists.txt b/Documentation/Examples/CMakeLists.txt
index 237e7a0b..7faa2f43 100644
--- a/Documentation/Examples/CMakeLists.txt
+++ b/Documentation/Examples/CMakeLists.txt
@@ -18,6 +18,7 @@ add_executable (staticsequence_example staticsequence_example.cc)
target_link_libraries (staticsequence_example
CORSIKAprocesssequence
CORSIKAunits
+ CORSIKAgeometry
CORSIKAlogging)
install (TARGETS staticsequence_example DESTINATION share/examples)
diff --git a/Documentation/Examples/staticsequence_example.cc b/Documentation/Examples/staticsequence_example.cc
index 8049983f..9a29d6f7 100644
--- a/Documentation/Examples/staticsequence_example.cc
+++ b/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,9 +88,14 @@ void modular() {
const auto sequence = m1 + m2 + m3 + m4;
DummyData p;
- DummyTrajectory t;
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));
+
const int n = 100000000;
for (int i = 0; i < n; ++i) { sequence.DoContinuous(p, t, s); }
diff --git a/Framework/Cascade/CMakeLists.txt b/Framework/Cascade/CMakeLists.txt
index 5dbb843b..2d67e801 100644
--- a/Framework/Cascade/CMakeLists.txt
+++ b/Framework/Cascade/CMakeLists.txt
@@ -9,10 +9,12 @@ set (
set (
CORSIKAcascade_HEADERS
Cascade.h
+ #TrackingStep.h
)
#set (
# CORSIKAcascade_SOURCES
+# TrackingStep.cc
# Cascade.cc
# )
diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h
index 0e920cc4..6beee09c 100644
--- a/Framework/Cascade/Cascade.h
+++ b/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,33 @@ 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::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);
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;
};
diff --git a/Framework/Cascade/testCascade.cc b/Framework/Cascade/testCascade.cc
index edb825bc..3ade3948 100644
--- a/Framework/Cascade/testCascade.cc
+++ b/Framework/Cascade/testCascade.cc
@@ -15,6 +15,7 @@
#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
@@ -34,11 +35,11 @@ public:
ProcessSplit() {}
template <typename Particle>
- double MinStepLength(Particle&) const {
- return 0;
+ void MinStepLength(Particle&, Trajectory& ) const {
+ //return 0;
}
- template <typename Particle, typename Trajectory, typename Stack>
+ template <typename Particle, typename Stack>
EProcessReturn DoContinuous(Particle&, Trajectory&, Stack&) const {
// corsika::utls::ignore(p);
return EProcessReturn::eOk;
@@ -58,19 +59,34 @@ public:
void Init() { fCount = 0; }
- int GetCount() { return fCount; }
+ int GetCount() const { return fCount; }
private:
};
+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::Line traj(p.GetPosition(), v);
+ return corsika::geometry::Trajectory<corsika::geometry::Line>(traj, 0_s, 100_ns);
+ }
+};
+
+
TEST_CASE("Cascade", "[Cascade]") {
- stack_inspector::StackInspector<setup::Stack, setup::Trajectory> p0(true);
+ NewtonTracking<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();
diff --git a/Framework/Geometry/CMakeLists.txt b/Framework/Geometry/CMakeLists.txt
index 580066ae..78f0dee4 100644
--- a/Framework/Geometry/CMakeLists.txt
+++ b/Framework/Geometry/CMakeLists.txt
@@ -15,7 +15,7 @@ set (
BaseVector.h
QuantityVector.h
Trajectory.h
- BaseTrajectory.h
+ # BaseTrajectory.h
)
set (
diff --git a/Framework/Geometry/Helix.h b/Framework/Geometry/Helix.h
index b647688d..8eb21582 100644
--- a/Framework/Geometry/Helix.h
+++ b/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);
}
};
diff --git a/Framework/Geometry/Line.h b/Framework/Geometry/Line.h
index 829c1604..5de0ca9c 100644
--- a/Framework/Geometry/Line.h
+++ b/Framework/Geometry/Line.h
@@ -34,9 +34,9 @@ namespace corsika::geometry {
return r0 + v0 * t;
}
- LengthType DistanceBetween(corsika::units::si::TimeType t1,
- corsika::units::si::TimeType t2) const {
- assert(t2 >= t1);
+ LengthType GetDistanceBetween(corsika::units::si::TimeType t1,
+ corsika::units::si::TimeType t2) const {
+ // assert(t2 >= t1);
return v0.norm() * (t2 - t1);
}
};
diff --git a/Framework/Geometry/Trajectory.h b/Framework/Geometry/Trajectory.h
index 97adfdd3..66e77a47 100644
--- a/Framework/Geometry/Trajectory.h
+++ b/Framework/Geometry/Trajectory.h
@@ -12,37 +12,43 @@
#ifndef _include_TRAJECTORY_H
#define _include_TRAJECTORY_H
-#include <corsika/geometry/BaseTrajectory.h>
+//#include <corsika/geometry/BaseTrajectory.h>
#include <corsika/units/PhysicalUnits.h>
namespace corsika::geometry {
template <typename T>
- class Trajectory : public BaseTrajectory {
+ class Trajectory : public T { // BaseTrajectory {
- T fTraj;
+ // T fTraj;
+ corsika::units::si::TimeType fTStart, fTEnd;
public:
+ using T::GetPosition;
+ using T::GetDistanceBetween;
+
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) {}
+ : T(theT)
+ , fTStart(pTStart)
+ , fTEnd(pTEnd) {}
+ //: BaseTrajectory(pTStart, pTEnd)
+ // , fTraj(theT) {}
- 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));
+ return T::GetPosition(fTEnd * u + fTStart * (1 - u));
}
-
+
+ /*
LengthType GetDistance(corsika::units::si::TimeType t1,
- corsika::units::si::TimeType t2) const {
+ corsika::units::si::TimeType t2) const {
return fTraj.DistanceBetween(t1, t2);
}
-
-
+ */
};
} // namespace corsika::geometry
diff --git a/Framework/Geometry/testGeometry.cc b/Framework/Geometry/testGeometry.cc
index 742938ff..274ec7d5 100644
--- a/Framework/Geometry/testGeometry.cc
+++ b/Framework/Geometry/testGeometry.cc
@@ -154,7 +154,7 @@ TEST_CASE("Trajectories") {
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") {
@@ -182,6 +182,6 @@ TEST_CASE("Trajectories") {
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));
}
}
diff --git a/Framework/ProcessSequence/CMakeLists.txt b/Framework/ProcessSequence/CMakeLists.txt
index c2456782..9e95f0c8 100644
--- a/Framework/ProcessSequence/CMakeLists.txt
+++ b/Framework/ProcessSequence/CMakeLists.txt
@@ -42,6 +42,7 @@ add_executable (
target_link_libraries (
testProcessSequence
+ CORSIKAgeometry
CORSIKAprocesssequence
CORSIKAthirdparty # for catch2
)
diff --git a/Framework/ProcessSequence/ProcessSequence.h b/Framework/ProcessSequence/ProcessSequence.h
index 4188fb0a..6c1cdf3b 100644
--- a/Framework/ProcessSequence/ProcessSequence.h
+++ b/Framework/ProcessSequence/ProcessSequence.h
@@ -17,8 +17,12 @@
#include <corsika/process/DiscreteProcess.h>
#include <corsika/process/ProcessReturn.h>
+#include <corsika/setup/SetupTrajectory.h>
+
//#include <type_traits> // still needed ?
+using corsika::setup::Trajectory;
+
namespace corsika::process {
/* namespace detail { */
@@ -143,7 +147,7 @@ namespace corsika::process {
// example for a trait-based call:
// void Hello() const { detail::CallHello<T1,T2>::Call(A, B); }
- template <typename Particle, typename Trajectory, typename Stack>
+ template <typename Particle, typename Stack>
inline EProcessReturn DoContinuous(Particle& p, Trajectory& t, Stack& s) const {
EProcessReturn ret = EProcessReturn::eOk;
if constexpr (!std::is_base_of<DiscreteProcess<T1>, T1>::value) {
@@ -155,9 +159,10 @@ namespace corsika::process {
return ret;
}
- template <typename D>
- inline double MinStepLength(D& d) const {
- return std::min(A.MinStepLength(d), B.MinStepLength(d));
+ template <typename Particle>
+ inline void MinStepLength(Particle& p, Trajectory& step) const {
+ A.MinStepLength(p, step);
+ B.MinStepLength(p, step);
}
/*
diff --git a/Framework/ProcessSequence/testProcessSequence.cc b/Framework/ProcessSequence/testProcessSequence.cc
index 2d4b46b1..5e768083 100644
--- a/Framework/ProcessSequence/testProcessSequence.cc
+++ b/Framework/ProcessSequence/testProcessSequence.cc
@@ -19,17 +19,26 @@
#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::units::si;
+
using namespace std;
using namespace corsika::process;
+
+static const int nData = 10;
+
+
class ContinuousProcess1 : public ContinuousProcess<ContinuousProcess1> {
public:
ContinuousProcess1() {}
void Init() { cout << "ContinuousProcess1::Init" << endl; }
- template <typename D, typename T, typename S>
- inline EProcessReturn DoContinuous(D& d, T&, S&) const {
+ template <typename D, typename S>
+ inline EProcessReturn DoContinuous(D& d, Trajectory&, S&) const {
cout << "ContinuousProcess1::DoContinuous" << endl;
- for (int i = 0; i < 10; ++i) d.p[i] += 0.933;
+ for (int i = 0; i < nData; ++i) d.p[i] += 0.933;
return EProcessReturn::eOk;
}
@@ -44,10 +53,10 @@ class ContinuousProcess2 : public ContinuousProcess<ContinuousProcess2> {
public:
ContinuousProcess2() {}
void Init() { cout << "ContinuousProcess2::Init" << endl; }
- template <typename D, typename T, typename S>
- inline EProcessReturn DoContinuous(D& d, T&, S&) const {
+ template <typename D, typename S>
+ inline EProcessReturn DoContinuous(D& d, Trajectory&, S&) const {
cout << "ContinuousProcess2::DoContinuous" << endl;
- for (int i = 0; i < 20; ++i) d.p[i] += 0.933;
+ for (int i = 0; i < nData; ++i) d.p[i] += 0.933;
return EProcessReturn::eOk;
}
@@ -62,9 +71,9 @@ class Process1 : public DiscreteProcess<Process1> {
public:
Process1() {}
void Init() { cout << "Process1::Init" << endl; }
- 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] += 1 + i;
+ template <typename D, typename S>
+ inline EProcessReturn DoContinuous(D& d, Trajectory&, S&) const {
+ for (int i = 0; i < nData; ++i) d.p[i] += 1 + i;
return EProcessReturn::eOk;
}
template <typename Particle, typename Stack>
@@ -78,9 +87,9 @@ class Process2 : public DiscreteProcess<Process2> {
public:
Process2() {}
void Init() { cout << "Process2::Init" << endl; }
- 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.7;
+ 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;
}
template <typename Particle, typename Stack>
@@ -94,9 +103,9 @@ class Process3 : public DiscreteProcess<Process3> {
public:
Process3() {}
void Init() { cout << "Process3::Init" << endl; }
- 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.933;
+ 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;
}
template <typename Particle, typename Stack>
@@ -110,9 +119,9 @@ class Process4 : public BaseProcess<Process4> {
public:
Process4() {}
void Init() { cout << "Process4::Init" << endl; }
- 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] /= 1.2;
+ 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;
return EProcessReturn::eOk;
}
// inline double MinStepLength(D& d) {
@@ -120,10 +129,9 @@ public:
};
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 {};
-struct DummyTrajectory {};
TEST_CASE("Cascade", "[Cascade]") {
@@ -142,12 +150,18 @@ TEST_CASE("Cascade", "[Cascade]") {
const auto sequence2 = cp1 + m2 + m3 + cp2;
DummyData p;
- DummyTrajectory t;
DummyStack s;
-
+
cout << "-->init" << endl;
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);
@@ -155,11 +169,11 @@ TEST_CASE("Cascade", "[Cascade]") {
sequence.Init();
- const int n = 100;
- cout << "Running loop with n=" << n << endl;
- for (int i = 0; i < n; ++i) { sequence.DoContinuous(p, t, s); }
-
- for (int i = 0; i < 10; i++) { cout << "data[" << i << "]=" << p.p[i] << endl; }
+ 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 < nData; i++) { cout << "data[" << i << "]=" << p.p[i] << endl; }
+ cout << "done" << endl;
}
SECTION("sectionThree") {}
diff --git a/Processes/StackInspector/StackInspector.cc b/Processes/StackInspector/StackInspector.cc
index 17acb325..7480d12d 100644
--- a/Processes/StackInspector/StackInspector.cc
+++ b/Processes/StackInspector/StackInspector.cc
@@ -14,6 +14,8 @@
#include <corsika/logging/Logger.h>
+#include <corsika/setup/SetupTrajectory.h>
+
#include <iostream>
using namespace std;
@@ -21,18 +23,16 @@ using namespace corsika;
using namespace corsika::units::si;
using namespace corsika::process::stack_inspector;
-template <typename Stack, typename Trajectory>
-StackInspector<Stack, Trajectory>::StackInspector(const bool aReport)
+template <typename Stack>
+StackInspector<Stack>::StackInspector(const bool aReport)
: fReport(aReport) {}
-template <typename Stack, typename Trajectory>
-StackInspector<Stack, Trajectory>::~StackInspector() {}
-
-template <typename Stack, typename Trajectory>
-process::EProcessReturn StackInspector<Stack, Trajectory>::DoContinuous(Particle&,
- Trajectory&,
- Stack& s) const {
+template <typename Stack>
+StackInspector<Stack>::~StackInspector() {}
+template <typename Stack>
+process::EProcessReturn StackInspector<Stack>::DoContinuous(
+ Particle&, corsika::setup::Trajectory&, Stack& s) const {
static int countStep = 0;
if (!fReport) return EProcessReturn::eOk;
[[maybe_unused]] int i = 0;
@@ -40,22 +40,26 @@ process::EProcessReturn StackInspector<Stack, Trajectory>::DoContinuous(Particle
for (auto& iterP : s) {
EnergyType E = iterP.GetEnergy();
Etot += E;
+ cout << "i=" << setw(5) << fixed << (i++)
+ << ", id=" << setw(30) << iterP.GetPID()
+ << " E=" << setw(15) << scientific << (E / 1_GeV) << " GeV, "
+ //<< " pos=" << iterP.GetPosition()
+ << endl;
}
countStep++;
cout << countStep << " " << s.GetSize() << " " << Etot / 1_GeV << " " << endl;
return EProcessReturn::eOk;
}
-template <typename Stack, typename Trajectory>
-double StackInspector<Stack, Trajectory>::MinStepLength(Particle&) const {
- return 0;
+template <typename Stack>
+void StackInspector<Stack>::MinStepLength(Particle&,
+ corsika::setup::Trajectory&) const {
+ // return 0;
}
-template <typename Stack, typename Trajectory>
-void StackInspector<Stack, Trajectory>::Init() {}
+template <typename Stack>
+void StackInspector<Stack>::Init() {}
#include <corsika/setup/SetupStack.h>
-#include <corsika/setup/SetupTrajectory.h>
-template class corsika::process::stack_inspector::StackInspector<setup::Stack,
- setup::Trajectory>;
+template class corsika::process::stack_inspector::StackInspector<setup::Stack>;
diff --git a/Processes/StackInspector/StackInspector.h b/Processes/StackInspector/StackInspector.h
index 118c06d3..9d8d14f6 100644
--- a/Processes/StackInspector/StackInspector.h
+++ b/Processes/StackInspector/StackInspector.h
@@ -14,13 +14,19 @@
#include <corsika/process/ContinuousProcess.h>
+#include <corsika/setup/SetupTrajectory.h>
+
+//namespace corsika::setup {
+//class BaseTrajectory;
+//}
+
namespace corsika::process {
namespace stack_inspector {
- template <typename Stack, typename Trajectory>
+ template <typename Stack>
class StackInspector
- : public corsika::process::ContinuousProcess<StackInspector<Stack, Trajectory>> {
+ : public corsika::process::ContinuousProcess<StackInspector<Stack>> {
typedef typename Stack::ParticleType Particle;
@@ -31,10 +37,10 @@ namespace corsika::process {
void Init();
// template <typename Particle, typename Trajectory, typename Stack>
- EProcessReturn DoContinuous(Particle&, Trajectory&, Stack& s) const;
+ EProcessReturn DoContinuous(Particle&, corsika::setup::Trajectory&, Stack& s) const;
- // template <typename Particle>
- double MinStepLength(Particle&) const;
+ // template <typename Particle>
+ void MinStepLength(Particle&, corsika::setup::Trajectory&) const;
private:
bool fReport;
diff --git a/Setup/SetupTrajectory.h b/Setup/SetupTrajectory.h
index d7f5c482..2950929c 100644
--- a/Setup/SetupTrajectory.h
+++ b/Setup/SetupTrajectory.h
@@ -12,12 +12,22 @@
#ifndef _corsika_setup_setuptrajectory_h_
#define _corsika_setup_setuptrajectory_h_
+#include <corsika/geometry/Helix.h>
#include <corsika/geometry/Line.h>
#include <corsika/geometry/Trajectory.h>
+#include <variant>
+
namespace corsika::setup {
- typedef corsika::geometry::Trajectory<corsika::geometry::Line> Trajectory;
-}
+ using corsika::geometry::Helix;
+ using corsika::geometry::Line;
+
+ typedef std::variant<std::monostate,
+ corsika::geometry::Trajectory<Line>,
+ corsika::geometry::Trajectory<Helix>>
+ Trajectory;
+
+} // namespace corsika::setup
#endif
diff --git a/Stack/SuperStupidStack/SuperStupidStack.h b/Stack/SuperStupidStack/SuperStupidStack.h
index 816a8c5b..fa6bf4e2 100644
--- a/Stack/SuperStupidStack/SuperStupidStack.h
+++ b/Stack/SuperStupidStack/SuperStupidStack.h
@@ -30,6 +30,7 @@ namespace corsika::stack {
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;
@@ -62,6 +63,12 @@ namespace corsika::stack {
MomentumVector GetMomentum() const { return GetStackData().GetMomentum(GetIndex()); }
Point GetPosition() const { return GetStackData().GetPosition(GetIndex()); }
TimeType GetTime() const { return GetStackData().GetTime(GetIndex()); }
+
+#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); }
+
};
/**
--
GitLab