From 15fbf4144c629ddc25bae63d8f9180b3d4d8d257 Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Tue, 11 Dec 2018 09:50:24 +0100
Subject: [PATCH] this mostly resolves it
---
Documentation/Examples/cascade_example.cc | 18 +--
Framework/Cascade/Cascade.h | 48 ++++++--
Framework/Cascade/testCascade.cc | 43 ++++---
Framework/ProcessSequence/BaseProcess.h | 19 ++-
Framework/ProcessSequence/DiscreteProcess.h | 6 +
Framework/ProcessSequence/ProcessReturn.h | 3 +-
Framework/ProcessSequence/ProcessSequence.h | 110 +++++++++++++++---
.../ProcessSequence/testProcessSequence.cc | 24 ++++
Framework/StackInterface/Stack.h | 3 +
Processes/StackInspector/StackInspector.cc | 5 +-
Processes/StackInspector/StackInspector.h | 2 +-
Stack/SuperStupidStack/SuperStupidStack.h | 2 +-
12 files changed, 231 insertions(+), 52 deletions(-)
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index c22fa23b..09ae4815 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -36,13 +36,12 @@ using namespace std;
static int fCount = 0;
-class ProcessSplit : public corsika::process::BaseProcess<ProcessSplit> {
+class ProcessSplit : public corsika::process::DiscreteProcess<ProcessSplit> {
public:
ProcessSplit() {}
template <typename Particle, typename Track>
- double MinStepLength(Particle& p, Track&) const {
-
+ double GetInteractionLength(Particle& p, Track&) const {
// beam particles for sibyll : 1, 2, 3 for p, pi, k
// read from cross section code table
int kBeam = 1;
@@ -81,17 +80,20 @@ public:
std::cout << "ProcessSplit: "
<< "interaction length (g/cm2): " << int_length << std::endl;
// add exponential sampling
- int a = 0;
- const double next_step = -int_length * log(s_rndm_(a));
/*
what are the units of the output? slant depth or 3space length?
*/
- std::cout << "ProcessSplit: "
- << "next step (g/cm2): " << next_step << std::endl;
- return next_step;
+ return int_length;
+ //
+ //int a = 0;
+ //const double next_step = -int_length * log(s_rndm_(a));
+ //std::cout << "ProcessSplit: "
+ // << "next step (g/cm2): " << next_step << std::endl;
+ //return next_step;
}
+
template <typename Particle, typename Track, typename Stack>
EProcessReturn DoContinuous(Particle&, Track&, Stack&) const {
// corsika::utls::ignore(p);
diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h
index 592372f5..a1e6ff8a 100644
--- a/Framework/Cascade/Cascade.h
+++ b/Framework/Cascade/Cascade.h
@@ -14,11 +14,12 @@
#include <corsika/process/ProcessReturn.h>
#include <corsika/units/PhysicalUnits.h>
+#include <corsika/setup/SetupTrajectory.h>
+#include <corsika/random/RNGManager.h>
#include <type_traits>
-#include <corsika/setup/SetupTrajectory.h>
-
+using namespace corsika;
using namespace corsika::units::si;
namespace corsika::cascade {
@@ -50,24 +51,53 @@ namespace corsika::cascade {
}
// do cascade equations, which can put new particles on Stack,
// thus, the double loop
- // DoCascadeEquations(); //
+ // DoCascadeEquations();
}
}
- void Step(Particle& particle) {
+ void Step(Particle& particle) {
+
corsika::setup::Trajectory step = fTracking.GetTrack(particle);
- fProcesseList.MinStepLength(particle, step);
-
+ const double total_inv_lambda = fProcesseList.GetTotalInverseInteractionLength(particle, step);
+
+ // sample random exponential step length
+ // ....
+ static corsika::random::RNG& rmng =
+ corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");
+ const double sample_step = rmng() / (double)rmng.max();
+ const double next_step = -log(sample_step) / total_inv_lambda;
+
+ const double min_step = fProcesseList.MinStepLength(particle, step);
+ // convert next_step from grammage to length
+ // Environment::GetDistance(step, next_step);
+ // ....
+
+ // update step with actual min(min_step, next_step)
+ // ....
+
/// here the particle is actually moved along the trajectory to new position:
// std::visit(corsika::setup::ParticleUpdate<Particle>{particle}, step);
particle.SetPosition(step.GetPosition(1));
-
+
corsika::process::EProcessReturn status =
fProcesseList.DoContinuous(particle, step, fStack);
+
if (status == corsika::process::EProcessReturn::eParticleAbsorbed) {
- fStack.Delete(particle); // TODO: check if this is really needed
+ // fStack.Delete(particle); // TODO: check if this is really needed
} else {
- fProcesseList.DoDiscrete(particle, fStack);
+
+ // check if min_step < random_step -> skip discrete if rndm>min_step/random_step
+ if (min_step<next_step) {
+ const double p_next = min_step/next_step;
+ const double sample_skip = rmng() / (double)rmng.max();
+ if (sample_skip>p_next) {
+ return;// corsika::process::EProcessReturn::eOk;
+ }
+ }
+ const double sample_process = rmng() / (double)rmng.max();
+ double inv_lambda_count = 0;
+ fProcesseList.SelectDiscrete(particle, fStack, total_inv_lambda,
+ sample_process, inv_lambda_count);
}
}
diff --git a/Framework/Cascade/testCascade.cc b/Framework/Cascade/testCascade.cc
index a218cb78..f315cc29 100644
--- a/Framework/Cascade/testCascade.cc
+++ b/Framework/Cascade/testCascade.cc
@@ -17,6 +17,8 @@
#include <corsika/stack/super_stupid/SuperStupidStack.h>
+#include <corsika/particles/ParticleProperties.h>
+
#include <corsika/geometry/Point.h>
#include <corsika/geometry/RootCoordinateSystem.h>
#include <corsika/geometry/Vector.h>
@@ -42,15 +44,18 @@ static int fCount = 0;
class ProcessSplit : public corsika::process::BaseProcess<ProcessSplit> {
public:
ProcessSplit() {}
-
- template <typename Particle>
- void MinStepLength(Particle&, setup::Trajectory&) const {}
-
- template <typename Particle, typename Stack>
- EProcessReturn DoContinuous(Particle&, setup::Trajectory&, Stack&) const {
+
+ template <typename Particle, typename T>
+ double MinStepLength(Particle&, T&) const { return 1; }
+
+ template <typename Particle, typename T>
+ double GetInteractionLength(Particle&, T&) const { return 1; }
+
+ template <typename Particle, typename T, typename Stack>
+ EProcessReturn DoContinuous(Particle&, T&, Stack&) const {
return EProcessReturn::eOk;
}
-
+
template <typename Particle, typename Stack>
void DoDiscrete(Particle& p, Stack& s) const {
EnergyType E = p.GetEnergy();
@@ -60,7 +65,9 @@ public:
} else {
p.SetEnergy(E / 2);
auto pnew = s.NewParticle();
- // s.Copy(p, pnew);
+ // s.Copy(p, pnew); fix that .... todo
+ pnew.SetPID(p.GetPID());
+ pnew.SetTime(p.GetTime());
pnew.SetEnergy(E / 2);
pnew.SetPosition(p.GetPosition());
pnew.SetMomentum(p.GetMomentum());
@@ -76,27 +83,34 @@ private:
TEST_CASE("Cascade", "[Cascade]") {
- tracking_line::TrackingLine<setup::Stack> tracking;
+ corsika::random::RNGManager& rmng = corsika::random::RNGManager::GetInstance();
+ ;
+ const std::string str_name = "s_rndm";
+ rmng.RegisterRandomStream(str_name);
+ 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(tracking, sequence, stack);
-
CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
-
+
stack.Clear();
auto particle = stack.NewParticle();
EnergyType E0 = 100_GeV;
+ particle.SetPID(particles::Code::Electron);
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}));
+ particle.SetTime(0_ns);
EAS.Init();
EAS.Run();
-
+
+ /*
SECTION("sectionTwo") {
for (int i = 0; i < 0; ++i) {
stack.Clear();
@@ -105,8 +119,9 @@ TEST_CASE("Cascade", "[Cascade]") {
particle.SetEnergy(E0);
EAS.Init();
EAS.Run();
-
+
// cout << "Result: E0=" << E0 / 1_GeV << "GeV, count=" << p1.GetCount() << endl;
}
}
+ */
}
diff --git a/Framework/ProcessSequence/BaseProcess.h b/Framework/ProcessSequence/BaseProcess.h
index 34dcb198..b17305d5 100644
--- a/Framework/ProcessSequence/BaseProcess.h
+++ b/Framework/ProcessSequence/BaseProcess.h
@@ -29,14 +29,27 @@ namespace corsika::process {
struct BaseProcess {
derived& GetRef() { return static_cast<derived&>(*this); }
const derived& GetRef() const { return static_cast<const derived&>(*this); }
-
+
template <typename Particle, typename Stack>
inline EProcessReturn DoDiscrete(Particle&, Stack&) const; // {}
-
+
template <typename Particle, typename Track, typename Stack>
inline EProcessReturn DoContinuous(Particle&, Track&, Stack&) const; // {}
- };
+ template <typename Particle, typename Track>
+ inline double GetInverseInteractionLength(Particle& p, Track& t) const {
+ return 1./GetRef().GetInteractionLength(p, t);
+ }
+};
+
+ /*
+ template<template<typename, typename> class T, typename A, typename B>
+ typename BaseProcess< T<A, B> >::is_process_sequence
+ {
+ static const bool value = true;
+ };
+ */
+
/*
template <typename T>
struct is_base {
diff --git a/Framework/ProcessSequence/DiscreteProcess.h b/Framework/ProcessSequence/DiscreteProcess.h
index a540c92c..d0464e4e 100644
--- a/Framework/ProcessSequence/DiscreteProcess.h
+++ b/Framework/ProcessSequence/DiscreteProcess.h
@@ -37,6 +37,12 @@ namespace corsika::process {
// -> enforce derived to implement DoDiscrete...
template <typename Particle, typename Stack>
inline EProcessReturn DoDiscrete(Particle&, Stack&) const; // {}
+
+ template <typename Particle, typename Track>
+ inline double GetInverseInteractionLength(Particle& p, Track& t) const {
+ return 1./GetRef().GetInteractionLength(p, t);
+ }
+
};
} // namespace corsika::process
diff --git a/Framework/ProcessSequence/ProcessReturn.h b/Framework/ProcessSequence/ProcessReturn.h
index 82cc816c..039a5622 100644
--- a/Framework/ProcessSequence/ProcessReturn.h
+++ b/Framework/ProcessSequence/ProcessReturn.h
@@ -23,7 +23,8 @@ namespace corsika::process {
enum class EProcessReturn {
eOk = 1,
eParticleAbsorbed = 2,
- };
+ eInteracted = 3,
+ };
} // namespace corsika::process
#endif
diff --git a/Framework/ProcessSequence/ProcessSequence.h b/Framework/ProcessSequence/ProcessSequence.h
index 54ac90e7..a24356d8 100644
--- a/Framework/ProcessSequence/ProcessSequence.h
+++ b/Framework/ProcessSequence/ProcessSequence.h
@@ -17,6 +17,9 @@
#include <corsika/process/DiscreteProcess.h>
#include <corsika/process/ProcessReturn.h>
+#include <limits>
+#include <cmath>
+
//#include <corsika/setup/SetupTrajectory.h>
// using corsika::setup::Trajectory;
//#include <variant>
@@ -25,7 +28,7 @@
namespace corsika::process {
- /* namespace detail { */
+ // namespace detail {
/* /\* template<typename TT1, typename TT2, typename Type = void> *\/ */
/* /\* struct CallHello { *\/ */
@@ -41,7 +44,7 @@ namespace corsika::process {
/* /\* static void Call(const TT1&, const TT2&) { *\/ */
/* /\* std::cout << "special" << std::endl; *\/ */
/* /\* } *\/ */
- /* /\* }; *\/ */
+ /* }; */
/* template<typename T1, typename T2, typename Particle, typename Trajectory, typename
* Stack> //, typename Type = void> */
@@ -121,7 +124,7 @@ namespace corsika::process {
/* } */
/* }; */
/* *\/ */
- /* } // end namespace detail */
+ //} // end namespace detail
/**
\class ProcessSequence
@@ -134,8 +137,24 @@ namespace corsika::process {
https://en.wikipedia.org/wiki/Curiously_recurring_template_pattern
*/
+ //template <typename T1, typename T2>
+ //class ProcessSequence : public BaseProcess<ProcessSequence<T1, T2> >;
+
+ // this is a marker to track which BaseProcess is also a ProcessSequence
+ template <typename T>
+ struct is_process_sequence {
+ static const bool value = false;
+ };
+
+
template <typename T1, typename T2>
class ProcessSequence : public BaseProcess<ProcessSequence<T1, T2> > {
+
+ /* template <>
+ struct BaseProcess<ProcessSequence<T1, T2> >::is_process_sequence<true> {
+ static const bool value = true;
+ };*/
+
public:
const T1& A;
const T2& B;
@@ -162,20 +181,38 @@ namespace corsika::process {
}
template <typename Particle, typename Track>
- inline void MinStepLength(Particle& p, Track& track) const {
- A.MinStepLength(p, track);
- B.MinStepLength(p, track);
+ inline double MinStepLength(Particle& p, Track& track) const {
+ double min_length = std::numeric_limits<double>::infinity();
+ if constexpr (!std::is_base_of<DiscreteProcess<T1>, T1>::value) {
+ min_length = std::min(min_length, A.MinStepLength(p,track));
+ }
+ if constexpr (!std::is_base_of<DiscreteProcess<T2>, T2>::value) {
+ min_length = std::min(min_length, B.MinStepLength(p,track));
+ }
+ return min_length;
+ }
+
+ template <typename Particle, typename Track>
+ inline double GetTotalInteractionLength(Particle& p, Track& t) const {
+ return 1. / GetInverseInteractionLength(p, t);
+ }
+
+ template <typename Particle, typename Track>
+ inline double GetTotalInverseInteractionLength(Particle& p, Track& t) const {
+ return GetInverseInteractionLength(p, t);
}
- /*
template <typename Particle, typename Track>
- inline Track Transport(Particle& p, double& length) const {
- A.Transport(p, length); // todo: maybe check (?) if there is more than one Transport
- // process implemented??
- return B.Transport(
- p, length); // need to do this also to decide which Track to return!!!!
+ inline double GetInverseInteractionLength(Particle& p, Track& t) const {
+ double tot = 0;
+ if constexpr (!std::is_base_of<ContinuousProcess<T1>, T1>::value) {
+ tot += A.GetInverseInteractionLength(p, t);
+ }
+ if constexpr (!std::is_base_of<ContinuousProcess<T2>, T2>::value) {
+ tot += B.GetInverseInteractionLength(p, t);
+ }
+ return tot;
}
- */
template <typename Particle, typename Stack>
inline EProcessReturn DoDiscrete(Particle& p, Stack& s) const {
@@ -188,6 +225,47 @@ namespace corsika::process {
return EProcessReturn::eOk;
}
+ template <typename Particle, typename Stack>
+ inline EProcessReturn SelectDiscrete(Particle& p, Stack& s,
+ const double lambda_inv_tot,
+ const double rndm_select,
+ double& lambda_inv_count) const {
+ if constexpr (is_process_sequence<T1>::value) {
+ // if A is a process sequence --> check inside
+ const EProcessReturn ret = A.SelectDiscrete(p, s, lambda_inv_count, rndm_select, lambda_inv_count);
+ // if A did suceed, stop routine
+ if (ret != EProcessReturn::eOk) {
+ return ret;
+ }
+ } else if constexpr (!std::is_base_of<ContinuousProcess<T1>, T1>::value) {
+ // if this is not a ContinuousProcess --> evaluate probability
+ lambda_inv_count += A.GetInverseInteractionLength(p, s);
+ // check if we should execute THIS process and then EXIT
+ if (rndm_select<lambda_inv_count/lambda_inv_tot) {
+ A.DoDiscrete(p, s);
+ return EProcessReturn::eInteracted;
+ }
+ } // end branch A
+
+ if constexpr (is_process_sequence<T2>::value) {
+ // if A is a process sequence --> check inside
+ const EProcessReturn ret = B.SelectDiscrete(p, s, lambda_inv_count, rndm_select, lambda_inv_count);
+ // if A did suceed, stop routine
+ if (ret != EProcessReturn::eOk) {
+ return ret;
+ }
+ } else if constexpr (!std::is_base_of<ContinuousProcess<T2>, T2>::value) {
+ // if this is not a ContinuousProcess --> evaluate probability
+ lambda_inv_count += B.GetInverseInteractionLength(p, s);
+ // check if we should execute THIS process and then EXIT
+ if (rndm_select<lambda_inv_count/lambda_inv_tot) {
+ B.DoDiscrete(p, s);
+ return EProcessReturn::eInteracted;
+ }
+ } // end branch A
+ return EProcessReturn::eOk;
+ }
+
/// TODO the const_cast is not nice, think about the constness here
inline void Init() const {
const_cast<T1*>(&A)->Init();
@@ -289,6 +367,12 @@ namespace corsika::process {
}
*/
+ template<template<typename, typename> class T, typename A, typename B>
+ struct is_process_sequence< T<A,B> >
+ {
+ static const bool value = true;
+ };
+
} // namespace corsika::process
#endif
diff --git a/Framework/ProcessSequence/testProcessSequence.cc b/Framework/ProcessSequence/testProcessSequence.cc
index 7094cacb..008201b7 100644
--- a/Framework/ProcessSequence/testProcessSequence.cc
+++ b/Framework/ProcessSequence/testProcessSequence.cc
@@ -100,6 +100,11 @@ public:
cout << "Process2::DoDiscrete" << endl;
return EProcessReturn::eOk;
}
+ template <typename Particle, typename Track>
+ inline double GetInteractionLength(Particle&, Track&) const {
+ cout << "Process2::GetInteractionLength" << endl;
+ return 3;
+ }
};
class Process3 : public DiscreteProcess<Process3> {
@@ -118,6 +123,11 @@ public:
cout << "Process3::DoDiscrete" << endl;
return EProcessReturn::eOk;
}
+ template <typename Particle, typename Track>
+ inline double GetInteractionLength(Particle&, Track&) const {
+ cout << "Process3::GetInteractionLength" << endl;
+ return 1.;
+ }
};
class Process4 : public BaseProcess<Process4> {
@@ -169,6 +179,20 @@ TEST_CASE("Process Sequence", "[Process Sequence]") {
// REQUIRE_THROWS(sequence_wrong.Init());
}
+ SECTION("interaction length") {
+ ContinuousProcess1 cp1(0);
+ Process2 m2(1);
+ Process3 m3(2);
+
+ DummyStack s;
+ DummyTrajectory t;
+
+ const auto sequence2 = cp1 + m2 + m3;
+ double tot = sequence2.GetTotalInteractionLength(s, t);
+ double tot_inv = sequence2.GetTotalInverseInteractionLength(s, t);
+ cout << "lambda_tot=" << tot << " lambda_tot_inv=" << tot_inv << endl;
+ }
+
SECTION("sectionTwo") {
ContinuousProcess1 cp1(0);
diff --git a/Framework/StackInterface/Stack.h b/Framework/StackInterface/Stack.h
index b4660144..e71880cc 100644
--- a/Framework/StackInterface/Stack.h
+++ b/Framework/StackInterface/Stack.h
@@ -67,6 +67,9 @@ namespace corsika::stack {
IncrementSize();
return StackIterator(*this, GetSize() - 1);
}
+ void Copy(StackIterator& a, StackIterator& b) {
+ Copy(a.GetIndex(), b.GetIndex());
+ }
/// delete this particle
void Delete(StackIterator& p) {
if (GetSize() == 0) { /*error*/
diff --git a/Processes/StackInspector/StackInspector.cc b/Processes/StackInspector/StackInspector.cc
index b0acefc4..abf47afd 100644
--- a/Processes/StackInspector/StackInspector.cc
+++ b/Processes/StackInspector/StackInspector.cc
@@ -17,6 +17,7 @@
#include <corsika/setup/SetupTrajectory.h>
+#include <limits>
#include <iostream>
using namespace std;
@@ -56,8 +57,8 @@ process::EProcessReturn StackInspector<Stack>::DoContinuous(Particle&, setup::Tr
}
template <typename Stack>
-void StackInspector<Stack>::MinStepLength(Particle&, setup::Trajectory&) const {
- // return 0;
+double StackInspector<Stack>::MinStepLength(Particle&, setup::Trajectory&) const {
+ return std::numeric_limits<double>::infinity();
}
template <typename Stack>
diff --git a/Processes/StackInspector/StackInspector.h b/Processes/StackInspector/StackInspector.h
index 7b68d92d..38d5cfc0 100644
--- a/Processes/StackInspector/StackInspector.h
+++ b/Processes/StackInspector/StackInspector.h
@@ -36,7 +36,7 @@ namespace corsika::process {
EProcessReturn DoContinuous(Particle&, corsika::setup::Trajectory&, Stack& s) const;
// template <typename Particle>
- void MinStepLength(Particle&, corsika::setup::Trajectory&) const;
+ double MinStepLength(Particle&, corsika::setup::Trajectory&) const;
private:
bool fReport;
diff --git a/Stack/SuperStupidStack/SuperStupidStack.h b/Stack/SuperStupidStack/SuperStupidStack.h
index d1c8bc2b..8c799344 100644
--- a/Stack/SuperStupidStack/SuperStupidStack.h
+++ b/Stack/SuperStupidStack/SuperStupidStack.h
@@ -128,7 +128,7 @@ namespace corsika::stack {
void Swap(const int i1, const int i2) {
std::swap(fDataPID[i2], fDataPID[i1]);
std::swap(fDataE[i2], fDataE[i1]);
- std::swap(fMomentum[i2], fMomentum[i1]); // should be Momentum !!!!
+ std::swap(fMomentum[i2], fMomentum[i1]);
std::swap(fPosition[i2], fPosition[i1]);
std::swap(fTime[i2], fTime[i1]);
}
--
GitLab