diff --git a/CMakeLists.txt b/CMakeLists.txt
index b22b7c9c0d18e7a389974c9d4e9496188d2a1d21..62d488fa2713a44964460d939aa6a5742b1840d8 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)
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index c3a8c95191effe7928ac470281617c4f34d486dd..49ab47047d476934a45de114f015cbfeea2c260b 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -46,12 +46,12 @@ using namespace corsika::units::si;
static int fCount = 0;
-class ProcessSplit : public corsika::process::BaseProcess<ProcessSplit> {
+class ProcessSplit : public corsika::process::InteractionProcess<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;
@@ -90,15 +90,17 @@ 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>
@@ -107,8 +109,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())) {
@@ -133,11 +135,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 3fac21b57739e5a343baf9e826e7abd67b79984d..7cd9d077de0bc6f84174789771714978d8ac3f94 100644
--- a/Framework/Cascade/Cascade.h
+++ b/Framework/Cascade/Cascade.h
@@ -13,11 +13,14 @@
#define _include_corsika_cascade_Cascade_h_
#include <corsika/process/ProcessReturn.h>
+#include <corsika/random/RNGManager.h>
+#include <corsika/setup/SetupTrajectory.h>
#include <corsika/units/PhysicalUnits.h>
#include <type_traits>
-#include <corsika/setup/SetupTrajectory.h>
+using namespace corsika;
+using namespace corsika::units::si;
namespace corsika::cascade {
@@ -48,24 +51,92 @@ namespace corsika::cascade {
}
// do cascade equations, which can put new particles on Stack,
// thus, the double loop
- // DoCascadeEquations(); //
+ // DoCascadeEquations();
}
}
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);
- fProcesseList.MinStepLength(particle, step);
+
+ // determine combined total interaction length (inverse)
+ const double total_inv_lambda =
+ fProcesseList.GetTotalInverseInteractionLength(particle, step);
+
+ // sample random exponential step length
+ const double sample_step = rmng() / (double)rmng.max();
+ 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);
+
+ // 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;
+ // ....
+
+ // 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
+ // fStack.Delete(particle); // TODO: check if this is really needed
} else {
- fProcesseList.DoDiscrete(particle, fStack);
+
+ 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 6bb0a10ddd39956df11e469739faf4dce3cfef7a..423940c156782b2978ed402524f209cf489178fb 100644
--- a/Framework/Cascade/testCascade.cc
+++ b/Framework/Cascade/testCascade.cc
@@ -20,6 +20,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>
@@ -43,20 +45,17 @@ using namespace corsika::units::si;
static int fCount = 0;
-class ProcessSplit : public corsika::process::BaseProcess<ProcessSplit> {
+class ProcessSplit : public corsika::process::ContinuousProcess<ProcessSplit> {
public:
ProcessSplit() {}
- template <typename Particle>
- void MinStepLength(Particle&, setup::Trajectory&) const {}
-
- template <typename Particle, typename Stack>
- EProcessReturn DoContinuous(Particle&, setup::Trajectory&, Stack&) const {
- return EProcessReturn::eOk;
+ template <typename Particle, typename T>
+ 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();
@@ -64,7 +63,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());
@@ -79,6 +80,12 @@ private:
};
TEST_CASE("Cascade", "[Cascade]") {
+
+ corsika::random::RNGManager& rmng = corsika::random::RNGManager::GetInstance();
+ ;
+ const std::string str_name = "s_rndm";
+ rmng.RegisterRandomStream(str_name);
+
corsika::environment::Environment env; // dummy environment
auto& universe = *(env.GetUniverse());
auto const radius = 1_m * std::numeric_limits<double>::infinity();
@@ -95,19 +102,21 @@ TEST_CASE("Cascade", "[Cascade]") {
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();
@@ -120,4 +129,5 @@ TEST_CASE("Cascade", "[Cascade]") {
// cout << "Result: E0=" << E0 / 1_GeV << "GeV, count=" << p1.GetCount() << endl;
}
}
+ */
}
diff --git a/Framework/ProcessSequence/BaseProcess.h b/Framework/ProcessSequence/BaseProcess.h
index 34dcb198e4d26f409837d253f3270cc0f488098b..9cc29196be0d136598e2f32480acd35dd16b158e 100644
--- a/Framework/ProcessSequence/BaseProcess.h
+++ b/Framework/ProcessSequence/BaseProcess.h
@@ -29,14 +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
+ {
+ static const bool value = true;
+ };
+ */
+
/*
template <typename T>
struct is_base {
diff --git a/Framework/ProcessSequence/CMakeLists.txt b/Framework/ProcessSequence/CMakeLists.txt
index c1f5a45862faffc92aba606e8fc1abde94220d89..6c3e5ff8cc26e571955b937920dc086fd925869c 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 0000000000000000000000000000000000000000..8a49803ed49913a1c8cf5c217a5c866de307f4bd
--- /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 a540c92c783dcad931f0d2bf8278c057736bda79..7ce50810f7fccbbbcb0c707d4124ec1ad03f0244 100644
--- a/Framework/ProcessSequence/DiscreteProcess.h
+++ b/Framework/ProcessSequence/DiscreteProcess.h
@@ -36,7 +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);
+ }
};
} // namespace corsika::process
diff --git a/Framework/ProcessSequence/InteractionProcess.h b/Framework/ProcessSequence/InteractionProcess.h
new file mode 100644
index 0000000000000000000000000000000000000000..40264e606a7de8ecf607fa8d5173190e9e4ca0e4
--- /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 82cc816c5ae9bd9d731f00c33f9bf0321a02d650..afd75d690a3b5c1ee54d7e0884a536f42221d2a2 100644
--- a/Framework/ProcessSequence/ProcessReturn.h
+++ b/Framework/ProcessSequence/ProcessReturn.h
@@ -23,6 +23,8 @@ namespace corsika::process {
enum class EProcessReturn {
eOk = 1,
eParticleAbsorbed = 2,
+ eInteracted = 3,
+ eDecayed = 4,
};
} // namespace corsika::process
diff --git a/Framework/ProcessSequence/ProcessSequence.h b/Framework/ProcessSequence/ProcessSequence.h
index 54ac90e7a18f7c40c1e5a53c82ebf854e5abe486..8517d4345202f8fb667f830eab3c093937b2ad05 100644
--- a/Framework/ProcessSequence/ProcessSequence.h
+++ b/Framework/ProcessSequence/ProcessSequence.h
@@ -14,9 +14,14 @@
#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 <cmath>
+#include <limits>
+
//#include <corsika/setup/SetupTrajectory.h>
// using corsika::setup::Trajectory;
//#include <variant>
@@ -25,7 +30,7 @@
namespace corsika::process {
- /* namespace detail { */
+ // namespace detail {
/* /\* template<typename TT1, typename TT2, typename Type = void> *\/ */
/* /\* struct CallHello { *\/ */
@@ -41,7 +46,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 +126,7 @@ namespace corsika::process {
/* } */
/* }; */
/* *\/ */
- /* } // end namespace detail */
+ //} // end namespace detail
/**
\class ProcessSequence
@@ -134,8 +139,15 @@ namespace corsika::process {
https://en.wikipedia.org/wiki/Curiously_recurring_template_pattern
*/
+ // 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> > {
+
public:
const T1& A;
const T2& B;
@@ -150,41 +162,154 @@ 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 void MinStepLength(Particle& p, Track& track) const {
- A.MinStepLength(p, track);
- B.MinStepLength(p, track);
+ 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>
- 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 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 double GetInverseInteractionLength(Particle& p, Track& t) const {
+ double tot = 0;
+ 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<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 tot;
+ }
+
+ // select decay process
+ template <typename Particle, typename Stack>
+ 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.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.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;
}
@@ -195,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:
@@ -207,87 +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;
- }
- */
+ 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 7094cacb94e8da03bcfce1b9190160fb1b8c1898..af9bfeeade6c2d149f88b803f0ddd2ebb105e6aa 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,13 +96,18 @@ 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>
+ inline double GetInteractionLength(Particle&, Track&) const {
+ cout << "Process2::GetInteractionLength" << endl;
+ return 3;
+ }
};
-class Process3 : public DiscreteProcess<Process3> {
+class Process3 : public InteractionProcess<Process3> {
int fV = 0;
public:
@@ -114,10 +119,15 @@ 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>
+ inline double GetInteractionLength(Particle&, Track&) const {
+ cout << "Process3::GetInteractionLength" << endl;
+ return 1.;
+ }
};
class Process4 : public BaseProcess<Process4> {
@@ -138,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;
}
};
@@ -169,6 +200,34 @@ 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("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);
@@ -189,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 b4660144df60754df6ff8e78e45d0920f5a58cf8..f47c15d3daf3071a8c71da1b81b71e575a620784 100644
--- a/Framework/StackInterface/Stack.h
+++ b/Framework/StackInterface/Stack.h
@@ -67,6 +67,7 @@ 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 1986c83bb4471b1adee8b2b13653cb64c7ef5167..266005cde54295a790a8748e91614bf6f9d58b37 100644
--- a/Processes/StackInspector/StackInspector.cc
+++ b/Processes/StackInspector/StackInspector.cc
@@ -19,6 +19,7 @@
#include <corsika/setup/SetupTrajectory.h>
#include <iostream>
+#include <limits>
using namespace std;
using namespace corsika;
@@ -57,8 +58,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>::MaxStepLength(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 7b68d92d3fa6f5527de8bc42b76059017a79221f..58b9ec5ebd7652b127f719beafcedacbbf79d503 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 MaxStepLength(Particle&, corsika::setup::Trajectory&) const;
private:
bool fReport;
diff --git a/Stack/SuperStupidStack/SuperStupidStack.h b/Stack/SuperStupidStack/SuperStupidStack.h
index d1c8bc2b38622698b359afe8805d22403d99a4ed..8c79934422420791c1d33395810e5b90afa4038e 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]);
}