diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index 09ae4815c8fb763c48f5871cc9eafab52777c8a2..d7e673ee574802daff457fdd00a4962bc1338458 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 a1e6ff8a27f162b6a6cdd4fcc1940e01b3188068..7cd9d077de0bc6f84174789771714978d8ac3f94 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 f315cc292d66f043526befa3d4276592edfc4187..6502dff5c2edaefd0c0df75b9a1b42abd284ff50 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 b17305d5ac6b89a53eba83454db77b1e05d980f9..9cc29196be0d136598e2f32480acd35dd16b158e 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 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 d0464e4ebc3e8f3b72930e9b7c562ff49fcf4b2a..7ce50810f7fccbbbcb0c707d4124ec1ad03f0244 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 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 039a562221ebd3eb0b3e235857476f4359ab3d3b..afd75d690a3b5c1ee54d7e0884a536f42221d2a2 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 a24356d8eac1cc72b9f55a21cde06c7141cd7f32..8517d4345202f8fb667f830eab3c093937b2ad05 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 008201b7e5dd50826e3df3f8cb850c6f35087bcf..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,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 e71880ccf2cff35f16abd517c03bfcaede8aa01d..f47c15d3daf3071a8c71da1b81b71e575a620784 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 abf47afd1dcd3d0815375e2c771864b0f9488987..1076803b88205f54570c447bc1c232e07fa7ccfc 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 38d5cfc016fb2eda8e1eba9ce30ab0e9dad7adf3..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>
- double MinStepLength(Particle&, corsika::setup::Trajectory&) const;
+ double MaxStepLength(Particle&, corsika::setup::Trajectory&) const;
private:
bool fReport;