diff --git a/corsika/detail/framework/core/Cascade.inl b/corsika/detail/framework/core/Cascade.inl
index 8cd5b2e39f3f87e303ee51aecf9de64ecdda9131..0e9161259de194737b19c46504645b79233c5b3d 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/FourVector.inl b/corsika/detail/framework/geometry/FourVector.inl
index 6c0a039e5466ef46d2d3014b7e91a6f0621c922c..135790bfcfa83efaf1a78f218228d0c2557f5853 100644
--- a/corsika/detail/framework/geometry/FourVector.inl
+++ b/corsika/detail/framework/geometry/FourVector.inl
@@ -95,7 +95,7 @@ namespace corsika {
template <typename TTimeType, typename TSpaceVecType>
typename FourVector<TTimeType, TSpaceVecType>::norm_type
- FourVector<TTimeType, TSpaceVecType>::operator*(FourVector const& b) {
+ FourVector<TTimeType, TSpaceVecType>::operator*(FourVector const& b) {
if constexpr (std::is_same<time_type, decltype(std::declval<space_type>() / meter *
second)>::value)
return timeLike_ * b.timeLike_ * constants::cSquared - spaceLike_.norm();
diff --git a/corsika/detail/framework/geometry/Plane.inl b/corsika/detail/framework/geometry/Plane.inl
index e158e04d36a59bbec08e1b3e2b7e6aebc76ee7ab..d32167cf138254a3325b4fe9503dd4884bdbde79 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 6c550c1ce716deb8e3004f775a81b9bebb182d75..45c6900f58ff152400e9cc44c7419aee8efbcf01 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 << " (ref:" << fmt::ptr(p.getCoordinateSystem()) << ")";
+ return os;
+ }
+
} // namespace corsika
diff --git a/corsika/detail/framework/geometry/QuantityVector.inl b/corsika/detail/framework/geometry/QuantityVector.inl
index 9df35aab4827ecf034697b9554498eb45f7fb428..6ba02e4f8d4bd8e8cec5231d47ffb23d87931430 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 cbfd597b3e556fc04c03b0829b4f949bcb459ae8..ff968693f8542ef2f95629c4d7ca7050b406a549 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 << " (ref:" << fmt::ptr(v.getCoordinateSystem()) << ")";
+ return os;
+ }
+
} // namespace corsika
diff --git a/corsika/detail/framework/process/ProcessSequence.inl b/corsika/detail/framework/process/ProcessSequence.inl
index 7759270ee0b0cfe0fc416d2c5923b3b9dc1498e3..19056ad6c696c8c74268742e7f1447bdca0be856 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>
+ template <typename TProcess1, typename TProcess2, int IndexStart, int IndexProcess1,
+ int IndexProcess2>
template <typename TParticle>
- ProcessReturn ProcessSequence<TProcess1, TProcess2>::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>
+ template <typename TProcess1, typename TProcess2, int IndexStart, int IndexProcess1,
+ int IndexProcess2>
template <typename TParticle, typename TTrack>
- ProcessReturn ProcessSequence<TProcess1, TProcess2>::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>
+ template <typename TProcess1, typename TProcess2, int IndexStart, int IndexProcess1,
+ int IndexProcess2>
template <typename TSecondaries>
- void ProcessSequence<TProcess1, TProcess2>::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>
- bool ProcessSequence<TProcess1, TProcess2>::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>
+ template <typename TProcess1, typename TProcess2, int IndexStart, int IndexProcess1,
+ int IndexProcess2>
template <typename TStack>
- void ProcessSequence<TProcess1, TProcess2>::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>
+ /**
+ * 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>::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>
+ template <typename TProcess1, typename TProcess2, int IndexStart, int IndexProcess1,
+ int IndexProcess2>
template <typename TParticle>
- InverseGrammageType ProcessSequence<TProcess1, TProcess2>::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>
+ template <typename TProcess1, typename TProcess2, int IndexStart, int IndexProcess1,
+ int IndexProcess2>
template <typename TSecondaryView>
- inline ProcessReturn ProcessSequence<TProcess1, TProcess2>::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>
+ template <typename TProcess1, typename TProcess2, int IndexStart, int IndexProcess1,
+ int IndexProcess2>
template <typename TParticle>
- inline InverseTimeType ProcessSequence<TProcess1, TProcess2>::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>
+ template <typename TProcess1, typename TProcess2, int IndexStart, int IndexProcess1,
+ int IndexProcess2>
// select decay process
template <typename TSecondaryView>
- inline ProcessReturn ProcessSequence<TProcess1, TProcess2>::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 62d33995afefe53aa3865df31969413d18a63a99..2eadc8782bf164ee7074c97db4f35d37db2abbd9 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,10 +27,14 @@
namespace corsika {
- template <typename TProcess1, typename TProcess2, typename TSelect>
+ template <typename TProcess1, typename TProcess2, typename TSelect, int IndexStart,
+ int IndexProcess1, int IndexProcess2>
template <typename TParticle, typename TVTNType>
- ProcessReturn SwitchProcessSequence<TProcess1, TProcess2, TSelect>::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>,
@@ -50,24 +56,27 @@ namespace corsika {
return ProcessReturn::Ok;
}
- template <typename TProcess1, typename TProcess2, typename TSelect>
+ template <typename TProcess1, typename TProcess2, typename TSelect, int IndexStart,
+ int IndexProcess1, int IndexProcess2>
template <typename TParticle, typename TTrack>
- inline ProcessReturn SwitchProcessSequence<TProcess1, TProcess2, TSelect>::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;
}
@@ -75,10 +84,12 @@ namespace corsika {
return ProcessReturn::Ok;
}
- template <typename TProcess1, typename TProcess2, typename TSelect>
+ template <typename TProcess1, typename TProcess2, typename TSelect, int IndexStart,
+ int IndexProcess1, int IndexProcess2>
template <typename TSecondaries>
- inline void SwitchProcessSequence<TProcess1, TProcess2, TSelect>::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: {
@@ -100,39 +111,42 @@ namespace corsika {
}
}
- template <typename TProcess1, typename TProcess2, typename TSelect>
+ template <typename TProcess1, typename TProcess2, typename TSelect, int IndexStart,
+ int IndexProcess1, int IndexProcess2>
template <typename TParticle, typename TTrack>
- inline LengthType SwitchProcessSequence<TProcess1, TProcess2,
- TSelect>::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>
+ template <typename TProcess1, typename TProcess2, typename TSelect, int IndexStart,
+ int IndexProcess1, int IndexProcess2>
template <typename TParticle>
- inline InverseGrammageType
- SwitchProcessSequence<TProcess1, TProcess2, TSelect>::getInverseInteractionLength(
- TParticle&& particle) {
+ inline InverseGrammageType SwitchProcessSequence<
+ TProcess1, TProcess2, TSelect, IndexStart, IndexProcess1,
+ IndexProcess2>::getInverseInteractionLength(TParticle&& particle) {
switch (select_(particle)) {
case SwitchResult::First: {
@@ -155,12 +169,14 @@ namespace corsika {
return 0 * meter * meter / gram; // default value
}
- template <typename TProcess1, typename TProcess2, typename TSelect>
+ template <typename TProcess1, typename TProcess2, typename TSelect, int IndexStart,
+ int IndexProcess1, int IndexProcess2>
template <typename TSecondaryView>
- inline ProcessReturn
- SwitchProcessSequence<TProcess1, TProcess2, TSelect>::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) {
@@ -203,11 +219,12 @@ namespace corsika {
return ProcessReturn::Ok;
}
- template <typename TProcess1, typename TProcess2, typename TSelect>
+ template <typename TProcess1, typename TProcess2, typename TSelect, int IndexStart,
+ int IndexProcess1, int IndexProcess2>
template <typename TParticle>
inline InverseTimeType
- SwitchProcessSequence<TProcess1, TProcess2, TSelect>::getInverseLifetime(
- TParticle&& particle) {
+ SwitchProcessSequence<TProcess1, TProcess2, TSelect, IndexStart, IndexProcess1,
+ IndexProcess2>::getInverseLifetime(TParticle&& particle) {
switch (select_(particle)) {
case SwitchResult::First: {
@@ -229,12 +246,15 @@ namespace corsika {
return 0 / second; // default value
}
- template <typename TProcess1, typename TProcess2, typename TSelect>
+ 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>::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) {
@@ -277,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/media/WeightProvider.inl b/corsika/detail/media/WeightProvider.inl
index af9849aa3ca75b9bdaeefc661d56232f1ab895f4..c075b236350127c9978fb7c08d87ed035292e5a8 100644
--- a/corsika/detail/media/WeightProvider.inl
+++ b/corsika/detail/media/WeightProvider.inl
@@ -20,7 +20,7 @@ namespace corsika {
template <class AConstIterator, class BConstIterator>
typename WeightProviderIterator<AConstIterator, BConstIterator>::value_type
- WeightProviderIterator<AConstIterator, BConstIterator>::operator*() const {
+ WeightProviderIterator<AConstIterator, BConstIterator>::operator*() const {
return ((*aIter_) * (*bIter_)).magnitude();
}
diff --git a/corsika/detail/modules/LongitudinalProfile.inl b/corsika/detail/modules/LongitudinalProfile.inl
index d18813c26f08a7fc42dfaa0042579943b3185d69..df6bb471dbb04c4dfeff3a9c62cc86f3501f2ec1 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 e266a357f3c89351d6126de767bc09978717a34f..74c29b336dc26a3837b036db5bf3af895feb96f4 100644
--- a/corsika/detail/modules/ObservationPlane.inl
+++ b/corsika/detail/modules/ObservationPlane.inl
@@ -7,6 +7,8 @@
*/
#include <corsika/modules/ObservationPlane.hpp>
+#include <corsika/framework/core/Logging.hpp>
+#include <corsika/setup/SetupTrajectory.hpp>
#include <fstream>
@@ -26,27 +28,22 @@ namespace corsika {
}
ProcessReturn ObservationPlane::doContinuous(
- corsika::setup::Stack::particle_type& particle,
- corsika::setup::Trajectory& trajectory) {
- TimeType const timeOfIntersection =
- (plane_.getCenter() - trajectory.getPosition(0)).dot(plane_.getNormal()) /
- trajectory.getVelocity(0).dot(plane_.getNormal());
+ corsika::setup::Stack::particle_type& particle, corsika::setup::Trajectory&,
+ bool const stepLimit) {
- if (timeOfIntersection < TimeType::zero()) { return ProcessReturn::Ok; }
+ /*
+ The current step did not yet reach the ObservationPlane, do nothing now and wait:
+ */
+ if (!stepLimit) { return ProcessReturn::Ok; }
- if (plane_.isAbove(trajectory.getPosition(0)) ==
- plane_.isAbove(trajectory.getPosition(1))) {
- return ProcessReturn::Ok;
- }
-
- const auto energy = particle.getEnergy();
- auto const displacement = trajectory.getPosition(1) - plane_.getCenter();
+ 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
- << (trajectory.getPosition(1) - plane_.getCenter()).getNorm() / 1_m
- << std::endl;
+ << (pointOfIntersection - plane_.getCenter()).getNorm() / 1_m << '\n';
if (deleteOnHit_) {
count_ground_++;
@@ -59,84 +56,37 @@ namespace corsika {
}
LengthType ObservationPlane::getMaxStepLength(
- corsika::setup::Stack::particle_type const& vParticle,
+ corsika::setup::Stack::particle_type const& particle,
corsika::setup::Trajectory const& trajectory) {
- int chargeNumber;
- if (is_nucleus(vParticle.getPID())) {
- chargeNumber = vParticle.getNuclearZ();
- } else {
- chargeNumber = get_charge_number(vParticle.getPID());
- }
- auto const* currentLogicalVolumeNode = vParticle.getNode();
- auto magneticfield = currentLogicalVolumeNode->getModelProperties().getMagneticField(
- vParticle.getPosition());
- auto direction = trajectory.getVelocity(0).normalized();
-
- if (chargeNumber != 0 &&
- abs(plane_.getNormal().dot(
- trajectory.getLine().getVelocity().cross(magneticfield))) *
- 1_s / 1_m / 1_T >
- 1e-6) {
- auto const* currentLogicalVolumeNode = vParticle.getNode();
- auto magneticfield =
- currentLogicalVolumeNode->getModelProperties().getMagneticField(
- vParticle.getPosition());
- auto k =
- chargeNumber * constants::c * 1_eV / (vParticle.getMomentum().getNorm() * 1_V);
-
- if (direction.dot(plane_.getNormal()) * direction.dot(plane_.getNormal()) -
- (plane_.getNormal().dot(trajectory.getPosition(0) - plane_.getCenter()) *
- plane_.getNormal().dot(direction.cross(magneticfield)) * 2 * k) <
- 0) {
- return std::numeric_limits<double>::infinity() * 1_m;
- }
-
- LengthType MaxStepLength1 =
- (sqrt(direction.dot(plane_.getNormal()) * direction.dot(plane_.getNormal()) -
- (plane_.getNormal().dot(trajectory.getPosition(0) - plane_.getCenter()) *
- plane_.getNormal().dot(direction.cross(magneticfield)) * 2 * k)) -
- direction.dot(plane_.getNormal()) / direction.getNorm()) /
- (plane_.getNormal().dot(direction.cross(magneticfield)) * k);
-
- LengthType MaxStepLength2 =
- (-sqrt(direction.dot(plane_.getNormal()) * direction.dot(plane_.getNormal()) -
- (plane_.getNormal().dot(trajectory.getPosition(0) - plane_.getCenter()) *
- plane_.getNormal().dot(direction.cross(magneticfield)) * 2 * k)) -
- direction.dot(plane_.getNormal()) / direction.getNorm()) /
- (plane_.getNormal().dot(direction.cross(magneticfield)) * k);
-
- if (MaxStepLength1 <= 0_m && MaxStepLength2 <= 0_m) {
- return std::numeric_limits<double>::infinity() * 1_m;
- } else if (MaxStepLength1 <= 0_m || MaxStepLength2 < MaxStepLength1) {
- std::cout << " steplength to obs plane 2: " << MaxStepLength2 << std::endl;
- return MaxStepLength2 *
- (direction + direction.cross(magneticfield) * MaxStepLength2 * k / 2)
- .getNorm() *
- 1.001;
- } else if (MaxStepLength2 <= 0_m || MaxStepLength1 < MaxStepLength2) {
- std::cout << " steplength to obs plane 1: " << MaxStepLength1 << std::endl;
- return MaxStepLength1 *
- (direction + direction.cross(magneticfield) * MaxStepLength2 * k / 2)
- .getNorm() *
- 1.001;
- }
- }
+ 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());
- TimeType const timeOfIntersection =
- (plane_.getCenter() - trajectory.getPosition(0)).dot(plane_.getNormal()) /
- trajectory.getVelocity(0).dot(plane_.getNormal());
+ TimeType const timeOfIntersection = intersection.getEntry();
+
+ 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;
+ }
double const fractionOfIntersection = timeOfIntersection / trajectory.getDuration();
auto const pointOfIntersection = trajectory.getPosition(fractionOfIntersection);
- auto dist = (trajectory.getPosition(0) - pointOfIntersection).getNorm() * 1.0001;
- CORSIKA_LOG_TRACE("ObservationPlane w/o magnetic field: 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 7a22b4d719beddef778acab9861b65ba7578217f..3fd882225fec411ced5ef823dff9b35fcea0ef56 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 27b47cf6571c23975842fdec279e582ed30ce775..c6f72edf437b4d49cf728cf3ab16b358176f78f9 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 89fe81c7bee6316b97c9cd78322dfef247eade3d..cbabb0892dc3c4880c74d910ddef99ce7036736c 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 8bfbd705b8bb0dfc19d1e7efe79252aa5a6a69be..26f604c998efa2360d4da94adab53e9fac4d28b9 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 2c284b4d01dddc3be3322dc305a84f6f836bd5e5..5869601a566d959ffac3432424b5fb9eaa4929cf 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/core/Logging.hpp b/corsika/framework/core/Logging.hpp
index 4d4c43a7ba2e4f7e9d2c90e87a905d542b554df8..29ef3e13c55470f95d933f4bcb6e998cdef68fef 100644
--- a/corsika/framework/core/Logging.hpp
+++ b/corsika/framework/core/Logging.hpp
@@ -35,9 +35,9 @@
// if this is a Debug build, include debug messages in objects
#ifdef DEBUG
// trace is the highest level of logging (ALL messages will be printed)
-#define SPDLOG_ACTIVE_LEVEL SPDLOG_LEVEL_DEBUG
+#define SPDLOG_ACTIVE_LEVEL SPDLOG_LEVEL_TRACE
#else // otherwise, remove everything but "error" and worse messages
-#define SPDLOG_ACTIVE_LEVEL SPDLOG_LEVEL_INFO
+#define SPDLOG_ACTIVE_LEVEL SPDLOG_LEVEL_DEBUG
#endif
#include <spdlog/fmt/ostr.h> // will output whenerver a streaming operator is found
diff --git a/corsika/framework/geometry/Plane.hpp b/corsika/framework/geometry/Plane.hpp
index 9b2eacade48f5d9913d5b5311f1ae595cbcde039..7218e6b3564e857946669743be61c6519f467f12 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 235237afe59b8ffbf19413e625e7bdd1178a18e8..eb6c2f11e2679d434e2dbf617875254a60af621f 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 e91651d12ea99cafd40718ae288a869852a5cf35..b8c96a54faa504a84c79272ce23f67534b6950b8 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
@@ -39,8 +43,28 @@ namespace corsika {
const TDerived& ref() const { return static_cast<const TDerived&>(*this); }
public:
+ //! Default number of processes ist just one, obviously
+ static unsigned int constexpr getNumberOfProcesses() { return 1; }
+
// Base processor type for use in other template classes
using process_type = TDerived;
};
+ /**
+ * ProcessTraits specialization
+ **/
+ template <typename TProcess>
+ struct is_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 {};
+
+ template <typename TProcess, int N>
+ struct count_processes<TProcess, N,
+ typename std::enable_if_t<is_process_v<TProcess> &&
+ !is_process_sequence_v<TProcess>>> {
+ static unsigned int constexpr count = N + 1;
+ };
+
} // namespace corsika
diff --git a/corsika/framework/process/ContinuousProcess.hpp b/corsika/framework/process/ContinuousProcess.hpp
index 1343c082272a8e0e49830e53c164f985abe5e860..51d73d332fbe00ea3c5845608a0c582f224d5b4c 100644
--- a/corsika/framework/process/ContinuousProcess.hpp
+++ b/corsika/framework/process/ContinuousProcess.hpp
@@ -8,8 +8,10 @@
#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>
namespace corsika {
@@ -28,13 +30,41 @@ 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&) const;
+ 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;
};
+ /**
+ * ProcessTraits specialization
+ **/
+ 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>>> {
+ 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 0000000000000000000000000000000000000000..bc2617d2bb6ed5366431d1369480ef9abf593426
--- /dev/null
+++ b/corsika/framework/process/ContinuousProcessIndex.hpp
@@ -0,0 +1,32 @@
+/*
+ * (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 {
+
+ /**
+ * To index individual processes (continuous processes) inside a
+ * ProcessSequence.
+ *
+ **/
+
+ 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 0000000000000000000000000000000000000000..32b921857aaf6762c584ac5bf68b6885c019e038
--- /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 70337e1099d09cd4945d7f557f2a1661dedc0d0b..d909d07151c4d15a1bdf040af11ce46e8a930f24 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,21 @@
namespace corsika {
+ template <typename TProcess, int N>
+ struct count_continuous<TProcess, N,
+ typename std::enable_if_t<is_process_sequence_v<TProcess>>> {
+ static unsigned int constexpr count =
+ N + std::decay_t<TProcess>::getNumberOfProcesses();
+ };
+
+ template <typename TProcess, int N>
+ struct count_processes<TProcess, N,
+ typename std::enable_if_t<is_process_v<TProcess> &&
+ is_process_sequence_v<TProcess>>> {
+ static unsigned int constexpr count =
+ N + std::decay_t<TProcess>::getNumberOfProcesses();
+ };
+
/**
*
* Definition of a static process list/sequence
@@ -39,13 +56,21 @@ 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 CRTP.
+ * (The sequence, and the processes use the CRTP, curiously recurring template
+ * pattern).
*
- * \todo There are several FIXME's in the ProcessSequence.inl due to
- * outstanding migration of SecondaryView::parent()
+ * 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
+ * - ProcessIndexOffset, IndexOfProcess1, IndexOfProcess2 are to count and index each
+ *ContinuousProcess in the entire process-chain
**/
- template <typename TProcess1, typename TProcess2 = NullModel>
+ template <typename TProcess1, typename TProcess2 = NullModel,
+ int ProcessIndexOffset = 0,
+ int IndexOfProcess1 = count_processes<
+ TProcess1, count_processes<TProcess2, ProcessIndexOffset>::count>::count,
+ int IndexOfProcess2 = count_processes<TProcess2, ProcessIndexOffset>::count>
class ProcessSequence : public BaseProcess<ProcessSequence<TProcess1, TProcess2>> {
using process1_type = typename std::decay_t<TProcess1>;
@@ -58,16 +83,13 @@ 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_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_process_v<process2_type>,
"can only use process derived from BaseProcess in "
"ProcessSequence, for Process 2");
- TProcess1 A_; /// process/list A, this is a reference, if possible
- TProcess2 B_; /// process/list B, this is a reference, if possible
-
public:
// resource management
ProcessSequence() = delete; // only initialized objects
@@ -96,7 +118,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,8 +140,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 ContinuousProcessStepLength getMaxStepLength(TParticle& particle,
+ TTrack& vTrack);
template <typename TParticle>
inline GrammageType getInteractionLength(TParticle&& particle) {
@@ -147,6 +181,17 @@ namespace corsika {
inline ProcessReturn selectDecay(
TSecondaryView& view, [[maybe_unused]] InverseTimeType decay_inv_select,
[[maybe_unused]] InverseTimeType decay_inv_sum = InverseTimeType::zero());
+
+ /**
+ * static counter to uniquely index (count) all ContinuousProcess in switch sequence.
+ **/
+ static unsigned int constexpr getNumberOfProcesses() { return numberOfProcesses_; }
+
+ private:
+ TProcess1 A_; /// process/list A, this is a reference, if possible
+ TProcess2 B_; /// process/list B, this is a reference, if possible
+
+ static unsigned int constexpr numberOfProcesses_ = IndexOfProcess1; // static counter
};
/**
@@ -170,12 +215,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_process_v<typename std::decay_t<TProcess1>>,
ProcessSequence<TProcess1, decltype(make_sequence(std::declval<TProcesses>()...))>>
make_sequence(TProcess1&& vA, TProcesses&&... vBs) {
return ProcessSequence<TProcess1,
@@ -192,12 +237,9 @@ 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>>
+ inline typename std::enable_if_t<is_process_v<typename std::decay_t<TProcess1>> &&
+ is_process_v<typename std::decay_t<TProcess2>>,
+ ProcessSequence<TProcess1, TProcess2>>
make_sequence(TProcess1&& vA, TProcess2&& vB) {
return ProcessSequence<TProcess1, TProcess2>(vA, vB);
}
@@ -211,10 +253,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_process_v<typename std::decay_t<TProcess>>,
+ ProcessSequence<TProcess, NullModel>>
make_sequence(TProcess&& vA) {
return ProcessSequence<TProcess, NullModel>(vA, NullModel());
}
@@ -226,11 +266,7 @@ 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_process_v<TProcess1> && is_process_v<TProcess2>, int>>
is_process_sequence() {}
};
diff --git a/corsika/framework/process/ProcessTraits.hpp b/corsika/framework/process/ProcessTraits.hpp
index 7b303d2d3ce3a9ce1363d2433e08832e80032f46..b08318093e7be4147e5408661b16c8feb1c3921e 100644
--- a/corsika/framework/process/ProcessTraits.hpp
+++ b/corsika/framework/process/ProcessTraits.hpp
@@ -16,6 +16,24 @@
namespace corsika {
+ /**
+ * A traits marker to identify BaseProcess, thus any type of process
+ */
+ template <typename TProcess, typename Enable = void>
+ struct is_process : std::false_type {};
+
+ template <typename TProcess>
+ bool constexpr is_process_v = is_process<TProcess>::value;
+
+ /**
+ * 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
**/
@@ -44,4 +62,20 @@ namespace corsika {
template <typename TClass>
bool constexpr contains_stack_process_v = contains_stack_process<TClass>::value;
+ /**
+ * traits class to count ContinuousProcess-es, general version
+ **/
+ template <typename TProcess, int N = 0, typename Enable = void>
+ struct count_continuous {
+ static unsigned int constexpr count = N;
+ };
+
+ /**
+ * traits class to count any type of Process, general version
+ **/
+ template <typename TProcess, int N = 0, typename Enable = void>
+ struct count_processes {
+ static unsigned int constexpr count = N;
+ };
+
} // namespace corsika
diff --git a/corsika/framework/process/SwitchProcessSequence.hpp b/corsika/framework/process/SwitchProcessSequence.hpp
index dab419ae5a5c66fda13eee9c0247c697e0cdcfd1..dcf4b214ec84f58bd69a44cf9cd47c2ae3f02a36 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>
@@ -57,10 +58,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
+ - IndexFirstProcess, IndexOfProcess1, IndexOfProcess2 are to count and index each
+ ContinuousProcess in the entire process-chain
+
See also class \sa ProcessSequence
**/
- template <typename TProcess1, typename TProcess2, typename TSelect>
+ template <typename TProcess1, typename TProcess2, typename TSelect,
+ int IndexFirstProcess = 0,
+ int IndexOfProcess1 = count_processes<TProcess1, IndexFirstProcess>::count,
+ int IndexOfProcess2 = count_processes<TProcess2, IndexOfProcess1>::count>
class SwitchProcessSequence
: public BaseProcess<SwitchProcessSequence<TProcess1, TProcess2, TSelect>> {
@@ -71,10 +81,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_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_process_v<process2_type>,
"can only use process derived from BaseProcess in "
"SwitchProcessSequence, for Process 2");
@@ -122,13 +132,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) {
@@ -158,12 +170,19 @@ namespace corsika {
TSecondaryView& view, [[maybe_unused]] InverseTimeType decay_inv_select,
[[maybe_unused]] InverseTimeType decay_inv_sum = InverseTimeType::zero());
+ /**
+ * static counter to uniquely index (count) all ContinuousProcess in switch sequence.
+ **/
+ static unsigned int constexpr getNumberOfProcesses() { return numberOfProcesses_; }
+
private:
TSelect select_; /// selector functor to switch between branch a and b, this is a
/// reference, if possible
TProcess1 A_; /// process branch a, this is a reference, if possible
TProcess2 B_; /// process branch b, this is a reference, if possible
+
+ static unsigned int constexpr numberOfProcesses_ = IndexOfProcess2; // static counter
};
/**
@@ -178,12 +197,9 @@ 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>>
+ inline typename std::enable_if_t<is_process_v<typename std::decay_t<TProcess1>> &&
+ is_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);
}
diff --git a/corsika/modules/LongitudinalProfile.hpp b/corsika/modules/LongitudinalProfile.hpp
index baaa5046824c142694097a18ba662ba97647799c..b06e2ef852541342f474881a0209a154cf85602e 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/ObservationPlane.hpp b/corsika/modules/ObservationPlane.hpp
index 7effcd93b6abfda4ae1cc8b2adca07b38cb6cdc6..8276c0d1235f44539bfc17468f20aff5eda902b5 100644
--- a/corsika/modules/ObservationPlane.hpp
+++ b/corsika/modules/ObservationPlane.hpp
@@ -31,7 +31,8 @@ namespace corsika {
bool = true);
ProcessReturn doContinuous(corsika::setup::Stack::particle_type& vParticle,
- corsika::setup::Trajectory& vTrajectory);
+ corsika::setup::Trajectory& vTrajectory,
+ bool const stepLimit);
LengthType getMaxStepLength(corsika::setup::Stack::particle_type const&,
corsika::setup::Trajectory const& vTrajectory);
diff --git a/corsika/modules/ParticleCut.hpp b/corsika/modules/ParticleCut.hpp
index 1678c528edbc9815bc9c1388bcd81d9c32748615..9842dd034530d7528e1359e63983d8adc85b25f0 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 c437bc7d15dfe4ada3ef269e874cccf2ce5a36a7..5bf9856b389fe80774ea629b27a4b4dc10f33c85 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 f8e8759376a3833bd3f7df438b2dbe4132c5e3c8..582a0c4e03046cc738105d32d2cfcdd7e94964a9 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 3b8976f3ac6bfe7a6d5251740782e4ca8fec89bc..0efbf12ecd6100ae3e6cc7e9bc777cacb2f2dfcb 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/em_shower.cpp b/examples/em_shower.cpp
index 2afddc84ee9744909d4378eb3cf2456267655308..e4d8b3cc432a1852f96dbababca6a0c024fc037b 100644
--- a/examples/em_shower.cpp
+++ b/examples/em_shower.cpp
@@ -141,11 +141,10 @@ int main(int argc, char** argv) {
// setup processes, decays and interactions
- // PROPOSAL processs proposal{...};
ParticleCut cut(10_GeV, 10_GeV, 100_PeV, 100_PeV, true);
- corsika::proposal::Interaction proposal(env);
- corsika::proposal::ContinuousProcess em_continuous(env);
- InteractionCounter proposalCounted(proposal);
+ corsika::proposal::Interaction emCascade(env);
+ corsika::proposal::ContinuousProcess emContinuous(env);
+ InteractionCounter emCascadeCounted(emCascade);
TrackWriter trackWriter("tracks.dat");
@@ -156,7 +155,7 @@ int main(int argc, char** argv) {
ObservationPlane observationLevel(obsPlane, DirectionVector(rootCS, {1., 0., 0.}),
"particles.dat");
- auto sequence = make_sequence(proposalCounted, em_continuous, longprof, cut,
+ auto sequence = make_sequence(emCascadeCounted, emContinuous, longprof, cut,
observationLevel, trackWriter);
// define air shower object, run simulation
setup::Tracking tracking;
@@ -169,18 +168,18 @@ int main(int argc, char** argv) {
EAS.run();
cut.showResults();
- em_continuous.showResults();
+ emContinuous.showResults();
observationLevel.showResults();
const HEPEnergyType Efinal = cut.getCutEnergy() + cut.getInvEnergy() +
- cut.getEmEnergy() + em_continuous.getEnergyLost() +
+ cut.getEmEnergy() + emContinuous.getEnergyLost() +
observationLevel.getEnergyGround();
cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl
<< "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
observationLevel.reset();
cut.reset();
- em_continuous.reset();
+ emContinuous.reset();
- auto const hists = proposalCounted.getHistogram();
+ auto const hists = emCascadeCounted.getHistogram();
save_hist(hists.labHist(), "inthist_lab_emShower.npz", true);
save_hist(hists.CMSHist(), "inthist_cms_emShower.npz", true);
longprof.save("longprof_emShower.txt");
diff --git a/examples/particle_list_example.cpp b/examples/particle_list_example.cpp
index 4b8d1d75770743d92e7ad2051ce05bb25d1156c9..8a10d9f094fcc5e4838778e2bc5526347b2bfba0 100644
--- a/examples/particle_list_example.cpp
+++ b/examples/particle_list_example.cpp
@@ -35,37 +35,39 @@ using namespace std;
int main() {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
+ corsika_logger->set_pattern("[%n:%^%-8l%$] %v");
logging::info(
+ "\n"
"------------------------------------------\n"
- "particles in CORSIKA\n"
- "------------------------------------------\n"
- "Name | "
+ " particles in CORSIKA\n"
+ "------------------------------------------\n");
+ int const width = 20 + 10 + 10 + 10 + 15 + 15 + 17;
+ logging::info(
+ "Name | "
"PDG-id | "
"SIBYLL-id | "
"QGSJETII-id| "
"PDG-mass (GeV) | "
- "SIBYLL-mass (GeV)|\n"
- "{:-}",
- "", 104);
+ "SIBYLL-mass (GeV)|");
+ logging::info("{:->{}}", ' ', width);
for (auto p : get_all_particles()) {
if (!is_nucleus(p)) {
corsika::sibyll::SibyllCode sib_id = corsika::sibyll::convertToSibyll(p);
auto const sib_mass = (sib_id != corsika::sibyll::SibyllCode::Unknown
? to_string(corsika::sibyll::getSibyllMass(p) / 1_GeV)
- : "--");
+ : "");
auto const qgs_id = corsika::qgsjetII::convertToQgsjetII(p);
- logging::info("{:20} | {:10} | {:10} | {:10.5} | {:18.5}", p,
+ logging::info("{:20} | {:10} | {:10} | {:10} | {:>15.5} | {:>15.5} |", p,
static_cast<int>(get_PDG(p)),
(sib_id != corsika::sibyll::SibyllCode::Unknown
? to_string(static_cast<int>(sib_id))
- : "--"),
+ : ""),
(qgs_id != corsika::qgsjetII::QgsjetIICode::Unknown
? to_string(static_cast<int>(qgs_id))
- : "--"),
+ : ""),
get_mass(p) / 1_GeV, sib_mass);
}
}
- logging::info("{:-104}", "");
+ logging::info("{:->{}}", ' ', width);
}
diff --git a/examples/staticsequence_example.cpp b/examples/staticsequence_example.cpp
index d7b9e06b5051d8b8a721b87bba5c9407f3090a70..3f6c5aff75c72ec9058773fabbb6347e0db15274 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 f012440e82deb19c3ab2b5220622c0255da9dbc5..04ddbb2221c9c753e7dbe9550dc07b9a3e45a77e 100644
--- a/examples/vertical_EAS.cpp
+++ b/examples/vertical_EAS.cpp
@@ -91,8 +91,8 @@ using MyExtraEnv = MediumPropertyModel<UniformMagneticField<T>>;
int main(int argc, char** argv) {
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
- logging::set_level(logging::level::info);
+ corsika_logger->set_pattern("[%n:%^%-8l%$] %s:%#: %v");
+ logging::set_level(logging::level::trace);
CORSIKA_LOG_INFO("vertical_EAS");
@@ -222,9 +222,9 @@ int main(int argc, char** argv) {
decaySibyll.printDecayConfig();
ParticleCut cut{60_GeV, 60_GeV, 60_GeV, 60_GeV, true};
- corsika::proposal::Interaction proposal(env);
- corsika::proposal::ContinuousProcess em_continuous(env);
- InteractionCounter proposalCounted(proposal);
+ corsika::proposal::Interaction emCascade(env);
+ corsika::proposal::ContinuousProcess emContinuous(env);
+ InteractionCounter emCascadeCounted(emCascade);
OnShellCheck reset_particle_mass(1.e-3, 1.e-1, false);
TrackWriter trackWriter("tracks.dat");
@@ -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, proposalCounted, em_continuous, cut,
- trackWriter, observationLevel, longprof);
+ auto sequence =
+ make_sequence(stackInspect, hadronSequence, reset_particle_mass, decaySequence,
+ emContinuous, cut, trackWriter, observationLevel, longprof);
// define air shower object, run simulation
setup::Tracking tracking;
@@ -268,19 +268,19 @@ int main(int argc, char** argv) {
EAS.run();
cut.showResults();
- em_continuous.showResults();
+ emContinuous.showResults();
observationLevel.showResults();
const HEPEnergyType Efinal = cut.getCutEnergy() + cut.getInvEnergy() +
- cut.getEmEnergy() + em_continuous.getEnergyLost() +
+ cut.getEmEnergy() + emContinuous.getEnergyLost() +
observationLevel.getEnergyGround();
cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl
<< "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
observationLevel.reset();
cut.reset();
- em_continuous.reset();
+ emContinuous.reset();
auto const hists = sibyllCounted.getHistogram() + sibyllNucCounted.getHistogram() +
- urqmdCounted.getHistogram() + proposalCounted.getHistogram();
+ urqmdCounted.getHistogram();
save_hist(hists.labHist(), "inthist_lab_verticalEAS.npz", true);
save_hist(hists.CMSHist(), "inthist_cms_verticalEAS.npz", true);
diff --git a/modules/conex b/modules/conex
index 73e3d1fa850d2d3602a15b3b948ac1789fbc968d..820f042b6a055276d465437c74160ef7c199b646 160000
--- a/modules/conex
+++ b/modules/conex
@@ -1 +1 @@
-Subproject commit 73e3d1fa850d2d3602a15b3b948ac1789fbc968d
+Subproject commit 820f042b6a055276d465437c74160ef7c199b646
diff --git a/modules/data b/modules/data
index afe83dc19c1464deb176f38c3ddab80a64e081f4..8b76a9ca2599cd0ce1f204b17362eb06bbcf5277 160000
--- a/modules/data
+++ b/modules/data
@@ -1 +1 @@
-Subproject commit afe83dc19c1464deb176f38c3ddab80a64e081f4
+Subproject commit 8b76a9ca2599cd0ce1f204b17362eb06bbcf5277
diff --git a/tests/framework/testGeometry.cpp b/tests/framework/testGeometry.cpp
index a39201c04de97d13ccc24343cd5b4a0148a2c5a5..e8847ee5e824242344147f46afeb37dfe0eb3708 100644
--- a/tests/framework/testGeometry.cpp
+++ b/tests/framework/testGeometry.cpp
@@ -26,18 +26,20 @@ using namespace corsika::testing;
double constexpr absMargin = 1.0e-8;
-TEST_CASE("transformations between CoordinateSystems") {
+TEST_CASE("Geometry CoordinateSystems") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
+ corsika_logger->set_pattern("[%n:%^%-8l%$] %v");
CoordinateSystemPtr rootCS = get_root_CoordinateSystem();
QuantityVector<length_d> const coordinates{0_m, 0_m, 0_m};
Point p1(rootCS, coordinates);
+ CORSIKA_LOG_INFO("Point p1={}", p1);
QuantityVector<magnetic_flux_density_d> components{1. * tesla, 0. * tesla, 0. * tesla};
Vector<magnetic_flux_density_d> v1(rootCS, components);
+ CORSIKA_LOG_INFO("Vector<magnetic_flux_density_d> v1={}", v1);
CHECK((p1.getCoordinates() - coordinates).getNorm().magnitude() ==
Approx(0).margin(absMargin));
@@ -46,6 +48,7 @@ TEST_CASE("transformations between CoordinateSystems") {
SECTION("translations") {
QuantityVector<length_d> const translationVector{0_m, 4_m, 0_m};
+ CORSIKA_LOG_INFO("QuantityVector<length_d> translationVector={}", translationVector);
CoordinateSystemPtr translatedCS = make_translation(rootCS, translationVector);
@@ -184,7 +187,7 @@ TEST_CASE("transformations between CoordinateSystems") {
}
}
-TEST_CASE("CoordinateSystem hirarchy") {
+TEST_CASE("Geometry CoordinateSystem-hirarchy") {
CoordinateSystemPtr rootCS = get_root_CoordinateSystem();
@@ -232,7 +235,7 @@ TEST_CASE("CoordinateSystem hirarchy") {
CHECK((p1 - p6).getNorm().magnitude() == Approx(0).margin(absMargin));
}
-TEST_CASE("Sphere") {
+TEST_CASE("Geometry Sphere") {
CoordinateSystemPtr const& rootCS = get_root_CoordinateSystem();
Point center(rootCS, {0_m, 3_m, 4_m});
Sphere sphere(center, 5_m);
@@ -251,7 +254,7 @@ TEST_CASE("Sphere") {
}
}
-TEST_CASE("Trajectories") {
+TEST_CASE("Geometry Trajectories") {
CoordinateSystemPtr rootCS = get_root_CoordinateSystem();
Point r0(rootCS, {0_m, 0_m, 0_m});
diff --git a/tests/framework/testProcessSequence.cpp b/tests/framework/testProcessSequence.cpp
index ac80847cabbd24913175b44a6a78e62af7e9db2c..d2612bb8860497c34aa1550b1eb227ba2d108c68 100644
--- a/tests/framework/testProcessSequence.cpp
+++ b/tests/framework/testProcessSequence.cpp
@@ -9,6 +9,8 @@
#include <corsika/framework/process/ProcessSequence.hpp>
#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>
@@ -17,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
@@ -31,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++;
@@ -115,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++;
@@ -139,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++;
@@ -163,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; }
@@ -187,7 +233,7 @@ private:
class Decay1 : public DecayProcess<Decay1> {
public:
- Decay1(const int) {
+ Decay1(int const) {
cout << "Decay1()" << endl;
globalCount++;
}
@@ -204,7 +250,7 @@ public:
class Decay2 : public DecayProcess<Decay2> {
public:
- Decay2(const int) {
+ Decay2(int const) {
cout << "Decay2()" << endl;
globalCount++;
}
@@ -221,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&) {
@@ -234,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) {}
@@ -247,10 +297,22 @@ struct DummyView {
DummyData& parent() { return p_; }
};
-TEST_CASE("Process Sequence", "[Process Sequence]") {
+TEST_CASE("ProcessSequence General", "ProcessSequence") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
+ corsika_logger->set_pattern("[%n:%^%-8l%$]: %v");
+
+ SECTION("BaseProcess") {
+
+ Process1 m1(0);
+ const Process4 m4(3);
+
+ CHECK(is_process_v<Process1>);
+ CHECK_FALSE(is_process_v<DummyData>);
+ CHECK(is_process_v<decltype(m4)>);
+ CHECK(is_process_v<decltype(Decay1(1))>);
+ CHECK(is_process_v<decltype(ContinuousProcess3{3, 3_m})>);
+ }
SECTION("Check construction") {
globalCount = 0;
@@ -264,10 +326,15 @@ TEST_CASE("Process Sequence", "[Process Sequence]") {
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_process_v<decltype(sequence1)>);
+ CHECK(is_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);
@@ -278,7 +345,7 @@ TEST_CASE("Process Sequence", "[Process Sequence]") {
SECTION("interaction length") {
globalCount = 0;
- ContinuousProcess1 cp1(0);
+ ContinuousProcess1 cp1(0, 1_m);
Process2 m2(1);
Process3 m3(2);
@@ -296,7 +363,7 @@ TEST_CASE("Process Sequence", "[Process Sequence]") {
SECTION("lifetime") {
globalCount = 0;
- ContinuousProcess1 cp1(0);
+ ContinuousProcess1 cp1(0, 1_m);
Process2 m2(1);
Process3 m3(2);
Decay1 d3(3);
@@ -313,30 +380,61 @@ TEST_CASE("Process Sequence", "[Process Sequence]") {
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;
@@ -355,14 +453,14 @@ TEST_CASE("Process Sequence", "[Process Sequence]") {
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);
@@ -373,39 +471,64 @@ TEST_CASE("Process Sequence", "[Process Sequence]") {
}
}
-TEST_CASE("Switch Process Sequence", "[Process Sequence]") {
+TEST_CASE("SwitchProcessSequence", "ProcessSequence") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
+ 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;
@@ -416,7 +539,7 @@ TEST_CASE("Switch Process Sequence", "[Process Sequence]") {
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);
@@ -427,7 +550,7 @@ TEST_CASE("Switch Process Sequence", "[Process Sequence]") {
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);
@@ -478,4 +601,168 @@ TEST_CASE("Switch Process Sequence", "[Process Sequence]") {
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("ProcessSequence Indexing", "ProcessSequence") {
+
+ logging::set_level(logging::level::info);
+ corsika_logger->set_pattern("[%n:%^%-8l%$]: %v");
+
+ 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;
+
+ CHECK(n0 == 0);
+ CHECK(n1 == 1);
+ CHECK(n1_b == 1);
+ CHECK(n1_c == 1);
+ CHECK(n2 == 2);
+ CHECK(n11_b == 11);
+ CHECK(n11_c == 11);
+ CHECK(n12 == 12);
+
+ std::cout << count_continuous<ContinuousProcess3>::count << std::endl;
+ std::cout << count_continuous<Process3>::count << 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;
+ }
+ };
+
+ 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, 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)::getNumberOfProcesses() == 3);
+ CHECK(count_continuous<decltype(sequence1)>::count == 3);
+ CHECK(count_continuous<decltype(sequence2)>::count == 4);
+ CHECK(switch_seq_n == 7);
+ CHECK(sequence3_n == 9);
+ CHECK(count_continuous<decltype(sequence4)>::count == 9);
+
+ 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 9b9be4656d447b76552d63ae8e6b83657effa3d2..90a12288bd0b24c4c3021adaf66630fd1265d688 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 0dede160e563c6abab6a9a3fb2b1f0318f3dc1c7..ef12f58e8489341ad64adcfd0d67194c064514cb 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>