From b12a7b73dd9cf9a849c8fdee06ded1a90325b365 Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Sun, 2 Dec 2018 22:36:17 +0100
Subject: [PATCH] first moving particles in corsika
---
Framework/Cascade/Cascade.h | 18 +++-------
Framework/Cascade/testCascade.cc | 38 ++++++++++++++++-----
Framework/Geometry/Trajectory.h | 36 +++++++++----------
Framework/Geometry/testGeometry.cc | 4 +--
Framework/ProcessSequence/ProcessSequence.h | 2 ++
Processes/StackInspector/StackInspector.cc | 2 +-
Setup/SetupTrajectory.h | 33 ++++++++++++++++--
Stack/SuperStupidStack/SuperStupidStack.h | 2 +-
8 files changed, 89 insertions(+), 46 deletions(-)
diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h
index 6beee09c..c97e3376 100644
--- a/Framework/Cascade/Cascade.h
+++ b/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) {
diff --git a/Framework/Cascade/testCascade.cc b/Framework/Cascade/testCascade.cc
index 6d6258d2..0cb36f52 100644
--- a/Framework/Cascade/testCascade.cc
+++ b/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();
diff --git a/Framework/Geometry/Trajectory.h b/Framework/Geometry/Trajectory.h
index 66e77a47..99b02b37 100644
--- a/Framework/Geometry/Trajectory.h
+++ b/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
diff --git a/Framework/Geometry/testGeometry.cc b/Framework/Geometry/testGeometry.cc
index 874c2ce0..c582b2d6 100644
--- a/Framework/Geometry/testGeometry.cc
+++ b/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());
diff --git a/Framework/ProcessSequence/ProcessSequence.h b/Framework/ProcessSequence/ProcessSequence.h
index 6c1cdf3b..f4e4eb80 100644
--- a/Framework/ProcessSequence/ProcessSequence.h
+++ b/Framework/ProcessSequence/ProcessSequence.h
@@ -19,6 +19,8 @@
#include <corsika/setup/SetupTrajectory.h>
+#include <variant>
+
//#include <type_traits> // still needed ?
using corsika::setup::Trajectory;
diff --git a/Processes/StackInspector/StackInspector.cc b/Processes/StackInspector/StackInspector.cc
index 94abaff7..546b578b 100644
--- a/Processes/StackInspector/StackInspector.cc
+++ b/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++;
diff --git a/Setup/SetupTrajectory.h b/Setup/SetupTrajectory.h
index 2950929c..bad174c4 100644
--- a/Setup/SetupTrajectory.h
+++ b/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
diff --git a/Stack/SuperStupidStack/SuperStupidStack.h b/Stack/SuperStupidStack/SuperStupidStack.h
index 27e97f4a..c157a6f6 100644
--- a/Stack/SuperStupidStack/SuperStupidStack.h
+++ b/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); }
};
--
GitLab