From 91d6b637c24267f9a899d9333b31c8c453eec77a Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Mon, 10 Dec 2018 08:32:18 +0100
Subject: [PATCH 1/3] switch on -O3 for Release build
---
CMakeLists.txt | 3 +++
1 file changed, 3 insertions(+)
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 2a5cae812..11d785415 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -38,6 +38,9 @@ if(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES)
"Debug" "Release" "MinSizeRel" "RelWithDebInfo")
endif()
+#set(CMAKE_CXX_FLAGS "-Wall -Wextra")
+set(CMAKE_CXX_FLAGS_DEBUG "-O0 -g")
+set(CMAKE_CXX_FLAGS_RELEASE "-O3")
# unit testing coverage, does not work yet
#include (CodeCoverage)
--
GitLab
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 2/3] 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 c22fa23bb..09ae4815c 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 592372f5d..a1e6ff8a2 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 a218cb788..f315cc292 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 34dcb198e..b17305d5a 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 a540c92c7..d0464e4eb 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 82cc816c5..039a56222 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 54ac90e7a..a24356d8e 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 7094cacb9..008201b7e 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 b4660144d..e71880ccf 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 b0acefc4a..abf47afd1 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 7b68d92d3..38d5cfc01 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 d1c8bc2b3..8c7993442 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
From beca7928b9cf8ff25332931076ba1a0274a396f1 Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Wed, 12 Dec 2018 00:55:05 +0100
Subject: [PATCH 3/3] added new test section, final work
---
Documentation/Examples/cascade_example.cc | 19 +-
Framework/Cascade/Cascade.h | 97 ++++--
Framework/Cascade/testCascade.cc | 30 +-
Framework/ProcessSequence/BaseProcess.h | 16 +-
Framework/ProcessSequence/CMakeLists.txt | 3 +-
Framework/ProcessSequence/DecayProcess.h | 51 +++
Framework/ProcessSequence/DiscreteProcess.h | 8 +-
.../ProcessSequence/InteractionProcess.h | 51 +++
Framework/ProcessSequence/ProcessReturn.h | 3 +-
Framework/ProcessSequence/ProcessSequence.h | 296 ++++++++----------
.../ProcessSequence/testProcessSequence.cc | 59 +++-
Framework/StackInterface/Stack.h | 4 +-
Processes/StackInspector/StackInspector.cc | 4 +-
Processes/StackInspector/StackInspector.h | 2 +-
14 files changed, 398 insertions(+), 245 deletions(-)
create mode 100644 Framework/ProcessSequence/DecayProcess.h
create mode 100644 Framework/ProcessSequence/InteractionProcess.h
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index 09ae4815c..d7e673ee5 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -36,7 +36,7 @@ using namespace std;
static int fCount = 0;
-class ProcessSplit : public corsika::process::DiscreteProcess<ProcessSplit> {
+class ProcessSplit : public corsika::process::InteractionProcess<ProcessSplit> {
public:
ProcessSplit() {}
@@ -86,14 +86,13 @@ public:
*/
return int_length;
//
- //int a = 0;
- //const double next_step = -int_length * log(s_rndm_(a));
- //std::cout << "ProcessSplit: "
+ // 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;
+ // return next_step;
}
-
template <typename Particle, typename Track, typename Stack>
EProcessReturn DoContinuous(Particle&, Track&, Stack&) const {
// corsika::utls::ignore(p);
@@ -101,8 +100,8 @@ public:
}
template <typename Particle, typename Stack>
- void DoDiscrete(Particle& p, Stack& s) const {
- cout << "DoDiscrete: " << p.GetPID() << " interaction? "
+ void DoInteraction(Particle& p, Stack& s) const {
+ cout << "DoInteraction: " << p.GetPID() << " interaction? "
<< process::sibyll::CanInteract(p.GetPID()) << endl;
if (process::sibyll::CanInteract(p.GetPID())) {
@@ -127,11 +126,11 @@ public:
int kTarget = 1; // p.GetPID();
std::cout << "ProcessSplit: "
- << " DoDiscrete: E(GeV):" << E / 1_GeV << " Ecm(GeV): " << Ecm / 1_GeV
+ << " DoInteraction: E(GeV):" << E / 1_GeV << " Ecm(GeV): " << Ecm / 1_GeV
<< std::endl;
if (E < 8.5_GeV || Ecm < 10_GeV) {
std::cout << "ProcessSplit: "
- << " DoDiscrete: dropping particle.." << std::endl;
+ << " DoInteraction: dropping particle.." << std::endl;
p.Delete();
fCount++;
} else {
diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h
index a1e6ff8a2..7cd9d077d 100644
--- a/Framework/Cascade/Cascade.h
+++ b/Framework/Cascade/Cascade.h
@@ -13,9 +13,9 @@
#define _include_corsika_cascade_Cascade_h_
#include <corsika/process/ProcessReturn.h>
-#include <corsika/units/PhysicalUnits.h>
-#include <corsika/setup/SetupTrajectory.h>
#include <corsika/random/RNGManager.h>
+#include <corsika/setup/SetupTrajectory.h>
+#include <corsika/units/PhysicalUnits.h>
#include <type_traits>
@@ -55,49 +55,88 @@ namespace corsika::cascade {
}
}
- void Step(Particle& particle) {
-
+ void Step(Particle& particle) {
+
+ // get access to random number generator
+ static corsika::random::RNG& rmng =
+ corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");
+
+ // determine geometric tracking
corsika::setup::Trajectory step = fTracking.GetTrack(particle);
- const double total_inv_lambda = fProcesseList.GetTotalInverseInteractionLength(particle, step);
-
+
+ // determine combined total interaction length (inverse)
+ 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 next_interact = -log(sample_step) / total_inv_lambda;
+ std::cout << "total_inv_lambda=" << total_inv_lambda
+ << ", next_interact=" << next_interact << std::endl;
+
+ // determine the maximum geometric step length
+ const double distance_max = fProcesseList.MaxStepLength(particle, step);
+ std::cout << "distance_max=" << distance_max << std::endl;
+
+ // determine combined total inverse decay time
+ const double total_inv_lifetime = fProcesseList.GetTotalInverseLifetime(particle);
- const double min_step = fProcesseList.MinStepLength(particle, step);
- // convert next_step from grammage to length
+ // sample random exponential decay time
+ const double sample_decay = rmng() / (double)rmng.max();
+ const double next_decay = -log(sample_decay) / total_inv_lifetime;
+ std::cout << "total_inv_lifetime=" << total_inv_lifetime
+ << ", next_decay=" << next_decay << std::endl;
+
+ // convert next_step from grammage to length [m]
// Environment::GetDistance(step, next_step);
+ const double distance_interact = next_interact;
// ....
-
- // update step with actual min(min_step, next_step)
+
+ // convert next_decay from time to length [m]
+ // Environment::GetDistance(step, next_decay);
+ const double distance_decay = next_decay;
// ....
-
+
+ // take minimum of geometry, interaction, decay for next step
+ const double distance_decay_interact = std::min(next_decay, next_interact);
+ const double distance_next = std::min(distance_decay_interact, distance_max);
+
/// 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));
-
+ // .... also update time, momentum, direction, ...
+
+ // apply all continuous processes on particle + track
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
} else {
- // 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);
+ std::cout << "select process " << (distance_decay_interact < distance_max)
+ << std::endl;
+ // check if geometry limits step, then quit this step
+ if (distance_decay_interact < distance_max) {
+ // check weather decay or interaction limits this step
+ if (distance_decay > distance_interact) {
+ std::cout << "collide" << std::endl;
+ const double actual_inv_length =
+ fProcesseList.GetTotalInverseInteractionLength(particle, step);
+ const double sample_process = rmng() / (double)rmng.max();
+ double inv_lambda_count = 0;
+ fProcesseList.SelectInteraction(particle, fStack, actual_inv_length,
+ sample_process, inv_lambda_count);
+ } else {
+ std::cout << "decay" << std::endl;
+ const double actual_decay_time =
+ fProcesseList.GetTotalInverseLifetime(particle);
+ const double sample_process = rmng() / (double)rmng.max();
+ double inv_decay_count = 0;
+ fProcesseList.SelectDecay(particle, fStack, actual_decay_time, sample_process,
+ inv_decay_count);
+ }
+ }
}
}
diff --git a/Framework/Cascade/testCascade.cc b/Framework/Cascade/testCascade.cc
index f315cc292..6502dff5c 100644
--- a/Framework/Cascade/testCascade.cc
+++ b/Framework/Cascade/testCascade.cc
@@ -41,23 +41,17 @@ using namespace std;
static int fCount = 0;
-class ProcessSplit : public corsika::process::BaseProcess<ProcessSplit> {
+class ProcessSplit : public corsika::process::ContinuousProcess<ProcessSplit> {
public:
ProcessSplit() {}
-
- 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;
+ double MaxStepLength(Particle&, T&) const {
+ return 1;
}
-
- template <typename Particle, typename Stack>
- void DoDiscrete(Particle& p, Stack& s) const {
+
+ template <typename Particle, typename T, typename Stack>
+ void DoContinuous(Particle& p, T&, Stack& s) const {
EnergyType E = p.GetEnergy();
if (E < 85_MeV) {
p.Delete();
@@ -89,15 +83,15 @@ TEST_CASE("Cascade", "[Cascade]") {
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;
@@ -109,7 +103,7 @@ TEST_CASE("Cascade", "[Cascade]") {
particle.SetTime(0_ns);
EAS.Init();
EAS.Run();
-
+
/*
SECTION("sectionTwo") {
for (int i = 0; i < 0; ++i) {
@@ -119,7 +113,7 @@ 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 b17305d5a..9cc29196b 100644
--- a/Framework/ProcessSequence/BaseProcess.h
+++ b/Framework/ProcessSequence/BaseProcess.h
@@ -29,27 +29,29 @@ 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
+ typename BaseProcess< T<A, B> >::is_process_sequence
{
static const bool value = true;
};
*/
-
+
/*
template <typename T>
struct is_base {
diff --git a/Framework/ProcessSequence/CMakeLists.txt b/Framework/ProcessSequence/CMakeLists.txt
index c1f5a4586..6c3e5ff8c 100644
--- a/Framework/ProcessSequence/CMakeLists.txt
+++ b/Framework/ProcessSequence/CMakeLists.txt
@@ -12,7 +12,8 @@ set (
CORSIKAprocesssequence_HEADERS
BaseProcess.h
ContinuousProcess.h
- DiscreteProcess.h
+ InteractionProcess.h
+ DecayProcess.h
ProcessSequence.h
ProcessReturn.h
)
diff --git a/Framework/ProcessSequence/DecayProcess.h b/Framework/ProcessSequence/DecayProcess.h
new file mode 100644
index 000000000..8a49803ed
--- /dev/null
+++ b/Framework/ProcessSequence/DecayProcess.h
@@ -0,0 +1,51 @@
+
+/**
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#ifndef _include_corsika_decayprocess_h_
+#define _include_corsika_decayprocess_h_
+
+#include <corsika/process/ProcessReturn.h> // for convenience
+#include <corsika/setup/SetupTrajectory.h>
+
+namespace corsika::process {
+
+ /**
+ \class DecayProcess
+
+ The structural base type of a process object in a
+ ProcessSequence. Both, the ProcessSequence and all its elements
+ are of type DecayProcess<T>
+
+ */
+
+ template <typename derived>
+ struct DecayProcess {
+
+ derived& GetRef() { return static_cast<derived&>(*this); }
+ const derived& GetRef() const { return static_cast<const derived&>(*this); }
+
+ /// here starts the interface-definition part
+ // -> enforce derived to implement DoDecay...
+ template <typename Particle, typename Stack>
+ inline EProcessReturn DoDecay(Particle&, Stack&) const;
+
+ template <typename Particle>
+ inline double GetLifetime(Particle& p) const;
+
+ template <typename Particle>
+ inline double GetInverseLifetime(Particle& p) const {
+ return 1. / GetRef().GetLifetime(p);
+ }
+ };
+
+} // namespace corsika::process
+
+#endif
diff --git a/Framework/ProcessSequence/DiscreteProcess.h b/Framework/ProcessSequence/DiscreteProcess.h
index d0464e4eb..7ce50810f 100644
--- a/Framework/ProcessSequence/DiscreteProcess.h
+++ b/Framework/ProcessSequence/DiscreteProcess.h
@@ -36,13 +36,15 @@ namespace corsika::process {
/// here starts the interface-definition part
// -> enforce derived to implement DoDiscrete...
template <typename Particle, typename Stack>
- inline EProcessReturn DoDiscrete(Particle&, Stack&) const; // {}
+ inline EProcessReturn DoDiscrete(Particle&, Stack&) const;
+
+ template <typename Particle, typename Track>
+ inline double GetInteractionLength(Particle& p, Track& t) const;
template <typename Particle, typename Track>
inline double GetInverseInteractionLength(Particle& p, Track& t) const {
- return 1./GetRef().GetInteractionLength(p, t);
+ return 1. / GetRef().GetInteractionLength(p, t);
}
-
};
} // namespace corsika::process
diff --git a/Framework/ProcessSequence/InteractionProcess.h b/Framework/ProcessSequence/InteractionProcess.h
new file mode 100644
index 000000000..40264e606
--- /dev/null
+++ b/Framework/ProcessSequence/InteractionProcess.h
@@ -0,0 +1,51 @@
+
+/**
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#ifndef _include_corsika_interactionprocess_h_
+#define _include_corsika_interactionprocess_h_
+
+#include <corsika/process/ProcessReturn.h> // for convenience
+#include <corsika/setup/SetupTrajectory.h>
+
+namespace corsika::process {
+
+ /**
+ \class InteractionProcess
+
+ The structural base type of a process object in a
+ ProcessSequence. Both, the ProcessSequence and all its elements
+ are of type InteractionProcess<T>
+
+ */
+
+ template <typename derived>
+ struct InteractionProcess {
+
+ derived& GetRef() { return static_cast<derived&>(*this); }
+ const derived& GetRef() const { return static_cast<const derived&>(*this); }
+
+ /// here starts the interface-definition part
+ // -> enforce derived to implement DoInteraction...
+ template <typename Particle, typename Stack>
+ inline EProcessReturn DoInteraction(Particle&, Stack&) const;
+
+ template <typename Particle, typename Track>
+ inline double GetInteractionLength(Particle& p, Track& t) const;
+
+ template <typename Particle, typename Track>
+ inline double GetInverseInteractionLength(Particle& p, Track& t) const {
+ return 1. / GetRef().GetInteractionLength(p, t);
+ }
+ };
+
+} // namespace corsika::process
+
+#endif
diff --git a/Framework/ProcessSequence/ProcessReturn.h b/Framework/ProcessSequence/ProcessReturn.h
index 039a56222..afd75d690 100644
--- a/Framework/ProcessSequence/ProcessReturn.h
+++ b/Framework/ProcessSequence/ProcessReturn.h
@@ -24,7 +24,8 @@ namespace corsika::process {
eOk = 1,
eParticleAbsorbed = 2,
eInteracted = 3,
- };
+ eDecayed = 4,
+ };
} // namespace corsika::process
#endif
diff --git a/Framework/ProcessSequence/ProcessSequence.h b/Framework/ProcessSequence/ProcessSequence.h
index a24356d8e..8517d4345 100644
--- a/Framework/ProcessSequence/ProcessSequence.h
+++ b/Framework/ProcessSequence/ProcessSequence.h
@@ -14,11 +14,13 @@
#include <corsika/process/BaseProcess.h>
#include <corsika/process/ContinuousProcess.h>
-#include <corsika/process/DiscreteProcess.h>
+#include <corsika/process/DecayProcess.h>
+//#include <corsika/process/DiscreteProcess.h>
+#include <corsika/process/InteractionProcess.h>
#include <corsika/process/ProcessReturn.h>
-#include <limits>
#include <cmath>
+#include <limits>
//#include <corsika/setup/SetupTrajectory.h>
// using corsika::setup::Trajectory;
@@ -137,24 +139,15 @@ 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 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;
@@ -169,27 +162,29 @@ namespace corsika::process {
template <typename Particle, typename Track, typename Stack>
inline EProcessReturn DoContinuous(Particle& p, Track& t, Stack& s) const {
EProcessReturn ret = EProcessReturn::eOk;
- if constexpr (!std::is_base_of<DiscreteProcess<T1>, T1>::value) {
- // A.DoContinuous(std::forward<Particle>(p), t, std::forward<Stack>(s));
+ if constexpr (std::is_base_of<ContinuousProcess<T1>, T1>::value ||
+ is_process_sequence<T1>::value) {
A.DoContinuous(p, t, s);
}
- if constexpr (!std::is_base_of<DiscreteProcess<T2>, T2>::value) {
- // B.DoContinuous(std::forward<Particle>(p), t, std::forward<Stack>(s));
+ if constexpr (std::is_base_of<ContinuousProcess<T2>, T2>::value ||
+ is_process_sequence<T2>::value) {
B.DoContinuous(p, t, s);
}
return ret;
}
template <typename Particle, typename 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;
+ inline double MaxStepLength(Particle& p, Track& track) const {
+ double max_length = std::numeric_limits<double>::infinity();
+ if constexpr (std::is_base_of<ContinuousProcess<T1>, T1>::value ||
+ is_process_sequence<T1>::value) {
+ max_length = std::min(max_length, A.MaxStepLength(p, track));
+ }
+ if constexpr (std::is_base_of<ContinuousProcess<T2>, T2>::value ||
+ is_process_sequence<T2>::value) {
+ max_length = std::min(max_length, B.MaxStepLength(p, track));
+ }
+ return max_length;
}
template <typename Particle, typename Track>
@@ -205,64 +200,116 @@ namespace corsika::process {
template <typename Particle, typename Track>
inline double GetInverseInteractionLength(Particle& p, Track& t) const {
double tot = 0;
- if constexpr (!std::is_base_of<ContinuousProcess<T1>, T1>::value) {
+ if constexpr (std::is_base_of<InteractionProcess<T1>, T1>::value ||
+ is_process_sequence<T1>::value) {
tot += A.GetInverseInteractionLength(p, t);
}
- if constexpr (!std::is_base_of<ContinuousProcess<T2>, T2>::value) {
+ if constexpr (std::is_base_of<InteractionProcess<T2>, T2>::value ||
+ is_process_sequence<T2>::value) {
tot += B.GetInverseInteractionLength(p, t);
}
return tot;
}
template <typename Particle, typename Stack>
- inline EProcessReturn DoDiscrete(Particle& p, Stack& s) const {
- if constexpr (!std::is_base_of<ContinuousProcess<T1>, T1>::value) {
- A.DoDiscrete(p, s);
+ inline EProcessReturn SelectInteraction(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.SelectInteraction(p, s, lambda_inv_count, rndm_select, lambda_inv_count);
+ // if A did succeed, stop routine
+ if (ret != EProcessReturn::eOk) { return ret; }
+ } else if constexpr (std::is_base_of<InteractionProcess<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.DoInteraction(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.SelectInteraction(p, s, lambda_inv_count, rndm_select, lambda_inv_count);
+ // if A did succeed, stop routine
+ if (ret != EProcessReturn::eOk) { return ret; }
+ } else if constexpr (std::is_base_of<InteractionProcess<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.DoInteraction(p, s);
+ return EProcessReturn::eInteracted;
+ }
+ } // end branch A
+ return EProcessReturn::eOk;
+ }
+
+ template <typename Particle>
+ inline double GetTotalLifetime(Particle& p) const {
+ return 1. / GetInverseLifetime(p);
+ }
+
+ template <typename Particle>
+ inline double GetTotalInverseLifetime(Particle& p) const {
+ return GetInverseLifetime(p);
+ }
+
+ template <typename Particle>
+ inline double GetInverseLifetime(Particle& p) const {
+ double tot = 0;
+ if constexpr (std::is_base_of<DecayProcess<T1>, T1>::value ||
+ is_process_sequence<T1>::value) {
+ tot += A.GetInverseLifetime(p);
}
- if constexpr (!std::is_base_of<ContinuousProcess<T2>, T2>::value) {
- B.DoDiscrete(p, s);
+ if constexpr (std::is_base_of<DecayProcess<T2>, T2>::value ||
+ is_process_sequence<T2>::value) {
+ tot += B.GetInverseLifetime(p);
}
- return EProcessReturn::eOk;
+ return tot;
}
+ // select decay process
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 {
+ inline EProcessReturn SelectDecay(Particle& p, Stack& s, const double decay_inv_tot,
+ const double rndm_select,
+ double& decay_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 A is a process sequence --> check inside
+ const EProcessReturn ret =
+ A.SelectDecay(p, s, decay_inv_count, rndm_select, decay_inv_count);
+ // if A did succeed, stop routine
+ if (ret != EProcessReturn::eOk) { return ret; }
+ } else if constexpr (std::is_base_of<DecayProcess<T1>, T1>::value) {
+ // if this is not a ContinuousProcess --> evaluate probability
+ decay_inv_count += A.GetInverseLifetime(p);
+ // check if we should execute THIS process and then EXIT
+ if (rndm_select < decay_inv_count / decay_inv_tot) {
+ A.DoDecay(p, s);
+ return EProcessReturn::eDecayed;
+ }
+ } // 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
+ // if A is a process sequence --> check inside
+ const EProcessReturn ret =
+ B.SelectDecay(p, s, decay_inv_count, rndm_select, decay_inv_count);
+ // if A did succeed, stop routine
+ if (ret != EProcessReturn::eOk) { return ret; }
+ } else if constexpr (std::is_base_of<DecayProcess<T2>, T2>::value) {
+ // if this is not a ContinuousProcess --> evaluate probability
+ decay_inv_count += B.GetInverseLifetime(p);
+ // check if we should execute THIS process and then EXIT
+ if (rndm_select < decay_inv_count / decay_inv_tot) {
+ B.DoDecay(p, s);
+ return EProcessReturn::eDecayed;
+ }
+ } // end branch B
return EProcessReturn::eOk;
}
@@ -273,8 +320,8 @@ namespace corsika::process {
}
};
- /// the +operator assembles many BaseProcess, ContinuousProcess, and
- /// DiscreteProcess objects into a ProcessSequence, all combinatorics
+ /// the + operator assembles many BaseProcess, ContinuousProcess, and
+ /// InteractionProcess objects into a ProcessSequence, all combinatorics
/// must be allowed, this is why we define a macro to define all
/// combinations here:
@@ -285,93 +332,26 @@ namespace corsika::process {
}
OPSEQ(BaseProcess, BaseProcess)
- OPSEQ(BaseProcess, DiscreteProcess)
+ OPSEQ(BaseProcess, InteractionProcess)
OPSEQ(BaseProcess, ContinuousProcess)
+ OPSEQ(BaseProcess, DecayProcess)
OPSEQ(ContinuousProcess, BaseProcess)
- OPSEQ(ContinuousProcess, DiscreteProcess)
+ OPSEQ(ContinuousProcess, InteractionProcess)
OPSEQ(ContinuousProcess, ContinuousProcess)
- OPSEQ(DiscreteProcess, BaseProcess)
- OPSEQ(DiscreteProcess, DiscreteProcess)
- OPSEQ(DiscreteProcess, ContinuousProcess)
-
- /*
- template <typename T1>
- struct depth_lhs
- {
- static const int num = 0;
- };
-
-
-
- // terminating condition
- template <typename T1, typename T2>
- struct depth_lhs< Sequence<T1,T2> >
- {
- // try to expand the left node (T1) which might be a Sequence type
- static const int num = 1 + depth_lhs<T1>::num;
- };
- */
-
- /*
- template <typename T1>
- struct mat_ptrs
- {
- static const int num = 0;
-
- inline static void
- get_ptrs(const Process** ptrs, const T1& X)
- {
- ptrs[0] = reinterpret_cast<const Process*>(&X);
- }
- };
-
-
- template <typename T1, typename T2>
- struct mat_ptrs< Sequence<T1,T2> >
- {
- static const int num = 1 + mat_ptrs<T1>::num;
-
- inline static void
- get_ptrs(const Process** in_ptrs, const Sequence<T1,T2>& X)
- {
- // traverse the left node
- mat_ptrs<T1>::get_ptrs(in_ptrs, X.A);
- // get address of the matrix on the right node
- in_ptrs[num] = reinterpret_cast<const Process*>(&X.B);
- }
- };
- */
-
- /*
- template<typename T1, typename T2>
- const Process&
- Process::operator=(const Sequence<T1,T2>& X)
- {
- int N = 1 + depth_lhs< Sequence<T1,T2> >::num;
- const Process* ptrs[N];
- mat_ptrs< Sequence<T1,T2> >::get_ptrs(ptrs, X);
- int r = ptrs[0]->rows;
- int c = ptrs[0]->cols;
- // ... check that all matrices have the same size ...
- set_size(r, c);
- for(int j=0; j<r*c; ++j)
- {
- double sum = ptrs[0]->data[j];
- for(int i=1; i<N; ++i)
- {
- sum += ptrs[i]->data[j];
- }
- data[j] = sum;
- }
- return *this;
- }
- */
-
- template<template<typename, typename> class T, typename A, typename B>
- struct is_process_sequence< T<A,B> >
- {
- static const bool value = true;
- };
+ OPSEQ(ContinuousProcess, DecayProcess)
+ OPSEQ(InteractionProcess, BaseProcess)
+ OPSEQ(InteractionProcess, InteractionProcess)
+ OPSEQ(InteractionProcess, ContinuousProcess)
+ OPSEQ(InteractionProcess, DecayProcess)
+ OPSEQ(DecayProcess, BaseProcess)
+ OPSEQ(DecayProcess, InteractionProcess)
+ OPSEQ(DecayProcess, ContinuousProcess)
+ OPSEQ(DecayProcess, DecayProcess)
+
+ 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
diff --git a/Framework/ProcessSequence/testProcessSequence.cc b/Framework/ProcessSequence/testProcessSequence.cc
index 008201b7e..af9bfeead 100644
--- a/Framework/ProcessSequence/testProcessSequence.cc
+++ b/Framework/ProcessSequence/testProcessSequence.cc
@@ -66,7 +66,7 @@ public:
}
};
-class Process1 : public DiscreteProcess<Process1> {
+class Process1 : public InteractionProcess<Process1> {
public:
Process1(const int v)
: fV(v) {}
@@ -76,7 +76,7 @@ public:
globalCount++;
}
template <typename D, typename S>
- inline EProcessReturn DoDiscrete(D& d, S&) const {
+ inline EProcessReturn DoInteraction(D& d, S&) const {
for (int i = 0; i < nData; ++i) d.p[i] += 1 + i;
return EProcessReturn::eOk;
}
@@ -84,7 +84,7 @@ public:
int fV;
};
-class Process2 : public DiscreteProcess<Process2> {
+class Process2 : public InteractionProcess<Process2> {
int fV = 0;
public:
@@ -96,8 +96,8 @@ public:
globalCount++;
}
template <typename Particle, typename Stack>
- inline EProcessReturn DoDiscrete(Particle&, Stack&) const {
- cout << "Process2::DoDiscrete" << endl;
+ inline EProcessReturn DoInteraction(Particle&, Stack&) const {
+ cout << "Process2::DoInteraction" << endl;
return EProcessReturn::eOk;
}
template <typename Particle, typename Track>
@@ -107,7 +107,7 @@ public:
}
};
-class Process3 : public DiscreteProcess<Process3> {
+class Process3 : public InteractionProcess<Process3> {
int fV = 0;
public:
@@ -119,8 +119,8 @@ public:
globalCount++;
}
template <typename Particle, typename Stack>
- inline EProcessReturn DoDiscrete(Particle&, Stack&) const {
- cout << "Process3::DoDiscrete" << endl;
+ inline EProcessReturn DoInteraction(Particle&, Stack&) const {
+ cout << "Process3::DoInteraction" << endl;
return EProcessReturn::eOk;
}
template <typename Particle, typename Track>
@@ -148,7 +148,28 @@ public:
}
// inline double MinStepLength(D& d) {
template <typename Particle, typename Stack>
- EProcessReturn DoDiscrete(Particle&, Stack&) const {
+ EProcessReturn DoInteraction(Particle&, Stack&) const {
+ return EProcessReturn::eOk;
+ }
+};
+
+class Decay1 : public DecayProcess<Decay1> {
+ int fV = 0;
+
+public:
+ Decay1(const int v)
+ : fV(v) {}
+ void Init() {
+ cout << "Decay1::Init" << endl;
+ assert(globalCount == fV);
+ globalCount++;
+ }
+ template <typename Particle>
+ double GetLifetime(Particle&) const {
+ return 1;
+ }
+ template <typename Particle, typename Stack>
+ EProcessReturn DoDecay(Particle&, Stack&) const {
return EProcessReturn::eOk;
}
};
@@ -192,7 +213,21 @@ TEST_CASE("Process Sequence", "[Process Sequence]") {
double tot_inv = sequence2.GetTotalInverseInteractionLength(s, t);
cout << "lambda_tot=" << tot << " lambda_tot_inv=" << tot_inv << endl;
}
-
+
+ SECTION("lifetime") {
+ ContinuousProcess1 cp1(0);
+ Process2 m2(1);
+ Process3 m3(2);
+ Decay1 d3(2);
+
+ DummyStack s;
+
+ const auto sequence2 = cp1 + m2 + m3 + d3;
+ double tot = sequence2.GetTotalLifetime(s);
+ double tot_inv = sequence2.GetTotalInverseLifetime(s);
+ cout << "lambda_tot=" << tot << " lambda_tot_inv=" << tot_inv << endl;
+ }
+
SECTION("sectionTwo") {
ContinuousProcess1 cp1(0);
@@ -213,14 +248,14 @@ TEST_CASE("Process Sequence", "[Process Sequence]") {
sequence2.DoContinuous(p, t, s);
cout << "-->dodisc" << endl;
- sequence2.DoDiscrete(p, s);
+ // sequence2.DoInteraction(p, s);
cout << "-->done" << endl;
const int nLoop = 5;
cout << "Running loop with n=" << nLoop << endl;
for (int i = 0; i < nLoop; ++i) {
sequence2.DoContinuous(p, t, s);
- sequence2.DoDiscrete(p, s);
+ // sequence2.DoInteraction(p, s);
}
for (int i = 0; i < nData; i++) { cout << "data[" << i << "]=" << p.p[i] << endl; }
cout << "done" << endl;
diff --git a/Framework/StackInterface/Stack.h b/Framework/StackInterface/Stack.h
index e71880ccf..f47c15d3d 100644
--- a/Framework/StackInterface/Stack.h
+++ b/Framework/StackInterface/Stack.h
@@ -67,9 +67,7 @@ namespace corsika::stack {
IncrementSize();
return StackIterator(*this, GetSize() - 1);
}
- void Copy(StackIterator& a, StackIterator& b) {
- Copy(a.GetIndex(), b.GetIndex());
- }
+ 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 abf47afd1..1076803b8 100644
--- a/Processes/StackInspector/StackInspector.cc
+++ b/Processes/StackInspector/StackInspector.cc
@@ -17,8 +17,8 @@
#include <corsika/setup/SetupTrajectory.h>
-#include <limits>
#include <iostream>
+#include <limits>
using namespace std;
using namespace corsika;
@@ -57,7 +57,7 @@ process::EProcessReturn StackInspector<Stack>::DoContinuous(Particle&, setup::Tr
}
template <typename Stack>
-double StackInspector<Stack>::MinStepLength(Particle&, setup::Trajectory&) const {
+double StackInspector<Stack>::MaxStepLength(Particle&, setup::Trajectory&) const {
return std::numeric_limits<double>::infinity();
}
diff --git a/Processes/StackInspector/StackInspector.h b/Processes/StackInspector/StackInspector.h
index 38d5cfc01..58b9ec5eb 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>
- double MinStepLength(Particle&, corsika::setup::Trajectory&) const;
+ double MaxStepLength(Particle&, corsika::setup::Trajectory&) const;
private:
bool fReport;
--
GitLab