From a8cf60a9f138030348917c1e52a217e5aaae4cb2 Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Mon, 18 Jan 2021 17:25:48 +0100
Subject: [PATCH] almost works, only switch-pro-seq not
---
corsika/detail/framework/core/Cascade.inl | 10 +-
corsika/detail/framework/geometry/Plane.inl | 8 +
corsika/detail/framework/geometry/Point.inl | 6 +
.../framework/geometry/QuantityVector.inl | 2 +-
corsika/detail/framework/geometry/Vector.inl | 8 +
.../framework/process/ProcessSequence.inl | 144 ++++---
.../process/SwitchProcessSequence.inl | 125 +++---
.../detail/modules/LongitudinalProfile.inl | 2 +-
corsika/detail/modules/ObservationPlane.inl | 73 ++--
corsika/detail/modules/ParticleCut.inl | 2 +-
corsika/detail/modules/TrackWriter.inl | 3 +-
.../modules/energy_loss/BetheBlochPDG.inl | 4 +-
.../modules/proposal/ContinuousProcess.inl | 2 +-
.../modules/tracking/TrackingStraight.inl | 2 +
corsika/framework/geometry/Plane.hpp | 4 +
corsika/framework/geometry/QuantityVector.hpp | 4 +-
corsika/framework/process/BaseProcess.hpp | 4 +
.../framework/process/ContinuousProcess.hpp | 41 +-
.../process/ContinuousProcessIndex.hpp | 26 ++
.../process/ContinuousProcessStepLength.hpp | 43 ++
corsika/framework/process/ProcessSequence.hpp | 99 +++--
corsika/framework/process/ProcessTraits.hpp | 63 +--
.../process/SwitchProcessSequence.hpp | 69 +++-
corsika/modules/LongitudinalProfile.hpp | 4 +-
corsika/modules/ParticleCut.hpp | 6 +-
corsika/modules/TrackWriter.hpp | 2 +-
corsika/modules/energy_loss/BetheBlochPDG.hpp | 12 +-
.../modules/proposal/ContinuousProcess.hpp | 5 +-
examples/staticsequence_example.cpp | 10 +-
examples/vertical_EAS.cpp | 6 +-
modules/conex | 2 +-
modules/data | 2 +-
tests/framework/testProcessSequence.cpp | 385 ++++++++++++++----
tests/modules/testObservationPlane.cpp | 56 ++-
tests/modules/testParticleCut.cpp | 1 +
35 files changed, 895 insertions(+), 340 deletions(-)
create mode 100644 corsika/framework/process/ContinuousProcessIndex.hpp
create mode 100644 corsika/framework/process/ContinuousProcessStepLength.hpp
diff --git a/corsika/detail/framework/core/Cascade.inl b/corsika/detail/framework/core/Cascade.inl
index 8cd5b2e39..0e9161259 100644
--- a/corsika/detail/framework/core/Cascade.inl
+++ b/corsika/detail/framework/core/Cascade.inl
@@ -10,6 +10,8 @@
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/process/ProcessReturn.hpp>
+#include <corsika/framework/process/ContinuousProcessStepLength.hpp>
+#include <corsika/framework/process/ContinuousProcessIndex.hpp>
#include <corsika/framework/random/ExponentialDistribution.hpp>
#include <corsika/framework/random/RNGManager.hpp>
#include <corsika/framework/random/UniformRealDistribution.hpp>
@@ -118,7 +120,10 @@ namespace corsika {
next_interact);
// determine the maximum geometric step length
- LengthType const continuous_max_dist = sequence_.getMaxStepLength(vParticle, step);
+ ContinuousProcessStepLength const continuousMaxStep =
+ sequence_.getMaxStepLength(vParticle, step);
+ LengthType const continuous_max_dist = continuousMaxStep;
+ ContinuousProcessIndex const limitingId = continuousMaxStep;
// take minimum of geometry, interaction, decay for next step
auto const min_distance =
@@ -142,7 +147,8 @@ namespace corsika {
vParticle.setTime(vParticle.getTime() + step.getDuration());
// apply all continuous processes on particle + track
- if (sequence_.doContinuous(vParticle, step) == ProcessReturn::ParticleAbsorbed) {
+ if (sequence_.doContinuous(vParticle, step, limitingId) ==
+ ProcessReturn::ParticleAbsorbed) {
CORSIKA_LOG_DEBUG("Cascade: delete absorbed particle PID={} E={} GeV",
vParticle.getPID(), vParticle.getEnergy() / 1_GeV);
if (!vParticle.isErased()) vParticle.erase();
diff --git a/corsika/detail/framework/geometry/Plane.inl b/corsika/detail/framework/geometry/Plane.inl
index e158e04d3..d32167cf1 100644
--- a/corsika/detail/framework/geometry/Plane.inl
+++ b/corsika/detail/framework/geometry/Plane.inl
@@ -12,6 +12,8 @@
#include <corsika/framework/geometry/Point.hpp>
#include <corsika/framework/geometry/Vector.hpp>
+#include <sstream>
+
namespace corsika {
inline bool Plane::isAbove(Point const& vP) const {
@@ -26,4 +28,10 @@ namespace corsika {
inline Plane::DimLessVec const& Plane::getNormal() const { return normal_; }
+ inline std::string Plane::asString() const {
+ std::ostringstream txt;
+ txt << "center=" << center_ << ", normal=" << normal_;
+ return txt.str();
+ }
+
} // namespace corsika
diff --git a/corsika/detail/framework/geometry/Point.inl b/corsika/detail/framework/geometry/Point.inl
index 6c550c1ce..e424812f0 100644
--- a/corsika/detail/framework/geometry/Point.inl
+++ b/corsika/detail/framework/geometry/Point.inl
@@ -95,4 +95,10 @@ namespace corsika {
return Vector<length_d>(cs, getCoordinates() - pB.getCoordinates(cs));
}
+ inline std::ostream& operator<<(std::ostream& os, corsika::Point const& p) {
+ auto const& qv = p.getCoordinates();
+ os << qv << " m";
+ return os;
+ }
+
} // namespace corsika
diff --git a/corsika/detail/framework/geometry/QuantityVector.inl b/corsika/detail/framework/geometry/QuantityVector.inl
index 9df35aab4..6ba02e4f8 100644
--- a/corsika/detail/framework/geometry/QuantityVector.inl
+++ b/corsika/detail/framework/geometry/QuantityVector.inl
@@ -145,7 +145,7 @@ namespace corsika {
template <typename TDimension>
inline std::ostream& operator<<(std::ostream& os,
- corsika::QuantityVector<TDimension> const qv) {
+ corsika::QuantityVector<TDimension> const& qv) {
using quantity_type = phys::units::quantity<TDimension, double>;
os << '(' << qv.eigenVector_(0) << ' ' << qv.eigenVector_(1) << ' '
diff --git a/corsika/detail/framework/geometry/Vector.inl b/corsika/detail/framework/geometry/Vector.inl
index cbfd597b3..2ecc82e65 100644
--- a/corsika/detail/framework/geometry/Vector.inl
+++ b/corsika/detail/framework/geometry/Vector.inl
@@ -228,4 +228,12 @@ namespace corsika {
bareResult);
}
+ template <typename TDimension>
+ inline std::ostream& operator<<(std::ostream& os,
+ corsika::Vector<TDimension> const& v) {
+ auto const& qv = v.getComponents();
+ os << qv;
+ return os;
+ }
+
} // namespace corsika
diff --git a/corsika/detail/framework/process/ProcessSequence.inl b/corsika/detail/framework/process/ProcessSequence.inl
index 2ca682b4a..19056ad6c 100644
--- a/corsika/detail/framework/process/ProcessSequence.inl
+++ b/corsika/detail/framework/process/ProcessSequence.inl
@@ -12,6 +12,8 @@
#include <corsika/framework/process/BaseProcess.hpp>
#include <corsika/framework/process/BoundaryCrossingProcess.hpp>
#include <corsika/framework/process/ContinuousProcess.hpp>
+#include <corsika/framework/process/ContinuousProcessStepLength.hpp>
+#include <corsika/framework/process/ContinuousProcessIndex.hpp>
#include <corsika/framework/process/DecayProcess.hpp>
#include <corsika/framework/process/InteractionProcess.hpp>
#include <corsika/framework/process/ProcessReturn.hpp>
@@ -24,11 +26,14 @@
namespace corsika {
- template <typename TProcess1, typename TProcess2, int IndexStart, int IndexProcess1, int IndexProcess2>
+ template <typename TProcess1, typename TProcess2, int IndexStart, int IndexProcess1,
+ int IndexProcess2>
template <typename TParticle>
- ProcessReturn ProcessSequence<TProcess1, TProcess2, IndexStart, IndexProcess1, IndexProcess2>::doBoundaryCrossing(
- TParticle& particle, typename TParticle::node_type const& from,
- typename TParticle::node_type const& to) {
+ ProcessReturn ProcessSequence<
+ TProcess1, TProcess2, IndexStart, IndexProcess1,
+ IndexProcess2>::doBoundaryCrossing(TParticle& particle,
+ typename TParticle::node_type const& from,
+ typename TParticle::node_type const& to) {
ProcessReturn ret = ProcessReturn::Ok;
if constexpr (std::is_base_of_v<BoundaryCrossingProcess<process1_type>,
@@ -46,25 +51,38 @@ namespace corsika {
return ret;
}
- template <typename TProcess1, typename TProcess2, int IndexStart, int IndexProcess1, int IndexProcess2>
+ template <typename TProcess1, typename TProcess2, int IndexStart, int IndexProcess1,
+ int IndexProcess2>
template <typename TParticle, typename TTrack>
- ProcessReturn ProcessSequence<TProcess1, TProcess2, IndexStart, IndexProcess1, IndexProcess2>::doContinuous(TParticle& particle,
- TTrack& vT) {
+ ProcessReturn
+ ProcessSequence<TProcess1, TProcess2, IndexStart, IndexProcess1, IndexProcess2>::
+ doContinuous(TParticle& particle, TTrack& vT,
+ [[maybe_unused]] ContinuousProcessIndex const limitId) {
ProcessReturn ret = ProcessReturn::Ok;
- if constexpr (std::is_base_of_v<ContinuousProcess<process1_type>, process1_type> ||
- t1ProcSeq) {
- ret |= A_.doContinuous(particle, vT);
+ if constexpr (t1ProcSeq) {
+ ret |= A_.doContinuous(particle, vT, limitId);
+ } else if constexpr (is_continuous_process_v<process1_type>) {
+ ret |=
+ A_.doContinuous(particle, vT, limitId == ContinuousProcessIndex(IndexProcess1));
}
- if constexpr (std::is_base_of_v<ContinuousProcess<process2_type>, process2_type> ||
- t2ProcSeq) {
- if (!isAbsorbed(ret)) { ret |= B_.doContinuous(particle, vT); }
+
+ if constexpr (t2ProcSeq) {
+ if (!isAbsorbed(ret)) { ret |= B_.doContinuous(particle, vT, limitId); }
+ } else if constexpr (is_continuous_process_v<process2_type>) {
+ if (!isAbsorbed(ret)) {
+ ret |= B_.doContinuous(particle, vT,
+ limitId == ContinuousProcessIndex(IndexProcess2));
+ }
}
+
return ret;
}
- template <typename TProcess1, typename TProcess2, int IndexStart, int IndexProcess1, int IndexProcess2>
+ template <typename TProcess1, typename TProcess2, int IndexStart, int IndexProcess1,
+ int IndexProcess2>
template <typename TSecondaries>
- void ProcessSequence<TProcess1, TProcess2, IndexStart, IndexProcess1, IndexProcess2>::doSecondaries(TSecondaries& vS) {
+ void ProcessSequence<TProcess1, TProcess2, IndexStart, IndexProcess1,
+ IndexProcess2>::doSecondaries(TSecondaries& vS) {
if constexpr (std::is_base_of_v<SecondariesProcess<process1_type>, process1_type> ||
t1ProcSeq) {
A_.doSecondaries(vS);
@@ -75,8 +93,10 @@ namespace corsika {
}
}
- template <typename TProcess1, typename TProcess2, int IndexStart, int IndexProcess1, int IndexProcess2>
- bool ProcessSequence<TProcess1, TProcess2, IndexStart, IndexProcess1, IndexProcess2>::checkStep() {
+ template <typename TProcess1, typename TProcess2, int IndexStart, int IndexProcess1,
+ int IndexProcess2>
+ bool ProcessSequence<TProcess1, TProcess2, IndexStart, IndexProcess1,
+ IndexProcess2>::checkStep() {
bool ret = false;
if constexpr (std::is_base_of_v<StackProcess<process1_type>, process1_type> ||
(t1ProcSeq && !t1SwitchProcSeq)) {
@@ -89,9 +109,11 @@ namespace corsika {
return ret;
}
- template <typename TProcess1, typename TProcess2, int IndexStart, int IndexProcess1, int IndexProcess2>
+ template <typename TProcess1, typename TProcess2, int IndexStart, int IndexProcess1,
+ int IndexProcess2>
template <typename TStack>
- void ProcessSequence<TProcess1, TProcess2, IndexStart, IndexProcess1, IndexProcess2>::doStack(TStack& stack) {
+ void ProcessSequence<TProcess1, TProcess2, IndexStart, IndexProcess1,
+ IndexProcess2>::doStack(TStack& stack) {
if constexpr (std::is_base_of_v<StackProcess<process1_type>, process1_type> ||
(t1ProcSeq && !t1SwitchProcSeq)) {
if (A_.checkStep()) { A_.doStack(stack); }
@@ -102,30 +124,52 @@ namespace corsika {
}
}
- template <typename TProcess1, typename TProcess2, int IndexStart, int IndexProcess1, int IndexProcess2>
+ /**
+ * Calculate the maximum allowed length of the next tracking step, based on all
+ * ContinuousProcess-es
+ *
+ * The maximum allowed step length is the minimum of the allowed track lenght over all
+ * ContinuousProcess-es in the ProcessSequence.
+ *
+ * \return: ContinuousProcessStepLength which contains the step length itself in
+ * LengthType, and a unique identifier of the related ContinuousProcess.
+ **/
+ template <typename TProcess1, typename TProcess2, int IndexStart, int IndexProcess1,
+ int IndexProcess2>
template <typename TParticle, typename TTrack>
- LengthType ProcessSequence<TProcess1, TProcess2, IndexStart, IndexProcess1, IndexProcess2>::getMaxStepLength(TParticle& particle,
- TTrack& vTrack) {
- LengthType max_length = // if no other process in the sequence implements it
- std::numeric_limits<double>::infinity() * meter;
+ ContinuousProcessStepLength
+ ProcessSequence<TProcess1, TProcess2, IndexStart, IndexProcess1,
+ IndexProcess2>::getMaxStepLength(TParticle& particle, TTrack& vTrack) {
+ // if no other process in the sequence implements it
+ ContinuousProcessStepLength max_length(std::numeric_limits<double>::infinity() *
+ meter);
- if constexpr (std::is_base_of_v<ContinuousProcess<process1_type>, process1_type> ||
- t1ProcSeq) {
- LengthType const len = A_.getMaxStepLength(particle, vTrack);
- max_length = std::min(max_length, len);
+ if constexpr (t1ProcSeq) {
+ ContinuousProcessStepLength const step = A_.getMaxStepLength(particle, vTrack);
+ max_length = std::min(max_length, step);
+ } else if constexpr (is_continuous_process_v<process1_type>) {
+ ContinuousProcessStepLength const step(A_.getMaxStepLength(particle, vTrack),
+ ContinuousProcessIndex(IndexProcess1));
+ max_length = std::min(max_length, step);
}
- if constexpr (std::is_base_of_v<ContinuousProcess<process2_type>, process2_type> ||
- t2ProcSeq) {
- LengthType const len = B_.getMaxStepLength(particle, vTrack);
- max_length = std::min(max_length, len);
+ // return ContinuousProcessStepLength(std::min(max_length, len), IndexProcess1);
+ if constexpr (t2ProcSeq) {
+ ContinuousProcessStepLength const step = B_.getMaxStepLength(particle, vTrack);
+ max_length = std::min(max_length, step);
+ } else if constexpr (is_continuous_process_v<process2_type>) {
+ ContinuousProcessStepLength const step(B_.getMaxStepLength(particle, vTrack),
+ ContinuousProcessIndex(IndexProcess2));
+ max_length = std::min(max_length, step);
}
return max_length;
}
- template <typename TProcess1, typename TProcess2, int IndexStart, int IndexProcess1, int IndexProcess2>
+ template <typename TProcess1, typename TProcess2, int IndexStart, int IndexProcess1,
+ int IndexProcess2>
template <typename TParticle>
- InverseGrammageType ProcessSequence<TProcess1, TProcess2, IndexStart, IndexProcess1, IndexProcess2>::getInverseInteractionLength(
- TParticle&& particle) {
+ InverseGrammageType
+ ProcessSequence<TProcess1, TProcess2, IndexStart, IndexProcess1,
+ IndexProcess2>::getInverseInteractionLength(TParticle&& particle) {
InverseGrammageType tot = 0 * meter * meter / gram; // default value
@@ -140,11 +184,14 @@ namespace corsika {
return tot;
}
- template <typename TProcess1, typename TProcess2, int IndexStart, int IndexProcess1, int IndexProcess2>
+ template <typename TProcess1, typename TProcess2, int IndexStart, int IndexProcess1,
+ int IndexProcess2>
template <typename TSecondaryView>
- inline ProcessReturn ProcessSequence<TProcess1, TProcess2, IndexStart, IndexProcess1, IndexProcess2>::selectInteraction(
- TSecondaryView& view, [[maybe_unused]] InverseGrammageType lambda_inv_select,
- [[maybe_unused]] InverseGrammageType lambda_inv_sum) {
+ inline ProcessReturn
+ ProcessSequence<TProcess1, TProcess2, IndexStart, IndexProcess1, IndexProcess2>::
+ selectInteraction(TSecondaryView& view,
+ [[maybe_unused]] InverseGrammageType lambda_inv_select,
+ [[maybe_unused]] InverseGrammageType lambda_inv_sum) {
// TODO: add check for lambda_inv_select>lambda_inv_tot
@@ -182,10 +229,12 @@ namespace corsika {
return ProcessReturn::Ok;
}
- template <typename TProcess1, typename TProcess2, int IndexStart, int IndexProcess1, int IndexProcess2>
+ template <typename TProcess1, typename TProcess2, int IndexStart, int IndexProcess1,
+ int IndexProcess2>
template <typename TParticle>
- inline InverseTimeType ProcessSequence<TProcess1, TProcess2, IndexStart, IndexProcess1, IndexProcess2>::getInverseLifetime(
- TParticle&& particle) {
+ inline InverseTimeType
+ ProcessSequence<TProcess1, TProcess2, IndexStart, IndexProcess1,
+ IndexProcess2>::getInverseLifetime(TParticle&& particle) {
InverseTimeType tot = 0 / second; // default value
@@ -200,12 +249,15 @@ namespace corsika {
return tot;
}
- template <typename TProcess1, typename TProcess2, int IndexStart, int IndexProcess1, int IndexProcess2>
+ template <typename TProcess1, typename TProcess2, int IndexStart, int IndexProcess1,
+ int IndexProcess2>
// select decay process
template <typename TSecondaryView>
- inline ProcessReturn ProcessSequence<TProcess1, TProcess2, IndexStart, IndexProcess1, IndexProcess2>::selectDecay(
- TSecondaryView& view, [[maybe_unused]] InverseTimeType decay_inv_select,
- [[maybe_unused]] InverseTimeType decay_inv_sum) {
+ inline ProcessReturn ProcessSequence<
+ TProcess1, TProcess2, IndexStart, IndexProcess1,
+ IndexProcess2>::selectDecay(TSecondaryView& view,
+ [[maybe_unused]] InverseTimeType decay_inv_select,
+ [[maybe_unused]] InverseTimeType decay_inv_sum) {
// TODO: add check for decay_inv_select>decay_inv_tot
diff --git a/corsika/detail/framework/process/SwitchProcessSequence.inl b/corsika/detail/framework/process/SwitchProcessSequence.inl
index 413caae4a..2eadc8782 100644
--- a/corsika/detail/framework/process/SwitchProcessSequence.inl
+++ b/corsika/detail/framework/process/SwitchProcessSequence.inl
@@ -12,6 +12,8 @@
#include <corsika/framework/process/ProcessTraits.hpp>
#include <corsika/framework/process/BoundaryCrossingProcess.hpp>
#include <corsika/framework/process/ContinuousProcess.hpp>
+#include <corsika/framework/process/ContinuousProcessStepLength.hpp>
+#include <corsika/framework/process/ContinuousProcessIndex.hpp>
#include <corsika/framework/process/DecayProcess.hpp>
#include <corsika/framework/process/InteractionProcess.hpp>
#include <corsika/framework/process/ProcessReturn.hpp>
@@ -25,11 +27,14 @@
namespace corsika {
- template <typename TProcess1, typename TProcess2, typename TSelect,
- int IndexStart, int IndexProcess1, int IndexProcess2>
+ template <typename TProcess1, typename TProcess2, typename TSelect, int IndexStart,
+ int IndexProcess1, int IndexProcess2>
template <typename TParticle, typename TVTNType>
- ProcessReturn SwitchProcessSequence<TProcess1, TProcess2, TSelect, IndexStart, IndexProcess1, IndexProcess2>::doBoundaryCrossing(
- TParticle& particle, TVTNType const& from, TVTNType const& to) {
+ ProcessReturn
+ SwitchProcessSequence<TProcess1, TProcess2, TSelect, IndexStart, IndexProcess1,
+ IndexProcess2>::doBoundaryCrossing(TParticle& particle,
+ TVTNType const& from,
+ TVTNType const& to) {
switch (select_(particle)) {
case SwitchResult::First: {
if constexpr (std::is_base_of_v<BoundaryCrossingProcess<process1_type>,
@@ -51,24 +56,27 @@ namespace corsika {
return ProcessReturn::Ok;
}
- template <typename TProcess1, typename TProcess2, typename TSelect, int IndexStart, int IndexProcess1, int IndexProcess2>
+ template <typename TProcess1, typename TProcess2, typename TSelect, int IndexStart,
+ int IndexProcess1, int IndexProcess2>
template <typename TParticle, typename TTrack>
- inline ProcessReturn SwitchProcessSequence<TProcess1, TProcess2, TSelect, IndexStart, IndexProcess1, IndexProcess2>::doContinuous(
- TParticle& particle, TTrack& vT) {
+ inline ProcessReturn SwitchProcessSequence<
+ TProcess1, TProcess2, TSelect, IndexStart, IndexProcess1,
+ IndexProcess2>::doContinuous(TParticle& particle, TTrack& vT,
+ ContinuousProcessIndex const idLimit) {
switch (select_(particle)) {
case SwitchResult::First: {
- if constexpr (std::is_base_of_v<ContinuousProcess<process1_type>,
- process1_type> ||
- t1ProcSeq) {
- return A_.doContinuous(particle, vT);
+ if constexpr (t1ProcSeq) { return A_.doContinuous(particle, vT, idLimit); }
+ if constexpr (is_continuous_process_v<process1_type>) {
+ return A_.doContinuous(particle, vT,
+ idLimit == ContinuousProcessIndex(IndexProcess1));
}
break;
}
case SwitchResult::Second: {
- if constexpr (std::is_base_of_v<ContinuousProcess<process2_type>,
- process2_type> ||
- t2ProcSeq) {
- return B_.doContinuous(particle, vT);
+ if constexpr (t2ProcSeq) { return B_.doContinuous(particle, vT, idLimit); }
+ if constexpr (is_continuous_process_v<process2_type>) {
+ return B_.doContinuous(particle, vT,
+ idLimit == ContinuousProcessIndex(IndexProcess2));
}
break;
}
@@ -76,10 +84,12 @@ namespace corsika {
return ProcessReturn::Ok;
}
- template <typename TProcess1, typename TProcess2, typename TSelect, int IndexStart, int IndexProcess1, int IndexProcess2>
+ template <typename TProcess1, typename TProcess2, typename TSelect, int IndexStart,
+ int IndexProcess1, int IndexProcess2>
template <typename TSecondaries>
- inline void SwitchProcessSequence<TProcess1, TProcess2, TSelect, IndexStart, IndexProcess1, IndexProcess2>::doSecondaries(
- TSecondaries& vS) {
+ inline void
+ SwitchProcessSequence<TProcess1, TProcess2, TSelect, IndexStart, IndexProcess1,
+ IndexProcess2>::doSecondaries(TSecondaries& vS) {
const auto& particle = vS.parent();
switch (select_(particle)) {
case SwitchResult::First: {
@@ -101,39 +111,42 @@ namespace corsika {
}
}
- template <typename TProcess1, typename TProcess2, typename TSelect, int IndexStart, int IndexProcess1, int IndexProcess2>
+ template <typename TProcess1, typename TProcess2, typename TSelect, int IndexStart,
+ int IndexProcess1, int IndexProcess2>
template <typename TParticle, typename TTrack>
- inline LengthType SwitchProcessSequence<TProcess1, TProcess2,
- TSelect, IndexStart, IndexProcess1, IndexProcess2>::getMaxStepLength(TParticle& particle,
- TTrack& vTrack) {
+ inline ContinuousProcessStepLength
+ SwitchProcessSequence<TProcess1, TProcess2, TSelect, IndexStart, IndexProcess1,
+ IndexProcess2>::getMaxStepLength(TParticle& particle,
+ TTrack& vTrack) {
switch (select_(particle)) {
case SwitchResult::First: {
- if constexpr (std::is_base_of_v<ContinuousProcess<process1_type>,
- process1_type> ||
- t1ProcSeq) {
- return A_.getMaxStepLength(particle, vTrack);
+ if constexpr (t1ProcSeq) { return A_.getMaxStepLength(particle, vTrack); }
+ if constexpr (is_continuous_process_v<process1_type>) {
+ return ContinuousProcessStepLength(A_.getMaxStepLength(particle, vTrack),
+ ContinuousProcessIndex(IndexProcess1));
}
break;
}
case SwitchResult::Second: {
- if constexpr (std::is_base_of_v<ContinuousProcess<process2_type>,
- process2_type> ||
- t2ProcSeq) {
- return B_.getMaxStepLength(particle, vTrack);
+ if constexpr (t2ProcSeq) { return B_.getMaxStepLength(particle, vTrack); }
+ if constexpr (is_continuous_process_v<process2_type>) {
+ return ContinuousProcessStepLength(B_.getMaxStepLength(particle, vTrack),
+ ContinuousProcessIndex(IndexProcess2));
}
break;
}
}
// if no other process in the sequence implements it
- return std::numeric_limits<double>::infinity() * meter;
+ return ContinuousProcessStepLength(std::numeric_limits<double>::infinity() * meter);
}
- template <typename TProcess1, typename TProcess2, typename TSelect, int IndexStart, int IndexProcess1, int IndexProcess2>
+ template <typename TProcess1, typename TProcess2, typename TSelect, int IndexStart,
+ int IndexProcess1, int IndexProcess2>
template <typename TParticle>
- inline InverseGrammageType
- SwitchProcessSequence<TProcess1, TProcess2, TSelect, IndexStart, IndexProcess1, IndexProcess2>::getInverseInteractionLength(
- TParticle&& particle) {
+ inline InverseGrammageType SwitchProcessSequence<
+ TProcess1, TProcess2, TSelect, IndexStart, IndexProcess1,
+ IndexProcess2>::getInverseInteractionLength(TParticle&& particle) {
switch (select_(particle)) {
case SwitchResult::First: {
@@ -156,12 +169,14 @@ namespace corsika {
return 0 * meter * meter / gram; // default value
}
- template <typename TProcess1, typename TProcess2, typename TSelect, int IndexStart, int IndexProcess1, int IndexProcess2>
+ template <typename TProcess1, typename TProcess2, typename TSelect, int IndexStart,
+ int IndexProcess1, int IndexProcess2>
template <typename TSecondaryView>
- inline ProcessReturn
- SwitchProcessSequence<TProcess1, TProcess2, TSelect, IndexStart, IndexProcess1, IndexProcess2>::selectInteraction(
- TSecondaryView& view, [[maybe_unused]] InverseGrammageType lambda_inv_select,
- [[maybe_unused]] InverseGrammageType lambda_inv_sum) {
+ inline ProcessReturn SwitchProcessSequence<TProcess1, TProcess2, TSelect, IndexStart,
+ IndexProcess1, IndexProcess2>::
+ selectInteraction(TSecondaryView& view,
+ [[maybe_unused]] InverseGrammageType lambda_inv_select,
+ [[maybe_unused]] InverseGrammageType lambda_inv_sum) {
switch (select_(view.parent())) {
case SwitchResult::First: {
if constexpr (t1ProcSeq) {
@@ -204,11 +219,12 @@ namespace corsika {
return ProcessReturn::Ok;
}
- template <typename TProcess1, typename TProcess2, typename TSelect, int IndexStart, int IndexProcess1, int IndexProcess2>
+ template <typename TProcess1, typename TProcess2, typename TSelect, int IndexStart,
+ int IndexProcess1, int IndexProcess2>
template <typename TParticle>
inline InverseTimeType
- SwitchProcessSequence<TProcess1, TProcess2, TSelect, IndexStart, IndexProcess1, IndexProcess2>::getInverseLifetime(
- TParticle&& particle) {
+ SwitchProcessSequence<TProcess1, TProcess2, TSelect, IndexStart, IndexProcess1,
+ IndexProcess2>::getInverseLifetime(TParticle&& particle) {
switch (select_(particle)) {
case SwitchResult::First: {
@@ -230,12 +246,15 @@ namespace corsika {
return 0 / second; // default value
}
- template <typename TProcess1, typename TProcess2, typename TSelect, int IndexStart, int IndexProcess1, int IndexProcess2>
+ template <typename TProcess1, typename TProcess2, typename TSelect, int IndexStart,
+ int IndexProcess1, int IndexProcess2>
// select decay process
template <typename TSecondaryView>
- inline ProcessReturn SwitchProcessSequence<TProcess1, TProcess2, TSelect, IndexStart, IndexProcess1, IndexProcess2>::selectDecay(
- TSecondaryView& view, [[maybe_unused]] InverseTimeType decay_inv_select,
- [[maybe_unused]] InverseTimeType decay_inv_sum) {
+ inline ProcessReturn SwitchProcessSequence<
+ TProcess1, TProcess2, TSelect, IndexStart, IndexProcess1,
+ IndexProcess2>::selectDecay(TSecondaryView& view,
+ [[maybe_unused]] InverseTimeType decay_inv_select,
+ [[maybe_unused]] InverseTimeType decay_inv_sum) {
switch (select_(view.parent())) {
case SwitchResult::First: {
if constexpr (t1ProcSeq) {
@@ -278,11 +297,23 @@ namespace corsika {
return ProcessReturn::Ok;
}
+ /*
+ /// traits marker to identify objectas ProcessSequence
+ template <typename TProcess1, typename TProcess2, typename TSelect>
+ struct is_process_sequence<ProcessSequence<typename std::decay_t<TProcess1>,
+ typename std::decay_t<TProcess2>,
+ typename std::decay_t<TSelect>>>
+ : std::true_type {};
+ */
/// traits marker to identify objectas ProcessSequence
template <typename TProcess1, typename TProcess2, typename TSelect>
struct is_process_sequence<SwitchProcessSequence<TProcess1, TProcess2, TSelect>>
: std::true_type {};
+ template <typename TProcess1, typename TProcess2, typename TSelect>
+ struct is_process_sequence<SwitchProcessSequence<TProcess1, TProcess2, TSelect>&>
+ : std::true_type {};
+
/// traits marker to identify objectas SwitchProcessSequence
template <typename TProcess1, typename TProcess2, typename TSelect>
struct is_switch_process_sequence<SwitchProcessSequence<TProcess1, TProcess2, TSelect>>
diff --git a/corsika/detail/modules/LongitudinalProfile.inl b/corsika/detail/modules/LongitudinalProfile.inl
index d18813c26..df6bb471d 100644
--- a/corsika/detail/modules/LongitudinalProfile.inl
+++ b/corsika/detail/modules/LongitudinalProfile.inl
@@ -27,7 +27,7 @@ namespace corsika {
template <typename TParticle, typename TTrack>
ProcessReturn LongitudinalProfile::doContinuous(TParticle const& vP,
- TTrack const& vTrack) {
+ TTrack const& vTrack, bool const) {
auto const pid = vP.getPID();
GrammageType const grammageStart = shower_axis_.getProjectedX(vTrack.getPosition(0));
diff --git a/corsika/detail/modules/ObservationPlane.inl b/corsika/detail/modules/ObservationPlane.inl
index 526fb273a..74c29b336 100644
--- a/corsika/detail/modules/ObservationPlane.inl
+++ b/corsika/detail/modules/ObservationPlane.inl
@@ -28,34 +28,30 @@ namespace corsika {
}
ProcessReturn ObservationPlane::doContinuous(
- corsika::setup::Stack::particle_type& particle,
- corsika::setup::Trajectory& trajectory,
+ corsika::setup::Stack::particle_type& particle, corsika::setup::Trajectory&,
bool const stepLimit) {
- auto const* volumeNode = particle.getNode();
-
- if (!stepLimit) {
- return ProcessReturn::Ok;
- }
-
- HEPEnergyType const energy = particle.getEnergy();
- TimeType const timeOfIntersection = particle.getTime();
- Point const pointOfIntersection = particle.getPosition();
- Vector const displacement = pointOfIntersection - plane_.getCenter();
-
- outputStream_ << static_cast<int>(get_PDG(particle.getPID())) << ' ' << energy / 1_eV
- << ' ' << displacement.dot(xAxis_) / 1_m << ' '
- << displacement.dot(yAxis_) / 1_m
- << (pointOfIntersection - plane_.getCenter()).getNorm() / 1_m
- << '\n';
-
- if (deleteOnHit_) {
- count_ground_++;
- energy_ground_ += energy;
- particle.erase();
- return ProcessReturn::ParticleAbsorbed;
- } else {
- return ProcessReturn::Ok;
+ /*
+ The current step did not yet reach the ObservationPlane, do nothing now and wait:
+ */
+ if (!stepLimit) { return ProcessReturn::Ok; }
+
+ HEPEnergyType const energy = particle.getEnergy();
+ Point const pointOfIntersection = particle.getPosition();
+ Vector const displacement = pointOfIntersection - plane_.getCenter();
+
+ outputStream_ << static_cast<int>(get_PDG(particle.getPID())) << ' ' << energy / 1_eV
+ << ' ' << displacement.dot(xAxis_) / 1_m << ' '
+ << displacement.dot(yAxis_) / 1_m
+ << (pointOfIntersection - plane_.getCenter()).getNorm() / 1_m << '\n';
+
+ if (deleteOnHit_) {
+ count_ground_++;
+ energy_ground_ += energy;
+ particle.erase();
+ return ProcessReturn::ParticleAbsorbed;
+ } else {
+ return ProcessReturn::Ok;
}
}
@@ -66,18 +62,22 @@ namespace corsika {
auto const& volumeNode = particle.getNode();
typedef typename std::remove_const_t<
- std::remove_reference_t<decltype(volumeNode->getModelProperties())>>
- medium_type;
-
- Intersections const intersection = setup::Tracking::intersect<corsika::setup::Stack::particle_type, medium_type>(particle, plane_,
- volumeNode->getModelProperties());
+ std::remove_reference_t<decltype(volumeNode->getModelProperties())>>
+ medium_type;
+ Intersections const intersection =
+ setup::Tracking::intersect<corsika::setup::Stack::particle_type, medium_type>(
+ particle, plane_, volumeNode->getModelProperties());
TimeType const timeOfIntersection = intersection.getEntry();
-
- if (timeOfIntersection < TimeType::zero()) {
+
+ CORSIKA_LOG_TRACE("particle={}, pos={}, dir={}, plane={}, timeOfIntersection={}",
+ particle.asString(), particle.getPosition(),
+ particle.getDirection(), plane_.asString(), timeOfIntersection);
+
+ if (timeOfIntersection < TimeType::zero()) {
return std::numeric_limits<double>::infinity() * 1_m;
- }
+ }
if (timeOfIntersection > trajectory.getDuration()) {
return std::numeric_limits<double>::infinity() * 1_m;
}
@@ -85,9 +85,8 @@ namespace corsika {
double const fractionOfIntersection = timeOfIntersection / trajectory.getDuration();
auto const pointOfIntersection = trajectory.getPosition(fractionOfIntersection);
- auto dist = (trajectory.getPosition(0) - pointOfIntersection).getNorm() * 1.001; // make max. step length a bit longer, to assure crossing
- CORSIKA_LOG_TRACE("ObservationPlane: getMaxStepLength l={} m",
- dist / 1_m);
+ auto dist = (trajectory.getPosition(0) - pointOfIntersection).getNorm();
+ CORSIKA_LOG_TRACE("ObservationPlane: getMaxStepLength l={} m", dist / 1_m);
return dist;
}
diff --git a/corsika/detail/modules/ParticleCut.inl b/corsika/detail/modules/ParticleCut.inl
index 7a22b4d71..3fd882225 100644
--- a/corsika/detail/modules/ParticleCut.inl
+++ b/corsika/detail/modules/ParticleCut.inl
@@ -145,7 +145,7 @@ namespace corsika {
}
ProcessReturn ParticleCut::doContinuous(corsika::setup::Stack::particle_type& particle,
- corsika::setup::Trajectory const&) {
+ corsika::setup::Trajectory const&, bool const) {
CORSIKA_LOG_TRACE("ParticleCut::DoContinuous");
if (checkCutParticle(particle)) {
CORSIKA_LOG_TRACE("removing during continuous");
diff --git a/corsika/detail/modules/TrackWriter.inl b/corsika/detail/modules/TrackWriter.inl
index 27b47cf65..c6f72edf4 100644
--- a/corsika/detail/modules/TrackWriter.inl
+++ b/corsika/detail/modules/TrackWriter.inl
@@ -31,7 +31,8 @@ namespace corsika {
}
template <typename TParticle, typename TTrack>
- ProcessReturn TrackWriter::doContinuous(const TParticle& vP, const TTrack& vT) {
+ ProcessReturn TrackWriter::doContinuous(const TParticle& vP, const TTrack& vT,
+ bool const) {
auto const start = vT.getPosition(0).getCoordinates();
auto const delta = vT.getPosition(1).getCoordinates() - start;
auto const pdg = static_cast<int>(get_PDG(vP.getPID()));
diff --git a/corsika/detail/modules/energy_loss/BetheBlochPDG.inl b/corsika/detail/modules/energy_loss/BetheBlochPDG.inl
index 89fe81c7b..cbabb0892 100644
--- a/corsika/detail/modules/energy_loss/BetheBlochPDG.inl
+++ b/corsika/detail/modules/energy_loss/BetheBlochPDG.inl
@@ -140,7 +140,7 @@ namespace corsika {
}
ProcessReturn BetheBlochPDG::doContinuous(setup::Stack::particle_type& p,
- setup::Trajectory const& t) {
+ setup::Trajectory const& t, bool const) {
if (p.getChargeNumber() == 0) return ProcessReturn::Ok;
GrammageType const dX =
@@ -174,7 +174,7 @@ namespace corsika {
// slightly smaller than emCut since, either this Step is limited
// by energy_lim, then the particle is stopped in a very short
// range (before doing anythin else) and is then removed
- // instantly. The exact position where it reaches emCut is not
+ // instantly. The exact 3D position where it reaches emCut is not very
// important, the important fact is that its E_kin is zero
// afterwards.
//
diff --git a/corsika/detail/modules/proposal/ContinuousProcess.inl b/corsika/detail/modules/proposal/ContinuousProcess.inl
index 8bfbd705b..26f604c99 100644
--- a/corsika/detail/modules/proposal/ContinuousProcess.inl
+++ b/corsika/detail/modules/proposal/ContinuousProcess.inl
@@ -88,7 +88,7 @@ namespace corsika::proposal {
template <>
ProcessReturn ContinuousProcess::doContinuous(setup::Stack::particle_type& vP,
- setup::Trajectory const& vT) {
+ setup::Trajectory const& vT, bool const) {
if (!canInteract(vP.getPID())) return ProcessReturn::Ok;
if (vT.getLength() == 0_m) return ProcessReturn::Ok;
diff --git a/corsika/detail/modules/tracking/TrackingStraight.inl b/corsika/detail/modules/tracking/TrackingStraight.inl
index 2c284b4d0..5869601a5 100644
--- a/corsika/detail/modules/tracking/TrackingStraight.inl
+++ b/corsika/detail/modules/tracking/TrackingStraight.inl
@@ -92,6 +92,8 @@ namespace corsika::tracking_line {
auto const n = plane.getNormal();
auto const c = n.dot(velocity);
+ CORSIKA_LOG_TRACE("c={}, delta={}, momentum={}", c, delta, particle.getMomentum());
+
return Intersections(c.magnitude() == 0
? std::numeric_limits<TimeType::value_type>::infinity() * 1_s
: n.dot(delta) / c);
diff --git a/corsika/framework/geometry/Plane.hpp b/corsika/framework/geometry/Plane.hpp
index 9b2eacade..7218e6b35 100644
--- a/corsika/framework/geometry/Plane.hpp
+++ b/corsika/framework/geometry/Plane.hpp
@@ -12,6 +12,8 @@
#include <corsika/framework/geometry/Point.hpp>
#include <corsika/framework/geometry/Vector.hpp>
+#include <string>
+
namespace corsika {
class Plane {
@@ -32,6 +34,8 @@ namespace corsika {
DimLessVec const& getNormal() const;
+ std::string asString() const;
+
public:
Point const center_;
DimLessVec const normal_;
diff --git a/corsika/framework/geometry/QuantityVector.hpp b/corsika/framework/geometry/QuantityVector.hpp
index 235237afe..eb6c2f11e 100644
--- a/corsika/framework/geometry/QuantityVector.hpp
+++ b/corsika/framework/geometry/QuantityVector.hpp
@@ -101,7 +101,7 @@ namespace corsika {
template <typename T>
friend class corsika::Vector;
template <typename TDim>
- friend std::ostream& operator<<(std::ostream& os, QuantityVector<TDim> qv);
+ friend std::ostream& operator<<(std::ostream& os, QuantityVector<TDim> const& qv);
protected:
Eigen::Vector3d
@@ -114,7 +114,7 @@ namespace corsika {
template <typename TDimension>
inline std::ostream& operator<<(std::ostream& os,
- corsika::QuantityVector<TDimension> const qv);
+ corsika::QuantityVector<TDimension> const& qv);
} // namespace corsika
diff --git a/corsika/framework/process/BaseProcess.hpp b/corsika/framework/process/BaseProcess.hpp
index e91651d12..55572291c 100644
--- a/corsika/framework/process/BaseProcess.hpp
+++ b/corsika/framework/process/BaseProcess.hpp
@@ -8,6 +8,10 @@
#pragma once
+//#include <corsika/framework/process/ProcessTraits.hpp>
+
+#include <type_traits>
+
namespace corsika {
class TDerived; // fwd decl
diff --git a/corsika/framework/process/ContinuousProcess.hpp b/corsika/framework/process/ContinuousProcess.hpp
index 065d4c234..521cb0ee2 100644
--- a/corsika/framework/process/ContinuousProcess.hpp
+++ b/corsika/framework/process/ContinuousProcess.hpp
@@ -8,6 +8,7 @@
#pragma once
+#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/process/BaseProcess.hpp>
#include <corsika/framework/process/ProcessReturn.hpp>
#include <corsika/framework/process/ProcessTraits.hpp>
@@ -29,11 +30,23 @@ namespace corsika {
protected:
public:
// here starts the interface part
- // -> enforce TDerived to implement DoContinuous...
+ /**
+ * Applies the effects of this ContinuousProcess on a Particle on a Track.
+ *
+ * Note, the stepLimit is a flag, if this particular process was responsible for the
+ * track-length limit. This can be used by the process to trigger activity.
+ *
+ * \todo -> enforce TDerived to implement doContinuous...
+ **/
template <typename TParticle, typename TTrack>
ProcessReturn doContinuous(TParticle&, TTrack const&, bool const stepLimit) const;
- // -> enforce TDerived to implement MaxStepLength...
+ /**
+ * Calculates/returns a possible step length limitation of this continuousprocess.
+ *
+ *
+ * \todo -> enforce TDerived to implement getMaxStepLength...
+ **/
template <typename TParticle, typename TTrack>
LengthType getMaxStepLength(TParticle const& p, TTrack const& track) const;
};
@@ -41,12 +54,20 @@ namespace corsika {
/**
* ProcessTraits specialization
**/
- /*
- template <typename TProcess, int N>
- struct count_continuous<TProcess, N,
- typename std::enable_if_t<std::is_base_of_v<ContinuousProcess<typename std::decay_t<TProcess>>,
- typename std::decay_t<TProcess>>>> {
- enum { count = N+1 };
- };
- */
+ template <typename TProcess>
+ struct is_continuous_process<
+ TProcess, std::enable_if_t<
+ std::is_base_of_v<ContinuousProcess<typename std::decay_t<TProcess>>,
+ typename std::decay_t<TProcess>>>>
+ : std::true_type {};
+
+ template <typename TProcess, int N>
+ struct count_continuous<TProcess, N,
+ typename std::enable_if_t<is_continuous_process_v<TProcess>>> {
+ // std::is_base_of_v<
+ // ContinuousProcess<typename std::decay_t<TProcess>>,
+ // typename std::decay_t<TProcess>>>> {
+ enum { count = N + 1 };
+ };
+
} // namespace corsika
diff --git a/corsika/framework/process/ContinuousProcessIndex.hpp b/corsika/framework/process/ContinuousProcessIndex.hpp
new file mode 100644
index 000000000..4cba88fde
--- /dev/null
+++ b/corsika/framework/process/ContinuousProcessIndex.hpp
@@ -0,0 +1,26 @@
+/*
+ * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * 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.
+ */
+
+#pragma once
+
+namespace corsika {
+
+ class ContinuousProcessIndex {
+ public:
+ ContinuousProcessIndex(int const id)
+ : id_(id) {}
+ void setIndex(int const id) { id_ = id; }
+ int getIndex() const { return id_; }
+ bool operator==(ContinuousProcessIndex const v) const { return id_ == v.id_; }
+ bool operator!=(ContinuousProcessIndex const v) const { return id_ != v.id_; }
+
+ private:
+ int id_;
+ };
+
+} // namespace corsika
diff --git a/corsika/framework/process/ContinuousProcessStepLength.hpp b/corsika/framework/process/ContinuousProcessStepLength.hpp
new file mode 100644
index 000000000..32b921857
--- /dev/null
+++ b/corsika/framework/process/ContinuousProcessStepLength.hpp
@@ -0,0 +1,43 @@
+/*
+ * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * 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.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/process/ContinuousProcessIndex.hpp>
+
+namespace corsika {
+
+ /**
+ * To store step length in LengthType and unique index in ProcessSequence of shortest
+ * step ContinuousProcess.
+ *
+ **/
+
+ class ContinuousProcessStepLength {
+ public:
+ ContinuousProcessStepLength(LengthType const step)
+ : step_(step)
+ , id_(0) {}
+ ContinuousProcessStepLength(LengthType const step, ContinuousProcessIndex const id)
+ : step_(step)
+ , id_(id) {}
+ void setLength(LengthType const step) { step_ = step; }
+ LengthType getLength() const { return step_; }
+ void setIndex(ContinuousProcessIndex const id) { id_ = id; }
+ ContinuousProcessIndex getIndex() const { return id_; }
+ operator ContinuousProcessIndex() const { return id_; }
+ operator LengthType() const { return step_; }
+ bool operator<(ContinuousProcessStepLength const& r) const { return step_ < r.step_; }
+
+ private:
+ LengthType step_;
+ ContinuousProcessIndex id_;
+ };
+
+} // namespace corsika
diff --git a/corsika/framework/process/ProcessSequence.hpp b/corsika/framework/process/ProcessSequence.hpp
index 594a6637b..d344673c6 100644
--- a/corsika/framework/process/ProcessSequence.hpp
+++ b/corsika/framework/process/ProcessSequence.hpp
@@ -16,6 +16,8 @@
#include <corsika/framework/process/ProcessTraits.hpp>
#include <corsika/framework/process/BoundaryCrossingProcess.hpp>
#include <corsika/framework/process/ContinuousProcess.hpp>
+#include <corsika/framework/process/ContinuousProcessIndex.hpp>
+#include <corsika/framework/process/ContinuousProcessStepLength.hpp>
#include <corsika/framework/process/DecayProcess.hpp>
#include <corsika/framework/process/InteractionProcess.hpp>
#include <corsika/framework/process/ProcessReturn.hpp>
@@ -26,6 +28,28 @@
namespace corsika {
+ /*
+ template <typename TProcess1, typename TProcess2 = NullModel, int IndexStart = 0,
+ int IndexProcess1 = count_continuous<TProcess1>::count,
+ int IndexProcess2 =
+ count_continuous<TProcess1, count_continuous<TProcess2>::count>::count>
+ class ProcessSequence;
+ */
+
+ /**
+ * traits class to count ContinuousProcess-es, specialized for ProcessSequence-es
+ **/
+ // template <typename TProcess1, typename TProcess2, int N>
+ // struct count_continuous<ProcessSequence<TProcess1, TProcess2>, N> {
+ // enum { count = N + ProcessSequence<TProcess1,TProcess2>::nContinuous };
+ //};
+
+ template <typename TProcess, int N>
+ struct count_continuous<TProcess, N,
+ typename std::enable_if_t<is_process_sequence_v<TProcess>>> {
+ enum { count = N + std::decay_t<TProcess>::nContinuous };
+ };
+
/**
*
* Definition of a static process list/sequence
@@ -39,14 +63,23 @@ namespace corsika {
* they are just classes. This allows us to handle both, rvalue as
* well as lvalue Processes in the ProcessSequence.
*
- * (The sequence, and the processes use the CRTP, curiously recurring template pattern).
+ * (The sequence, and the processes use the CRTP, curiously recurring template
+ * pattern).
*
+ * Template parameters:
+ * - TProcess1 is of type BaseProcess, either a dedicatd process, or a ProcessSequence
+ * - TProcess2 is of type BaseProcess, either a dedicatd process, or a ProcessSequence
+ * - IndexStart, IndexProcess1, IndexProcess2 are to count and index each
+ *ContinuousProcess in the entire process-chain
**/
- template <typename TProcess1, typename TProcess2 = NullModel,
- int IndexStart = 0,
- int IndexProcess1 = count_continuous<TProcess1>::count,
- int IndexProcess2 = count_continuous<TProcess1, count_continuous<TProcess2>::count>::count>
+ template <typename TProcess1, typename TProcess2 = NullModel, int IndexStart = 0,
+ int IndexProcess1 = count_continuous<
+ TProcess1, count_continuous<TProcess2, IndexStart>::count>::count,
+ int IndexProcess2 = count_continuous<TProcess2, IndexStart>::count>
+ // template <typename TProcess1, typename TProcess2, int IndexStart, int
+ // IndexProcess1,
+ // int IndexProcess2>
class ProcessSequence : public BaseProcess<ProcessSequence<TProcess1, TProcess2>> {
using process1_type = typename std::decay_t<TProcess1>;
@@ -59,15 +92,18 @@ namespace corsika {
static bool constexpr t2SwitchProcSeq = is_switch_process_sequence_v<process2_type>;
// make sure only BaseProcess types TProcess1/2 are passed
- static_assert(std::is_base_of_v<BaseProcess<process1_type>, process1_type>,
+ static_assert(is_base_process_v<process1_type>,
"can only use process derived from BaseProcess in "
"ProcessSequence, for Process 1");
- static_assert(std::is_base_of_v<BaseProcess<process2_type>, process2_type>,
+ static_assert(is_base_process_v<process2_type>,
"can only use process derived from BaseProcess in "
"ProcessSequence, for Process 2");
public:
- enum { nContinuous = IndexProcess2 }; // static counter to index continuous processes
+ /**
+ * static counter to uniquely index (count) all ContinuousProcess in sequence.
+ **/
+ enum { nContinuous = IndexProcess1 };
// resource management
ProcessSequence() = delete; // only initialized objects
@@ -96,7 +132,8 @@ namespace corsika {
typename TParticle::node_type const& to);
template <typename TParticle, typename TTrack>
- inline ProcessReturn doContinuous(TParticle& particle, TTrack& vT);
+ inline ProcessReturn doContinuous(TParticle& particle, TTrack& vT,
+ ContinuousProcessIndex const limitID);
template <typename TSecondaries>
inline void doSecondaries(TSecondaries& vS);
@@ -117,9 +154,19 @@ namespace corsika {
template <typename TStack>
inline void doStack(TStack& stack);
+ /**
+ * Determines the step-length limitation caused by ContinuousProcess-es in this
+ ProcessSequence.
+ *
+ * Returns a ContinuousProcessStepLength object, which contains both: the actual
+ * length (LengthType) as well as the index of the ContinuousProcess inside the
+ * ProcessSequence (for identification).
+
+ **/
+
template <typename TParticle, typename TTrack>
- inline LengthType getMaxStepLength(TParticle& particle, TTrack& vTrack);
- //inline std::pair<LengthType,int> getMaxStepLength(TParticle& particle, TTrack& vTrack);
+ inline ContinuousProcessStepLength getMaxStepLength(TParticle& particle,
+ TTrack& vTrack);
template <typename TParticle>
inline GrammageType getInteractionLength(TParticle&& particle) {
@@ -152,7 +199,6 @@ namespace corsika {
private:
TProcess1 A_; /// process/list A, this is a reference, if possible
TProcess2 B_; /// process/list B, this is a reference, if possible
-
};
/**
@@ -176,12 +222,12 @@ namespace corsika {
*
* \param vA needs to derive from BaseProcess or ProcessSequence
* \param vB paramter-pack, needs to derive BaseProcess or ProcessSequence
+ *
**/
template <typename... TProcesses, typename TProcess1>
inline typename std::enable_if_t<
- std::is_base_of_v<BaseProcess<typename std::decay_t<TProcess1>>,
- typename std::decay_t<TProcess1>>,
+ is_base_process_v<typename std::decay_t<TProcess1>>,
ProcessSequence<TProcess1, decltype(make_sequence(std::declval<TProcesses>()...))>>
make_sequence(TProcess1&& vA, TProcesses&&... vBs) {
return ProcessSequence<TProcess1,
@@ -198,13 +244,11 @@ namespace corsika {
* \param vB needs to derive BaseProcess or ProcessSequence
**/
template <typename TProcess1, typename TProcess2>
- inline typename std::enable_if_t<
- std::is_base_of_v<BaseProcess<typename std::decay_t<TProcess1>>,
- typename std::decay_t<TProcess1>> &&
- std::is_base_of_v<BaseProcess<typename std::decay_t<TProcess2>>,
- typename std::decay_t<TProcess2>>,
- ProcessSequence<TProcess1, TProcess2>>
- make_sequence(TProcess1&& vA, TProcess2&& vB) {
+ inline
+ typename std::enable_if_t<is_base_process_v<typename std::decay_t<TProcess1>> &&
+ is_base_process_v<typename std::decay_t<TProcess2>>,
+ ProcessSequence<TProcess1, TProcess2>>
+ make_sequence(TProcess1&& vA, TProcess2&& vB) {
return ProcessSequence<TProcess1, TProcess2>(vA, vB);
}
@@ -217,10 +261,8 @@ namespace corsika {
* \param vA needs to derive from BaseProcess or ProcessSequence
**/
template <typename TProcess>
- inline typename std::enable_if_t<
- std::is_base_of_v<BaseProcess<typename std::decay_t<TProcess>>,
- typename std::decay_t<TProcess>>,
- ProcessSequence<TProcess, NullModel>>
+ inline typename std::enable_if_t<is_base_process_v<typename std::decay_t<TProcess>>,
+ ProcessSequence<TProcess, NullModel>>
make_sequence(TProcess&& vA) {
return ProcessSequence<TProcess, NullModel>(vA, NullModel());
}
@@ -232,15 +274,10 @@ namespace corsika {
struct is_process_sequence<ProcessSequence<TProcess1, TProcess2>> : std::true_type {
// only switch on for BaseProcesses
template <typename std::enable_if_t<
- std::is_base_of_v<BaseProcess<typename std::decay_t<TProcess1>>,
- typename std::decay_t<TProcess1>> &&
- std::is_base_of_v<BaseProcess<typename std::decay_t<TProcess2>>,
- typename std::decay_t<TProcess2>>,
- int>>
+ is_base_process_v<TProcess1> && is_base_process_v<TProcess2>, int>>
is_process_sequence() {}
};
-
} // namespace corsika
#include <corsika/detail/framework/process/ProcessSequence.inl>
diff --git a/corsika/framework/process/ProcessTraits.hpp b/corsika/framework/process/ProcessTraits.hpp
index 0eb5e7aae..104a40711 100644
--- a/corsika/framework/process/ProcessTraits.hpp
+++ b/corsika/framework/process/ProcessTraits.hpp
@@ -12,12 +12,40 @@
* \file ProcessTraits.hpp
*/
-#include <corsika/framework/process/ProcessTraits.hpp>
+//#include <corsika/framework/process/BaseProcess.hpp>
+//#include <corsika/framework/process/ProcessSequence.hpp>
+//#include <corsika/framework/process/SwitchProcessSequence.hpp>
+//#include <corsika/framework/process/ContinuousProcess.hpp>
#include <type_traits>
namespace corsika {
+ /**
+ * A traits marker to identify BaseProcess
+ */
+ template <typename TProcess, typename Enable = void>
+ struct is_base_process : std::false_type {};
+
+ template <typename TProcess>
+ bool constexpr is_base_process_v = is_base_process<TProcess>::value;
+
+ template <typename TProcess>
+ struct is_base_process<
+ TProcess,
+ std::enable_if_t<std::is_base_of_v<BaseProcess<typename std::decay_t<TProcess>>,
+ typename std::decay_t<TProcess>>>>
+ : std::true_type {};
+
+ /**
+ * A traits marker to identify ContinuousProcess
+ */
+ template <typename TProcess, typename Enable = void>
+ struct is_continuous_process : std::false_type {};
+
+ template <typename TProcess>
+ bool constexpr is_continuous_process_v = is_continuous_process<TProcess>::value;
+
/**
* A traits marker to track which BaseProcess is also a ProcessSequence
**/
@@ -47,34 +75,11 @@ namespace corsika {
bool constexpr contains_stack_process_v = contains_stack_process<TClass>::value;
/**
- * traits class to count ContinuousProcess-es
- **/
- template <typename TProcess, int N=0, typename Enable=void>
- struct count_continuous {
- enum { count = N };
- };
-
-
- /**
- * traits class to count ContinuousProcess-es
- **/
- template <typename TProcessSequence, int N>
- struct count_continuous<TProcessSequence, N,
- typename std::enable_if_t<is_process_sequence_v<TProcessSequence>>> {
- enum { count = N+TProcessSequence::nContinuous };
- };
-
- /**
- * traits class to count ContinuousProcess-es
+ * traits class to count ContinuousProcess-es, general version
**/
-
- template <typename TSwitchProcessSequence, int N>
- struct count_continuous<TSwitchProcessSequence, N,
- typename std::enable_if_t<is_switch_process_sequence_v<TSwitchProcessSequence>>> {
- enum { count = N+TSwitchProcessSequence::nContinuous };
- };
-
-
-
+ template <typename TProcess, int N = 0, typename Enable = void>
+ struct count_continuous {
+ enum { count = N };
+ };
} // namespace corsika
diff --git a/corsika/framework/process/SwitchProcessSequence.hpp b/corsika/framework/process/SwitchProcessSequence.hpp
index 95f21a60d..df951154a 100644
--- a/corsika/framework/process/SwitchProcessSequence.hpp
+++ b/corsika/framework/process/SwitchProcessSequence.hpp
@@ -16,6 +16,7 @@
#include <corsika/framework/process/ProcessTraits.hpp>
#include <corsika/framework/process/BoundaryCrossingProcess.hpp>
#include <corsika/framework/process/ContinuousProcess.hpp>
+#include <corsika/framework/process/ContinuousProcessIndex.hpp>
#include <corsika/framework/process/DecayProcess.hpp>
#include <corsika/framework/process/InteractionProcess.hpp>
#include <corsika/framework/process/ProcessReturn.hpp>
@@ -29,6 +30,32 @@
namespace corsika {
+ /*
+ template <typename TProcess1, typename TProcess2, typename TSelect, int IndexStart = 0,
+ int IndexProcess1 = count_continuous<TProcess1>::count,
+ int IndexProcess2 = count_continuous<TProcess2>::count>
+ class SwitchProcessSequence;
+ */
+
+ /**
+ * traits class to count ContinuousProcess-es, specialized for ProcessSequence-es
+ **/
+
+ /* template <typename TProcess1, typename TProcess2, typename TSelect, int N>
+ struct count_continuous<SwitchProcessSequence<TProcess1, TProcess2, TSelect>, N> {
+ enum {
+ count = N + SwitchProcessSequence<TProcess1, TProcess2, TSelect>::nContinuous
+ };
+ };*/
+
+ /*
+ template <typename TSProcess, int N>
+ struct count_continuous<
+ TSProcess, N, typename std::enable_if_t<is_switch_process_sequence_v<TSProcess>>> {
+ enum { count = N + TSProcess::nContinuous };
+ };
+ */
+
/**
* enum for the process switch selection: identify if First or
* Second process branch should be used.
@@ -57,13 +84,19 @@ namespace corsika {
since this makes no sense. The StackProcess acts on an entire
particle stack and not on indiviidual particles.
+ Template parameters:
+ - TProcess1 is of type BaseProcess, either a dedicatd process, or a ProcessSequence
+ - TProcess2 is of type BaseProcess, either a dedicatd process, or a ProcessSequence
+ - IndexStart, IndexProcess1, IndexProcess2 are to count and index each
+ ContinuousProcess in the entire process-chain
+
+
See also class \sa ProcessSequence
**/
- template <typename TProcess1, typename TProcess2, typename TSelect,
- int IndexStart = 0,
- int IndexProcess1 = count_continuous<TProcess1>::count,
- int IndexProcess2 = count_continuous<TProcess1, count_continuous<TProcess2>::count>::count>
+ template <typename TProcess1, typename TProcess2, typename TSelect, int IndexStart = 0,
+ int IndexProcess1 = count_continuous<TProcess1, IndexStart>::count,
+ int IndexProcess2 = count_continuous<TProcess2, IndexProcess1>::count>
class SwitchProcessSequence
: public BaseProcess<SwitchProcessSequence<TProcess1, TProcess2, TSelect>> {
@@ -74,10 +107,10 @@ namespace corsika {
static bool constexpr t2ProcSeq = is_process_sequence_v<process2_type>;
// make sure only BaseProcess types TProcess1/2 are passed
- static_assert(std::is_base_of_v<BaseProcess<process1_type>, process1_type>,
+ static_assert(is_base_process_v<process1_type>,
"can only use process derived from BaseProcess in "
"SwitchProcessSequence, for Process 1");
- static_assert(std::is_base_of_v<BaseProcess<process2_type>, process2_type>,
+ static_assert(is_base_process_v<process2_type>,
"can only use process derived from BaseProcess in "
"SwitchProcessSequence, for Process 2");
@@ -97,7 +130,10 @@ namespace corsika {
"ProcessSequence 2");
public:
- enum { nContinuous = IndexProcess1+IndexProcess2 }; // static counter to index continuous processes
+ /**
+ * static counter to uniquely index (count) all ContinuousProcess in switch sequence.
+ **/
+ enum { nContinuous = IndexProcess2 };
// resource management
SwitchProcessSequence() = delete; // only initialized objects
@@ -127,13 +163,15 @@ namespace corsika {
TVTNType const& to);
template <typename TParticle, typename TTrack>
- inline ProcessReturn doContinuous(TParticle& particle, TTrack& vT);
+ inline ProcessReturn doContinuous(TParticle& particle, TTrack& vT,
+ ContinuousProcessIndex const limitId);
template <typename TSecondaries>
inline void doSecondaries(TSecondaries& vS);
template <typename TParticle, typename TTrack>
- inline LengthType getMaxStepLength(TParticle& particle, TTrack& vTrack);
+ inline ContinuousProcessStepLength getMaxStepLength(TParticle& particle,
+ TTrack& vTrack);
template <typename TParticle>
inline GrammageType getInteractionLength(TParticle&& particle) {
@@ -183,17 +221,14 @@ namespace corsika {
**/
template <typename TProcess1, typename TProcess2, typename TSelect>
- inline typename std::enable_if_t<
- std::is_base_of_v<BaseProcess<typename std::decay_t<TProcess1>>,
- typename std::decay_t<TProcess1>> &&
- std::is_base_of_v<BaseProcess<typename std::decay_t<TProcess2>>,
- typename std::decay_t<TProcess2>>,
- SwitchProcessSequence<TProcess1, TProcess2, TSelect>>
- make_select(TProcess1&& vA, TProcess2&& vB, TSelect selector) {
+ inline
+ typename std::enable_if_t<is_base_process_v<typename std::decay_t<TProcess1>> &&
+ is_base_process_v<typename std::decay_t<TProcess2>>,
+ SwitchProcessSequence<TProcess1, TProcess2, TSelect>>
+ make_select(TProcess1&& vA, TProcess2&& vB, TSelect selector) {
return SwitchProcessSequence<TProcess1, TProcess2, TSelect>(vA, vB, selector);
}
-
} // namespace corsika
#include <corsika/detail/framework/process/SwitchProcessSequence.inl>
diff --git a/corsika/modules/LongitudinalProfile.hpp b/corsika/modules/LongitudinalProfile.hpp
index baaa50468..b06e2ef85 100644
--- a/corsika/modules/LongitudinalProfile.hpp
+++ b/corsika/modules/LongitudinalProfile.hpp
@@ -41,7 +41,9 @@ namespace corsika {
GrammageType dX = 10_g / square(1_cm)); // profile binning);
template <typename TParticle, typename TTrack>
- ProcessReturn doContinuous(TParticle const&, TTrack const&);
+ ProcessReturn doContinuous(
+ TParticle const&, TTrack const&,
+ bool const flagLimit = false); // not needed for LongitudinalProfile
template <typename TParticle, typename TTrack>
LengthType getMaxStepLength(TParticle const&, TTrack const&) {
diff --git a/corsika/modules/ParticleCut.hpp b/corsika/modules/ParticleCut.hpp
index 1678c528e..9842dd034 100644
--- a/corsika/modules/ParticleCut.hpp
+++ b/corsika/modules/ParticleCut.hpp
@@ -50,8 +50,10 @@ namespace corsika {
bool const em, bool const inv);
void doSecondaries(corsika::setup::StackView&);
- ProcessReturn doContinuous(corsika::setup::Stack::particle_type& vParticle,
- corsika::setup::Trajectory const& vTrajectory);
+ ProcessReturn doContinuous(
+ corsika::setup::Stack::particle_type& vParticle,
+ corsika::setup::Trajectory const& vTrajectory,
+ const bool limitFlag = false); // this is not used for ParticleCut
LengthType getMaxStepLength(corsika::setup::Stack::particle_type const&,
corsika::setup::Trajectory const&) {
return meter * std::numeric_limits<double>::infinity();
diff --git a/corsika/modules/TrackWriter.hpp b/corsika/modules/TrackWriter.hpp
index c437bc7d1..5bf9856b3 100644
--- a/corsika/modules/TrackWriter.hpp
+++ b/corsika/modules/TrackWriter.hpp
@@ -22,7 +22,7 @@ namespace corsika {
TrackWriter(std::string const& filename);
template <typename TParticle, typename TTrack>
- ProcessReturn doContinuous(TParticle const&, TTrack const&);
+ ProcessReturn doContinuous(TParticle const&, TTrack const&, bool const limitFlag);
template <typename TParticle, typename TTrack>
LengthType getMaxStepLength(TParticle const&, TTrack const&);
diff --git a/corsika/modules/energy_loss/BetheBlochPDG.hpp b/corsika/modules/energy_loss/BetheBlochPDG.hpp
index f8e875937..582a0c4e0 100644
--- a/corsika/modules/energy_loss/BetheBlochPDG.hpp
+++ b/corsika/modules/energy_loss/BetheBlochPDG.hpp
@@ -43,7 +43,17 @@ namespace corsika {
public:
BetheBlochPDG(ShowerAxis const& showerAxis);
- ProcessReturn doContinuous(setup::Stack::particle_type&, setup::Trajectory const&);
+ /** clang-format-off
+ * Interface function of ContinuousProcess.
+ *
+ * \param particle The particle to process in its current state
+ * \param track The trajectory in space of this particle, on which doContinuous should
+ * act
+ * \param limitFlag flag to identify, if BetheBlochPDG::getMaxStepLength is the
+ * globally limiting factor (or not)
+ clang-format-on **/
+ ProcessReturn doContinuous(setup::Stack::particle_type& particle,
+ setup::Trajectory const& track, bool const limitFlag);
LengthType getMaxStepLength(setup::Stack::particle_type const&,
setup::Trajectory const&)
const; //! limited by the energy threshold! By default the limit is the particle
diff --git a/corsika/modules/proposal/ContinuousProcess.hpp b/corsika/modules/proposal/ContinuousProcess.hpp
index 3b8976f3a..0efbf12ec 100644
--- a/corsika/modules/proposal/ContinuousProcess.hpp
+++ b/corsika/modules/proposal/ContinuousProcess.hpp
@@ -67,8 +67,11 @@ namespace corsika::proposal {
//! If the particle if below the given energy threshold where it will be
//! considered stochastically, it will be absorbed.
//!
+ //! \param limitFlag is true, if the track was actually limited by
+ //! proposal::ContinuousProcess::getMaxStepLength
+ //!
template <typename TParticle, typename TTrack>
- ProcessReturn doContinuous(TParticle&, TTrack const&);
+ ProcessReturn doContinuous(TParticle&, TTrack const& track, bool const limitFlag);
//!
//! Calculates maximal step length of process.
diff --git a/examples/staticsequence_example.cpp b/examples/staticsequence_example.cpp
index d7b9e06b5..3f6c5aff7 100644
--- a/examples/staticsequence_example.cpp
+++ b/examples/staticsequence_example.cpp
@@ -24,7 +24,7 @@ class Process1 : public ContinuousProcess<Process1> {
public:
Process1() {}
template <typename D, typename T>
- ProcessReturn doContinuous(D& d, T&) const {
+ ProcessReturn doContinuous(D& d, T&, bool const) const {
for (int i = 0; i < nData; ++i) d.p[i] += 1;
return ProcessReturn::Ok;
}
@@ -35,7 +35,7 @@ public:
Process2() {}
template <typename D, typename T>
- inline ProcessReturn doContinuous(D& d, T&) const {
+ inline ProcessReturn doContinuous(D& d, T&, bool const) const {
for (int i = 0; i < nData; ++i) d.p[i] -= 0.1 * i;
return ProcessReturn::Ok;
}
@@ -46,7 +46,7 @@ public:
Process3() {}
template <typename D, typename T>
- inline ProcessReturn doContinuous(D&, T&) const {
+ inline ProcessReturn doContinuous(D&, T&, bool const) const {
return ProcessReturn::Ok;
}
};
@@ -56,7 +56,7 @@ public:
Process4(const double v)
: fV(v) {}
template <typename D, typename T>
- inline ProcessReturn doContinuous(D& d, T&) const {
+ inline ProcessReturn doContinuous(D& d, T&, bool const) const {
for (int i = 0; i < nData; ++i) d.p[i] *= fV;
return ProcessReturn::Ok;
}
@@ -90,7 +90,7 @@ void modular() {
const int nEv = 10;
for (int iEv = 0; iEv < nEv; ++iEv) {
- sequence.doContinuous(particle, track);
+ sequence.doContinuous(particle, track, ContinuousProcessIndex(0));
for (int i = 0; i < nData; ++i) {
check[i] += 1. - 0.1 * i;
check[i] *= 1.5;
diff --git a/examples/vertical_EAS.cpp b/examples/vertical_EAS.cpp
index 4ebffab2e..a661ca3ed 100644
--- a/examples/vertical_EAS.cpp
+++ b/examples/vertical_EAS.cpp
@@ -254,9 +254,9 @@ int main(int argc, char** argv) {
auto hadronSequence = make_select(
urqmdCounted, make_sequence(sibyllNucCounted, sibyllCounted), EnergySwitch(55_GeV));
auto decaySequence = make_sequence(decayPythia, decaySibyll);
- auto sequence = make_sequence(stackInspect, hadronSequence, reset_particle_mass,
- decaySequence, em_continuous, cut,
- trackWriter, observationLevel, longprof);
+ auto sequence =
+ make_sequence(stackInspect, hadronSequence, reset_particle_mass, decaySequence,
+ em_continuous, cut, trackWriter, observationLevel, longprof);
// define air shower object, run simulation
setup::Tracking tracking;
diff --git a/modules/conex b/modules/conex
index 73e3d1fa8..180215204 160000
--- a/modules/conex
+++ b/modules/conex
@@ -1 +1 @@
-Subproject commit 73e3d1fa850d2d3602a15b3b948ac1789fbc968d
+Subproject commit 180215204035fa17c5f6cacb42588dca85d171ca
diff --git a/modules/data b/modules/data
index afe83dc19..fb7577314 160000
--- a/modules/data
+++ b/modules/data
@@ -1 +1 @@
-Subproject commit afe83dc19c1464deb176f38c3ddab80a64e081f4
+Subproject commit fb7577314e5e3c837fa5a3006dc4f3c0146cabf7
diff --git a/tests/framework/testProcessSequence.cpp b/tests/framework/testProcessSequence.cpp
index 4dd0efc8b..5ae3239e4 100644
--- a/tests/framework/testProcessSequence.cpp
+++ b/tests/framework/testProcessSequence.cpp
@@ -10,6 +10,7 @@
#include <corsika/framework/process/SwitchProcessSequence.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/process/ProcessTraits.hpp>
+#include <corsika/framework/process/ContinuousProcessStepLength.hpp>
#include <catch2/catch.hpp>
@@ -18,10 +19,12 @@
#include <iostream>
#include <typeinfo>
+#include <boost/type_index.hpp>
+
using namespace corsika;
using namespace std;
-static const int nData = 10;
+static int const nData = 10;
int globalCount = 0; // simple counter
@@ -32,68 +35,110 @@ int checkCont = 0; // use this as a bit field
class ContinuousProcess1 : public ContinuousProcess<ContinuousProcess1> {
public:
- ContinuousProcess1(const int v)
- : v_(v) {
+ ContinuousProcess1(int const v, LengthType const step)
+ : v_(v)
+ , step_(step) {
cout << "globalCount: " << globalCount << ", v_: " << v_ << std::endl;
globalCount++;
}
+ void setStep(LengthType const v) { step_ = v; }
+
template <typename D, typename T>
- inline ProcessReturn doContinuous(D& d, T&) const {
+ inline ProcessReturn doContinuous(D& d, T&, bool const flag) const {
+ flag_ = flag;
cout << "ContinuousProcess1::DoContinuous" << endl;
checkCont |= 1;
for (int i = 0; i < nData; ++i) d.data_[i] += 0.933;
return ProcessReturn::Ok;
}
+ template <typename TParticle, typename TTrack>
+ inline LengthType getMaxStepLength(TParticle&, TTrack&) {
+ return step_;
+ }
+
+ bool getFlag() const { return flag_; }
+ void resetFlag() { flag_ = false; }
+
private:
int v_ = 0;
+ LengthType step_ = 0_m;
+ mutable bool flag_ = false;
};
class ContinuousProcess2 : public ContinuousProcess<ContinuousProcess2> {
public:
- ContinuousProcess2(const int v)
- : v_(v) {
+ ContinuousProcess2(int const v, LengthType const step)
+ : v_(v)
+ , step_(step) {
cout << "globalCount: " << globalCount << ", v_: " << v_ << std::endl;
globalCount++;
}
+ void setStep(LengthType const v) { step_ = v; }
+
template <typename D, typename T>
- inline ProcessReturn doContinuous(D& d, T&) const {
+ inline ProcessReturn doContinuous(D& d, T&, bool const flag) const {
+ flag_ = flag;
cout << "ContinuousProcess2::DoContinuous" << endl;
checkCont |= 2;
for (int i = 0; i < nData; ++i) d.data_[i] += 0.111;
return ProcessReturn::Ok;
}
+ template <typename TParticle, typename TTrack>
+ inline LengthType getMaxStepLength(TParticle&, TTrack&) {
+ return step_;
+ }
+
+ bool getFlag() const { return flag_; }
+ void resetFlag() { flag_ = false; }
+
private:
int v_ = 0;
+ LengthType step_ = 0_m;
+ mutable bool flag_ = false;
};
class ContinuousProcess3 : public ContinuousProcess<ContinuousProcess3> {
public:
- ContinuousProcess3(const int v)
- : v_(v) {
+ ContinuousProcess3(int const v, LengthType const step)
+ : v_(v)
+ , step_(step) {
cout << "globalCount: " << globalCount << ", v_: " << v_ << std::endl;
globalCount++;
}
+ void setStep(LengthType const v) { step_ = v; }
+
template <typename D, typename T>
- inline ProcessReturn doContinuous(D& d, T&) const {
+ inline ProcessReturn doContinuous(D& d, T&, bool const flag) const {
+ flag_ = flag;
cout << "ContinuousProcess3::DoContinuous" << endl;
checkCont |= 4;
for (int i = 0; i < nData; ++i) d.data_[i] += 0.333;
return ProcessReturn::Ok;
}
+ template <typename TParticle, typename TTrack>
+ inline LengthType getMaxStepLength(TParticle&, TTrack&) {
+ return step_;
+ }
+
+ bool getFlag() const { return flag_; }
+ void resetFlag() { flag_ = false; }
+
private:
int v_ = 0;
+ LengthType step_ = 0_m;
+ mutable bool flag_ = false;
};
class Process1 : public InteractionProcess<Process1> {
public:
- Process1(const int v)
+ Process1(int const v)
: v_(v) {
cout << "globalCount: " << globalCount << ", v_: " << v_ << std::endl;
globalCount++;
@@ -116,7 +161,7 @@ private:
class Process2 : public InteractionProcess<Process2> {
public:
- Process2(const int v)
+ Process2(int const v)
: v_(v) {
cout << "globalCount: " << globalCount << ", v_: " << v_ << std::endl;
globalCount++;
@@ -140,7 +185,7 @@ private:
class Process3 : public InteractionProcess<Process3> {
public:
- Process3(const int v)
+ Process3(int const v)
: v_(v) {
cout << "globalCount: " << globalCount << ", v_: " << v_ << std::endl;
globalCount++;
@@ -164,14 +209,14 @@ private:
class Process4 : public BaseProcess<Process4> {
public:
- Process4(const int v)
+ Process4(int const v)
: v_(v) {
cout << "globalCount: " << globalCount << ", v_: " << v_ << std::endl;
globalCount++;
}
template <typename D, typename T>
- inline ProcessReturn doContinuous(D& d, T&) const {
+ inline ProcessReturn doContinuous(D& d, T&, bool const) const {
std::cout << "Base::doContinuous" << std::endl;
checkCont |= 8;
for (int i = 0; i < nData; ++i) { d.data_[i] /= 1.2; }
@@ -188,7 +233,7 @@ private:
class Decay1 : public DecayProcess<Decay1> {
public:
- Decay1(const int) {
+ Decay1(int const) {
cout << "Decay1()" << endl;
globalCount++;
}
@@ -205,7 +250,7 @@ public:
class Decay2 : public DecayProcess<Decay2> {
public:
- Decay2(const int) {
+ Decay2(int const) {
cout << "Decay2()" << endl;
globalCount++;
}
@@ -222,7 +267,7 @@ public:
class Stack1 : public StackProcess<Stack1> {
public:
- Stack1(const int n)
+ Stack1(int const n)
: StackProcess(n) {}
template <typename TStack>
ProcessReturn doStack(TStack&) {
@@ -235,12 +280,16 @@ private:
int count_ = 0;
};
+// The stack is non-existent for this example
struct DummyStack {};
+// our data object (particle) is a simple arrary of doubles
struct DummyData {
double data_[nData] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
};
+// there is no real trajectory/track
struct DummyTrajectory {};
-
+// since there is no stack, there is also no view. This is a simplistic dummy object
+// sufficient here.
struct DummyView {
DummyView(DummyData& p)
: p_(p) {}
@@ -253,6 +302,18 @@ TEST_CASE("Process Sequence General", "ProcessSequence") {
logging::set_level(logging::level::info);
corsika_logger->set_pattern("[%n:%^%-8l%$]: %v");
+ SECTION("BaseProcess") {
+
+ Process1 m1(0);
+ const Process4 m4(3);
+
+ CHECK(is_base_process_v<Process1>);
+ CHECK_FALSE(is_base_process_v<DummyData>);
+ CHECK(is_base_process_v<decltype(m4)>);
+ CHECK(is_base_process_v<decltype(Decay1(1))>);
+ CHECK(is_base_process_v<decltype(ContinuousProcess3{3, 3_m})>);
+ }
+
SECTION("Check construction") {
globalCount = 0;
Process1 m1(0);
@@ -265,10 +326,15 @@ TEST_CASE("Process Sequence General", "ProcessSequence") {
CHECK(globalCount == 4);
auto sequence1 = make_sequence(m1, m2, m3, m4);
- CHECK(is_process_sequence_v<decltype(sequence1)> == true);
- CHECK(is_process_sequence_v<decltype(m2)> == false);
- CHECK(is_switch_process_sequence_v<decltype(sequence1)> == false);
- CHECK(is_switch_process_sequence_v<decltype(m2)> == false);
+ CHECK(is_base_process_v<decltype(sequence1)>);
+ CHECK(is_base_process_v<decltype(m2)>);
+ CHECK(is_process_sequence_v<decltype(sequence1)>);
+ CHECK_FALSE(is_process_sequence_v<decltype(m2)>);
+ CHECK_FALSE(is_switch_process_sequence_v<decltype(sequence1)>);
+ CHECK_FALSE(is_switch_process_sequence_v<decltype(m2)>);
+
+ CHECK_FALSE(is_process_sequence_v<decltype(Decay1(7))>);
+ CHECK_FALSE(is_switch_process_sequence_v<decltype(Decay1(7))>);
auto sequence2 = make_sequence(m1, m2, m3);
CHECK(is_process_sequence_v<decltype(sequence2)> == true);
@@ -279,7 +345,7 @@ TEST_CASE("Process Sequence General", "ProcessSequence") {
SECTION("interaction length") {
globalCount = 0;
- ContinuousProcess1 cp1(0);
+ ContinuousProcess1 cp1(0, 1_m);
Process2 m2(1);
Process3 m3(2);
@@ -297,7 +363,7 @@ TEST_CASE("Process Sequence General", "ProcessSequence") {
SECTION("lifetime") {
globalCount = 0;
- ContinuousProcess1 cp1(0);
+ ContinuousProcess1 cp1(0, 1_m);
Process2 m2(1);
Process3 m3(2);
Decay1 d3(3);
@@ -314,30 +380,61 @@ TEST_CASE("Process Sequence General", "ProcessSequence") {
globalCount = 0;
}
- SECTION("sectionTwo") {
+ SECTION("ContinousProcess") {
globalCount = 0;
- ContinuousProcess1 cp1(0); // += 0.933
- ContinuousProcess2 cp2(1); // += 0.111
- Process2 m2(2); // /= 1.1
- Process3 m3(3); // *= 1.01
+ ContinuousProcess1 cp1(0, 1_m); // += 0.933
+ ContinuousProcess2 cp2(1, 1.1_m); // += 0.111
+ Process2 m2(2); // /= 1.1
+ Process3 m3(3); // *= 1.01
auto sequence2 = make_sequence(cp1, m2, m3, cp2);
+ std::cout << boost::typeindex::type_id<decltype(sequence2)>().pretty_name()
+ << std::endl;
+
DummyData particle;
DummyTrajectory track;
+ cp1.resetFlag();
+ cp2.resetFlag();
+
+ ContinuousProcessStepLength const step1 = sequence2.getMaxStepLength(particle, track);
+ CHECK(LengthType(step1) == 1_m);
+ sequence2.doContinuous(particle, track, step1);
+ CHECK(cp1.getFlag());
+ CHECK_FALSE(cp2.getFlag());
+ CORSIKA_LOG_INFO("step1, l={}, i={}", LengthType(step1),
+ ContinuousProcessIndex(step1).getIndex());
+
+ cp1.resetFlag();
+ cp2.resetFlag();
+
+ cp1.setStep(10_m);
+ ContinuousProcessStepLength const step2 = sequence2.getMaxStepLength(particle, track);
+ CHECK(LengthType(step2) == 1.1_m);
+ CHECK(ContinuousProcessIndex(step1) != ContinuousProcessIndex(step2));
+ sequence2.doContinuous(particle, track, step2);
+ CHECK_FALSE(cp1.getFlag());
+ CHECK(cp2.getFlag());
+ CORSIKA_LOG_INFO("step2, l={}, i={}", LengthType(step2),
+ ContinuousProcessIndex(step2).getIndex());
+
cout << "-->init sequence2" << endl;
globalCount = 0;
- cout << "-->docont" << endl;
+ cout << "-->docontinuous" << endl;
// validation data
double test_data[nData] = {0};
- const int nLoop = 5;
+ // reset
+ particle = DummyData();
+ track = DummyTrajectory();
+
+ int const nLoop = 5;
cout << "Running loop with n=" << nLoop << endl;
for (int iLoop = 0; iLoop < nLoop; ++iLoop) {
for (int i = 0; i < nData; ++i) { test_data[i] += 0.933 + 0.111; }
- sequence2.doContinuous(particle, track);
+ sequence2.doContinuous(particle, track, ContinuousProcessIndex(1));
}
for (int i = 0; i < nData; i++) {
cout << "data_[" << i << "]=" << particle.data_[i] << endl;
@@ -356,14 +453,14 @@ TEST_CASE("Process Sequence General", "ProcessSequence") {
DummyStack stack;
- const int nLoop = 20;
+ int const nLoop = 20;
for (int i = 0; i < nLoop; ++i) { sequence1.doStack(stack); }
CHECK(s1.getCount() == 20);
CHECK(s2.getCount() == 10);
- ContinuousProcess2 cp2(1); // += 0.111
- Process2 m2(2); // /= 1.1
+ ContinuousProcess2 cp2(1, 2_m); // += 0.111
+ Process2 m2(2); // /= 1.1
auto sequence2 = make_sequence(cp2, m2);
auto sequence3 = make_sequence(cp2, m2, s1);
@@ -379,34 +476,59 @@ TEST_CASE("Switch Process Sequence", "ProcessSequence") {
logging::set_level(logging::level::info);
corsika_logger->set_pattern("[%n:%^%-8l%$]: %v");
+ /**
+ * In this example switching is done only by "data_[0]>0", where
+ * data in an arrray of doubles, DummyData.
+ */
+
+ struct SwitchSelect {
+ SwitchResult operator()(DummyData const& p) const {
+ if (p.data_[0] > 0) return SwitchResult::First;
+ return SwitchResult::Second;
+ }
+ };
+ SwitchSelect select1;
+
+ auto cp1 = ContinuousProcess1(0, 1_m);
+ auto cp2 = ContinuousProcess2(0, 2_m);
+ auto cp3 = ContinuousProcess3(0, 3_m);
+
+ auto sequence1 = make_sequence(Process1(0), cp2, Decay1(0));
+ auto sequence2 = make_sequence(cp3, Process2(0), Decay2(0));
+
+ auto sequence3 = make_sequence(cp1, Process3(0),
+ SwitchProcessSequence(sequence1, sequence2, select1));
+
SECTION("Check construction") {
- struct TestSelect {
- SwitchResult operator()(const DummyData& p) const {
- std::cout << "TestSelect data=" << p.data_[0] << std::endl;
- if (p.data_[0] > 0) return SwitchResult::First;
- return SwitchResult::Second;
- }
- };
- TestSelect select1;
+ auto sequence_alt =
+ make_sequence(cp1, Process3(0),
+ make_select(make_sequence(Process1(0), cp2, Decay1(0)),
+ make_sequence(cp3, Process2(0), Decay2(0)), select1));
- auto sequence1 = make_sequence(Process1(0), ContinuousProcess2(0), Decay1(0));
- auto sequence2 = make_sequence(ContinuousProcess3(0), Process2(0), Decay2(0));
+ auto switch_seq = SwitchProcessSequence(sequence1, sequence2, select1);
+ CHECK(is_process_sequence_v<decltype(switch_seq)>);
+ CHECK(is_switch_process_sequence_v<decltype(switch_seq)>);
+ // CHECK(is_switch_process_sequence_v<decltype(&switch_seq)>);
+ CHECK(is_switch_process_sequence_v<decltype(
+ SwitchProcessSequence(sequence1, sequence2, select1))>);
- auto sequence3 = make_sequence(ContinuousProcess1(0), Process3(0),
- SwitchProcessSequence(sequence1, sequence2, select1));
+ CHECK(is_process_sequence_v<decltype(sequence3)>);
+ CHECK_FALSE(is_switch_process_sequence_v<decltype(sequence3)>);
- auto sequence_alt = make_sequence(
- ContinuousProcess1(0), Process3(0),
- make_select(make_sequence(Process1(0), ContinuousProcess2(0), Decay1(0)),
- make_sequence(ContinuousProcess3(0), Process2(0), Decay2(0)),
- select1));
+ CHECK(is_process_sequence_v<decltype(
+ SwitchProcessSequence(sequence1, sequence2, select1))>);
+ CHECK(is_switch_process_sequence_v<decltype(
+ SwitchProcessSequence(sequence1, sequence2, select1))>);
// check that same process sequence can be build in different ways
CHECK(typeid(sequence3) == typeid(sequence_alt));
- CHECK(is_process_sequence_v<decltype(sequence3)> == true);
+ CHECK(is_process_sequence_v<decltype(sequence3)>);
CHECK(is_process_sequence_v<decltype(
- SwitchProcessSequence(sequence1, sequence2, select1))> == true);
+ SwitchProcessSequence(sequence1, sequence2, select1))>);
+ }
+
+ SECTION("Check interfaces") {
DummyData particle;
DummyTrajectory track;
@@ -417,7 +539,7 @@ TEST_CASE("Switch Process Sequence", "ProcessSequence") {
checkSec = 0;
checkCont = 0;
particle.data_[0] = 100; // data positive
- sequence3.doContinuous(particle, track);
+ sequence3.doContinuous(particle, track, ContinuousProcessIndex(1));
CHECK(checkInteract == 0);
CHECK(checkDecay == 0);
CHECK(checkCont == 0b011);
@@ -428,7 +550,7 @@ TEST_CASE("Switch Process Sequence", "ProcessSequence") {
checkSec = 0;
checkCont = 0;
particle.data_[0] = -100; // data negative
- sequence_alt.doContinuous(particle, track);
+ sequence3.doContinuous(particle, track, ContinuousProcessIndex(1));
CHECK(checkInteract == 0);
CHECK(checkDecay == 0);
CHECK(checkCont == 0b101);
@@ -479,6 +601,95 @@ TEST_CASE("Switch Process Sequence", "ProcessSequence") {
CHECK(checkCont == 0);
CHECK(checkSec == 0);
}
+
+ SECTION("Check ContinuousProcesses in SwitchProcessSequence") {
+
+ DummyData particle;
+ DummyTrajectory track;
+
+ particle.data_[0] =
+ 100; // data positive, selects particular branch on SwitchProcessSequence
+
+ cp1.setStep(10_m);
+ cp2.setStep(15_m);
+ cp3.setStep(100_m);
+
+ cp1.resetFlag();
+ cp2.resetFlag();
+ cp3.resetFlag();
+
+ ContinuousProcessStepLength const step1 = sequence3.getMaxStepLength(particle, track);
+ CHECK(LengthType(step1) == 10_m);
+ sequence3.doContinuous(particle, track, step1);
+ CHECK(cp1.getFlag());
+ CHECK_FALSE(cp2.getFlag());
+ CHECK_FALSE(cp3.getFlag());
+ CORSIKA_LOG_INFO("step1, l={}, i={}", LengthType(step1),
+ ContinuousProcessIndex(step1).getIndex());
+
+ particle.data_[0] =
+ 100; // data positive, selects particular branch on SwitchProcessSequence
+
+ cp1.setStep(50_m);
+ cp2.setStep(15_m);
+ cp3.setStep(100_m);
+
+ cp1.resetFlag();
+ cp2.resetFlag();
+ cp3.resetFlag();
+
+ ContinuousProcessStepLength const step2 = sequence3.getMaxStepLength(particle, track);
+ CHECK(LengthType(step2) == 15_m);
+ sequence3.doContinuous(particle, track, step2);
+ CHECK_FALSE(cp1.getFlag());
+ CHECK(cp2.getFlag());
+ CHECK_FALSE(cp3.getFlag());
+ CORSIKA_LOG_INFO("step2, len_cont={}, indexLimit={} type={}", LengthType(step2),
+ ContinuousProcessIndex(step2).getIndex(),
+ boost::typeindex::type_id<decltype(sequence3)>().pretty_name());
+
+ particle.data_[0] =
+ -100; // data positive, selects particular branch on SwitchProcessSequence
+
+ cp1.setStep(11_m);
+ cp2.setStep(15_m);
+ cp3.setStep(100_m);
+
+ cp1.resetFlag();
+ cp2.resetFlag();
+ cp3.resetFlag();
+
+ ContinuousProcessStepLength const step3 = sequence3.getMaxStepLength(particle, track);
+ CHECK(LengthType(step3) == 11_m);
+ sequence3.doContinuous(particle, track, step3);
+ CHECK(cp1.getFlag());
+ CHECK_FALSE(cp2.getFlag());
+ CHECK_FALSE(cp3.getFlag());
+ CORSIKA_LOG_INFO("step3, len_cont={}, indexLimit={} type={}", LengthType(step3),
+ ContinuousProcessIndex(step3).getIndex(),
+ boost::typeindex::type_id<decltype(sequence3)>().pretty_name());
+
+ particle.data_[0] =
+ -100; // data positive, selects particular branch on SwitchProcessSequence
+
+ cp1.setStep(11_m);
+ cp2.setStep(15_m);
+ cp3.setStep(2_m);
+
+ cp1.resetFlag();
+ cp2.resetFlag();
+ cp3.resetFlag();
+
+ ContinuousProcessStepLength const step4 = sequence3.getMaxStepLength(particle, track);
+ CHECK(LengthType(step4) == 2_m);
+ sequence3.doContinuous(particle, track, step4);
+ CHECK_FALSE(cp1.getFlag());
+ CHECK_FALSE(cp2.getFlag());
+ CHECK(cp3.getFlag());
+ CORSIKA_LOG_INFO("step4, len_cont={}, indexLimit={} type={}", LengthType(step4),
+ ContinuousProcessIndex(step4).getIndex(),
+ boost::typeindex::type_id<decltype(sequence3)>().pretty_name());
+ }
}
TEST_CASE("Continuous Process Indexing", "ProcessSequence") {
@@ -486,16 +697,24 @@ TEST_CASE("Continuous Process Indexing", "ProcessSequence") {
logging::set_level(logging::level::info);
corsika_logger->set_pattern("[%n:%^%-8l%$]: %v");
- SECTION("Counting") {
+ SECTION("Indexing") {
int const n0 = count_continuous<Decay2>::count;
int const n1 = count_continuous<ContinuousProcess3>::count;
- int const n2 = count_continuous<ContinuousProcess2, count_continuous<ContinuousProcess3>::count>::count;
- int const n1_b = count_continuous<Process2, count_continuous<ContinuousProcess3>::count>::count;
- int const n1_c = count_continuous<ContinuousProcess3, count_continuous<Process2>::count>::count;
- int const n12 = count_continuous<ContinuousProcess2, count_continuous<ContinuousProcess3, 10>::count>::count;
- int const n11_b = count_continuous<Process1, count_continuous<ContinuousProcess3, 10>::count>::count;
- int const n11_c = count_continuous<ContinuousProcess3, count_continuous<Process1, 10>::count>::count;
+ int const n2 = count_continuous<ContinuousProcess2,
+ count_continuous<ContinuousProcess3>::count>::count;
+ int const n1_b =
+ count_continuous<Process2, count_continuous<ContinuousProcess3>::count>::count;
+ int const n1_c =
+ count_continuous<ContinuousProcess3, count_continuous<Process2>::count>::count;
+ int const n12 =
+ count_continuous<ContinuousProcess2,
+ count_continuous<ContinuousProcess3, 10>::count>::count;
+ int const n11_b =
+ count_continuous<Process1,
+ count_continuous<ContinuousProcess3, 10>::count>::count;
+ int const n11_c = count_continuous<ContinuousProcess3,
+ count_continuous<Process1, 10>::count>::count;
CHECK(n0 == 0);
CHECK(n1 == 1);
@@ -506,30 +725,44 @@ TEST_CASE("Continuous Process Indexing", "ProcessSequence") {
CHECK(n11_c == 11);
CHECK(n12 == 12);
+ std::cout << count_continuous<ContinuousProcess3>::count << std::endl;
+ std::cout << count_continuous<Process3>::count << std::endl;
- struct TestSelect {
- SwitchResult operator()(const DummyData& p) const {
- std::cout << "TestSelect data=" << p.data_[0] << std::endl;
+ struct SwitchSelect {
+ SwitchResult operator()(DummyData const& p) const {
+ std::cout << "SwitchSelect data=" << p.data_[0] << std::endl;
if (p.data_[0] > 0) return SwitchResult::First;
return SwitchResult::Second;
}
};
- TestSelect select1;
- auto sequence1 = make_sequence(Process1(0), ContinuousProcess2(0), Decay1(0));
- auto sequence2 = make_sequence(ContinuousProcess3(0), Process2(0), Decay2(0));
+ auto sequence1 = make_sequence(Process1(0), ContinuousProcess2(0, 2_m), Decay1(0));
+ auto sequence2 = make_sequence(ContinuousProcess3(0, 3_m), Process2(0), Decay2(0),
+ ContinuousProcess1(0, 1_m));
+ SwitchSelect select1;
auto switch_seq = SwitchProcessSequence(sequence1, sequence2, select1);
- auto sequence3 = make_sequence(ContinuousProcess1(0), Process3(0),
- switch_seq);
+ auto sequence3 = make_sequence(ContinuousProcess1(0, 1_m), Process3(0), switch_seq);
+ auto sequence4 = make_sequence(ContinuousProcess1(0, 1_m), Process3(0),
+ SwitchProcessSequence(sequence1, sequence2, select1));
+ int const switch_seq_n = count_continuous<decltype(switch_seq)>::count;
+ int const sequence3_n = count_continuous<decltype(sequence3)>::count;
CHECK(decltype(sequence1)::nContinuous == 1);
CHECK(count_continuous<decltype(sequence1)>::count == 1);
- CHECK(count_continuous<decltype(sequence2)>::count == 1);
- CHECK(count_continuous<decltype(switch_seq)>::count == 2);
- CHECK(count_continuous<decltype(sequence3)>::count == 3);
-
+ CHECK(count_continuous<decltype(sequence2)>::count == 2);
+ CHECK(switch_seq_n == 3);
+ CHECK(sequence3_n == 4);
+ CHECK(count_continuous<decltype(sequence4)>::count == 4);
+
+ std::cout << "switch_seq "
+ << boost::typeindex::type_id<decltype(switch_seq)>().pretty_name()
+ << std::endl;
+
+ std::cout << "sequence3 "
+ << boost::typeindex::type_id<decltype(sequence3)>().pretty_name()
+ << std::endl;
}
}
diff --git a/tests/modules/testObservationPlane.cpp b/tests/modules/testObservationPlane.cpp
index 9b9be4656..90a12288b 100644
--- a/tests/modules/testObservationPlane.cpp
+++ b/tests/modules/testObservationPlane.cpp
@@ -23,10 +23,10 @@
using namespace corsika;
-TEST_CASE("ContinuousProcess interface", "[proccesses][observation_plane]") {
+TEST_CASE("ObservationPlane", "[proccesses][observation_plane]") {
- logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
+ logging::set_level(logging::level::trace);
+ corsika_logger->set_pattern("[%n:%^%-8l%$]: %v");
auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Oxygen);
auto const& cs = *csPtr;
@@ -34,48 +34,64 @@ TEST_CASE("ContinuousProcess interface", "[proccesses][observation_plane]") {
[[maybe_unused]] auto const& node_dummy = nodePtr;
/*
- Test with downward going 1_GeV neutrino, starting at 0,1_m,10m
+ Test with 1_GeV neutrino, starting at 0,0,0 travelling into +x direction (see
+ setup_stack).
- ObservationPlane has origin at 0,0,0
+ ObservationPlane has origin at 10,0,0 and a normal in x-direction
*/
auto [stack, viewPtr] =
setup::testing::setup_stack(Code::NuE, 0, 0, 1_GeV, nodePtr, cs);
[[maybe_unused]] setup::StackView& view = *viewPtr;
auto particle = stack->getNextParticle();
+ // dummy track. Not used for calculation!
Point const start(cs, {0_m, 1_m, 10_m});
VelocityVector vec(cs, 0_m / second, 0_m / second, -constants::c);
Line line(start, vec);
-
- setup::Trajectory track =
+ setup::Trajectory no_used_track =
setup::testing::make_track<setup::Trajectory>(line, 12_m / constants::c);
- particle.setPosition(Point(cs, {1_m, 1_m, 10_m})); // moving already along -z
-
SECTION("horizontal plane") {
- Plane const obsPlane(Point(cs, {0_m, 0_m, 0_m}), DirectionVector(cs, {0., 0., 1.}));
- ObservationPlane obs(obsPlane, DirectionVector(cs, {1., 0., 0.}), "particles.dat",
+ Plane const obsPlane(Point(cs, {10_m, 0_m, 0_m}), DirectionVector(cs, {1., 0., 0.}));
+ ObservationPlane obs(obsPlane, DirectionVector(cs, {0., 1., 0.}), "particles.dat",
true);
- LengthType const length = obs.getMaxStepLength(particle, track);
- ProcessReturn const ret = obs.doContinuous(particle, track);
+ LengthType const length = obs.getMaxStepLength(particle, no_used_track);
+ ProcessReturn const ret = obs.doContinuous(particle, no_used_track, true);
CHECK(length / 10_m == Approx(1).margin(1e-4));
CHECK(ret == ProcessReturn::ParticleAbsorbed);
}
- SECTION("inclined plane") {}
-
SECTION("transparent plane") {
- Plane const obsPlane(Point(cs, {0_m, 0_m, 0_m}), DirectionVector(cs, {0., 0., 1.}));
- ObservationPlane obs(obsPlane, DirectionVector(cs, {1., 0., 0.}), "particles.dat",
+ Plane const obsPlane(Point(cs, {1_m, 0_m, 0_m}), DirectionVector(cs, {1., 0., 0.}));
+ ObservationPlane obs(obsPlane, DirectionVector(cs, {0., 0., 1.}), "particles.dat",
false);
- LengthType const length = obs.getMaxStepLength(particle, track);
- ProcessReturn const ret = obs.doContinuous(particle, track);
+ LengthType const length = obs.getMaxStepLength(particle, no_used_track);
+ ProcessReturn const ret = obs.doContinuous(particle, no_used_track, true);
- CHECK(length / 10_m == Approx(1).margin(1e-4));
+ CHECK(length / 1_m == Approx(1).margin(1e-4));
CHECK(ret == ProcessReturn::Ok);
}
+
+ SECTION("inclined plane, inclined particle") {
+
+ MomentumVector const pnew = MomentumVector(cs, {1_GeV, 0.5_GeV, -0.4_GeV});
+ HEPEnergyType const enew = sqrt(pnew.dot(pnew));
+ particle.setMomentum(pnew);
+ particle.setEnergy(enew);
+
+ Plane const obsPlane(Point(cs, {10_m, 5_m, 5_m}),
+ DirectionVector(cs, {1, 0.1, -0.05}));
+ ObservationPlane obs(obsPlane, DirectionVector(cs, {0., 1., 0.}), "particles.dat",
+ true);
+
+ LengthType const length = obs.getMaxStepLength(particle, no_used_track);
+ ProcessReturn const ret = obs.doContinuous(particle, no_used_track, true);
+
+ CHECK(length / 10_m == Approx(1.1375).margin(1e-4));
+ CHECK(ret == ProcessReturn::ParticleAbsorbed);
+ }
}
diff --git a/tests/modules/testParticleCut.cpp b/tests/modules/testParticleCut.cpp
index 0dede160e..ef12f58e8 100644
--- a/tests/modules/testParticleCut.cpp
+++ b/tests/modules/testParticleCut.cpp
@@ -14,6 +14,7 @@
#include <corsika/framework/geometry/Vector.hpp>
#include <corsika/framework/utility/CorsikaFenv.hpp>
#include <corsika/media/Environment.hpp>
+#include <corsika/framework/process/ContinuousProcessIndex.hpp>
#include <SetupTestStack.hpp>
#include <SetupTestTrajectory.hpp>
--
GitLab