diff --git a/corsika/detail/framework/core/Cascade.inl b/corsika/detail/framework/core/Cascade.inl
index 07ec495469d8b1ef64f896f2b6facb3853ed60c0..64ca0a3aa5ff2f9fb0a95062d27ba4dcdba4807c 100644
--- a/corsika/detail/framework/core/Cascade.inl
+++ b/corsika/detail/framework/core/Cascade.inl
@@ -31,10 +31,12 @@ namespace corsika {
setNodes(); // put each particle on stack in correct environment volume
while (!stack_.isEmpty()) {
+
+ sequence_.initCascadeEquations();
+
while (!stack_.isEmpty()) {
CORSIKA_LOG_TRACE("Stack: {}", stack_.asString());
count_++;
-
auto pNext = stack_.getNextParticle();
CORSIKA_LOG_TRACE(
@@ -46,9 +48,10 @@ namespace corsika {
step(pNext);
sequence_.doStack(stack_);
}
+
// do cascade equations, which can put new particles on Stack,
// thus, the double loop
- // doCascadeEquations();
+ sequence_.doCascadeEquations(stack_);
}
}
diff --git a/corsika/detail/framework/process/CascadeEquationsProcess.hpp b/corsika/detail/framework/process/CascadeEquationsProcess.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..94e973c7888b15e1d68c209a2a0369cdd387618b
--- /dev/null
+++ b/corsika/detail/framework/process/CascadeEquationsProcess.hpp
@@ -0,0 +1,60 @@
+/*
+ * (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/process/ProcessTraits.hpp>
+#include <corsika/framework/utility/HasMethodSignature.hpp>
+
+/**
+ * @file CascadeEquationsProcess.hpp
+ */
+
+namespace corsika {
+
+ /**
+ * traits test for CascadeEquationsProcess::doCascadeEquations method.
+ */
+ template <class TProcess, typename TReturn, typename... TArg>
+ struct has_method_doCascadeEquations
+ : public detail::has_method_signature<TReturn, TArg...> {
+
+ //! method signature
+ using detail::has_method_signature<TReturn, TArg...>::testSignature;
+
+ //! the default value
+ template <class T>
+ static std::false_type test(...);
+
+ //! templated parameter option
+ template <class T>
+ static decltype(testSignature(&T::template doCascadeEquations<TArg...>)) test(
+ std::nullptr_t);
+
+ //! non templated parameter option
+ template <class T>
+ static decltype(testSignature(&T::doCascadeEquations)) test(std::nullptr_t);
+
+ public:
+ /**
+ * @name traits results
+ * @{
+ */
+ using type = decltype(test<std::decay_t<TProcess>>(nullptr));
+ static const bool value = type::value;
+ //! @}
+ };
+
+ /**
+ * value traits type.
+ */
+ template <class TProcess, typename TReturn, typename... TArg>
+ bool constexpr has_method_doCascadeEquations_v =
+ has_method_doCascadeEquations<TProcess, TReturn, TArg...>::value;
+
+} // namespace corsika
diff --git a/corsika/detail/framework/process/ProcessSequence.inl b/corsika/detail/framework/process/ProcessSequence.inl
index 5eb3d20ae785b7645101b362cb4697282a29ecdd..4d381abb16af98729d6defdfb53bfc0503d420f9 100644
--- a/corsika/detail/framework/process/ProcessSequence.inl
+++ b/corsika/detail/framework/process/ProcessSequence.inl
@@ -19,6 +19,7 @@
#include <corsika/framework/process/ProcessReturn.hpp>
#include <corsika/framework/process/SecondariesProcess.hpp>
#include <corsika/framework/process/StackProcess.hpp>
+#include <corsika/framework/process/CascadeEquationsProcess.hpp>
#include <cmath>
#include <limits>
@@ -331,6 +332,79 @@ namespace corsika {
return tot;
}
+ template <typename TProcess1, typename TProcess2, int IndexStart, int IndexProcess1,
+ int IndexProcess2>
+ template <typename TStack>
+ void ProcessSequence<TProcess1, TProcess2, IndexStart, IndexProcess1,
+ IndexProcess2>::doCascadeEquations(TStack& stack) {
+
+ if constexpr (is_process_v<process1_type>) { // to protect from further compiler
+ // errors if process1_type is invalid
+ if constexpr (process1_type::is_process_sequence &&
+ !process1_type::is_switch_process_sequence) {
+
+ A_.doCascadeEquations(stack);
+
+ } else if constexpr (is_cascade_equations_process_v<process1_type>) {
+
+ // interface checking on TProcess1
+ static_assert(has_method_doCascadeEquations_v<TProcess1,
+ void, // return type
+ TStack&>, // parameter
+ "TDerived has no method with correct signature \"void "
+ "doCascadeEquations(TStack&)\" required for "
+ "CascadeEquationsProcess<TDerived>. ");
+
+ A_.doCascadeEquations(stack);
+ }
+ }
+
+ if constexpr (is_process_v<process2_type>) { // to protect from further compiler
+ // errors if process2_type is invalid
+ if constexpr (process2_type::is_process_sequence &&
+ !process2_type::is_switch_process_sequence) {
+
+ B_.doCascadeEquations(stack);
+
+ } else if constexpr (is_cascade_equations_process_v<process2_type>) {
+
+ // interface checking on TProcess2
+ static_assert(has_method_doCascadeEquations_v<TProcess2,
+ void, // return type
+ TStack&>, // parameter
+ "TDerived has no method with correct signature \"void "
+ "doCascadeEquations(TStack&)\" required for "
+ "CascadeEquationsProcess<TDerived>. ");
+
+ B_.doCascadeEquations(stack);
+ }
+ }
+ }
+
+ template <typename TProcess1, typename TProcess2, int IndexStart, int IndexProcess1,
+ int IndexProcess2>
+ void ProcessSequence<TProcess1, TProcess2, IndexStart, IndexProcess1,
+ IndexProcess2>::initCascadeEquations() {
+
+ if constexpr (is_process_v<process1_type>) { // to protect from further compiler
+ // errors if process1_type is invalid
+ if constexpr ((process1_type::is_process_sequence &&
+ !process1_type::is_switch_process_sequence) ||
+ is_cascade_equations_process_v<process1_type>) {
+ A_.initCascadeEquations();
+ }
+ }
+
+ if constexpr (is_process_v<process2_type>) { // to protect from further compiler
+ // errors if process2_type is invalid
+ if constexpr ((process2_type::is_process_sequence &&
+ !process2_type::is_switch_process_sequence) ||
+ is_cascade_equations_process_v<process2_type>) {
+ B_.initCascadeEquations();
+ }
+ }
+ } // namespace corsika
+
template <typename TProcess1, typename TProcess2, int IndexStart, int IndexProcess1,
int IndexProcess2>
template <typename TSecondaryView>
@@ -485,8 +559,8 @@ namespace corsika {
}
/**
- * traits marker to identify objects containing any StackProcesses
- **/
+ * traits marker to identify objects containing any StackProcesses.
+ */
namespace detail {
// need helper alias to achieve this:
template <
diff --git a/corsika/detail/framework/process/SecondariesProcess.hpp b/corsika/detail/framework/process/SecondariesProcess.hpp
index ee76cce853d29f7007d9a82ed88b12e0826a31f6..f0055581466dcfd0cfa60f040fe2323e4fefe51a 100644
--- a/corsika/detail/framework/process/SecondariesProcess.hpp
+++ b/corsika/detail/framework/process/SecondariesProcess.hpp
@@ -14,8 +14,8 @@
namespace corsika {
/**
- traits test for SecondariesProcess::doSecondaries method
- */
+ * traits test for SecondariesProcess::doSecondaries method.
+ */
template <class TProcess, typename TReturn, typename... TArg>
struct has_method_doSecondaries
: public detail::has_method_signature<TReturn, TArg...> {
@@ -38,16 +38,19 @@ namespace corsika {
public:
/**
- @name traits results
- @{
- */
+ * @name traits results
+ * @{
+ */
using type = decltype(test<std::decay_t<TProcess>>(nullptr));
static const bool value = type::value;
//! @}
};
- //! @file DecayProcess.hpp
- //! value traits type
+ /**
+ * @file SecondariesProcess.hpp
+ *
+ * @brief value traits type.
+ */
template <class TProcess, typename TReturn, typename... TArg>
bool constexpr has_method_doSecondaries_v =
has_method_doSecondaries<TProcess, TReturn, TArg...>::value;
diff --git a/corsika/detail/modules/StackInspector.inl b/corsika/detail/modules/StackInspector.inl
index 24a356fcb749f90796db5920bef351dee2d10ac9..06705ac8f6be0892750e8c3942e4710d273b6dc8 100644
--- a/corsika/detail/modules/StackInspector.inl
+++ b/corsika/detail/modules/StackInspector.inl
@@ -61,28 +61,44 @@ namespace corsika {
if (dE < dE_threshold_) return;
double const progress = dE / E0_;
- double const eta_seconds = elapsed_seconds.count() / progress;
std::time_t const start_time = std::chrono::system_clock::to_time_t(StartTime_);
- std::time_t const eta_time = std::chrono::system_clock::to_time_t(
- StartTime_ + std::chrono::seconds((int)eta_seconds));
-
- int const yday0 = std::localtime(&start_time)->tm_yday;
- int const yday1 = std::localtime(&eta_time)->tm_yday;
- int const dyday = yday1 - yday0;
-
- CORSIKA_LOG_INFO(
- "StackInspector: "
- " time={}"
- ", running={} seconds"
- " ( {:.1f}%)"
- ", nStep={}"
- ", stackSize={}"
- ", Estack={} GeV"
- ", ETA={}{}",
- std::put_time(std::localtime(&now_time), "%T"), elapsed_seconds.count(),
- (progress * 100), getStep(), vS.getSize(), Etot / 1_GeV,
- (dyday == 0 ? "" : fmt::format("+{}d ", dyday)),
- std::put_time(std::localtime(&eta_time), "%T"));
+
+ if (progress > 0) {
+
+ double const eta_seconds = elapsed_seconds.count() / progress;
+ std::time_t const eta_time = std::chrono::system_clock::to_time_t(
+ StartTime_ + std::chrono::seconds((int)eta_seconds));
+
+ int const yday0 = std::localtime(&start_time)->tm_yday;
+ int const yday1 = std::localtime(&eta_time)->tm_yday;
+ int const dyday = yday1 - yday0;
+
+ CORSIKA_LOG_INFO(
+ "StackInspector: "
+ " time={}"
+ ", running={} seconds"
+ " ( {:.1f}%)"
+ ", nStep={}"
+ ", stackSize={}"
+ ", Estack={} GeV"
+ ", ETA={}{}",
+ std::put_time(std::localtime(&now_time), "%T"), elapsed_seconds.count(),
+ (progress * 100), getStep(), vS.getSize(), Etot / 1_GeV,
+ (dyday == 0 ? "" : fmt::format("+{}d ", dyday)),
+ std::put_time(std::localtime(&eta_time), "%T"));
+ } else {
+ CORSIKA_LOG_INFO(
+ "StackInspector: "
+ " time={}"
+ ", running={} seconds"
+ " ( {:.1f}%)"
+ ", nStep={}"
+ ", stackSize={}"
+ ", Estack={} GeV"
+ ", ETA={}{}",
+ std::put_time(std::localtime(&now_time), "%T"), elapsed_seconds.count(),
+ (progress * 100), getStep(), vS.getSize(), Etot / 1_GeV, "---", "---");
+ }
}
} // namespace corsika
diff --git a/corsika/detail/modules/conex/CONEXhybrid.inl b/corsika/detail/modules/conex/CONEXhybrid.inl
index 3a8109f631d41275567bda02e713e3d47b404a8f..00e79005e236ef415bbdc5174d0615dec90f5049 100644
--- a/corsika/detail/modules/conex/CONEXhybrid.inl
+++ b/corsika/detail/modules/conex/CONEXhybrid.inl
@@ -28,6 +28,9 @@ namespace corsika {
: center_{center}
, showerAxis_{showerAxis}
, groundDist_{groundDist}
+ , injectionHeight_{injectionHeight}
+ , primaryEnergy_{primaryEnergy}
+ , primaryPDG_{primaryPDG}
, showerCore_{showerAxis_.getStart() + showerAxis_.getDirection() * groundDist_}
, conexObservationCS_{std::invoke([&]() {
auto const& c8cs = center.getCoordinateSystem();
@@ -77,7 +80,9 @@ namespace corsika {
CORSIKA_LOG_DEBUG("showerCore (C8): {}",
showerCore_.getCoordinates(center.getCoordinateSystem()));
- int randomSeeds[3] = {1234, 0, 0}; // will be overwritten later??
+ int randomSeeds[3] = {1234, 0,
+ 0}; // SEEDS ARE NOT USED. All random numbers are obtained from
+ // the CORSIKA 8 stream "conex" and "epos"!
int heModel = eSibyll23;
int nShower = 1; // large to avoid final stats.
@@ -92,8 +97,9 @@ namespace corsika {
particleListMode,
#endif
configPath.c_str(), configPath.size());
+ }
- double eprima = primaryEnergy / 1_GeV;
+ inline void CONEXhybrid::initCascadeEquations() {
// set phi, theta
Vector<length_d> ez{conexObservationCS_, {0._m, 0._m, -1_m}};
@@ -111,16 +117,17 @@ namespace corsika {
"; phi (deg) = {}",
theta, phi);
- int ipart = static_cast<int>(primaryPDG);
- auto rng = RNGManager<>::getInstance().getRandomStream("conex");
+ int ipart = static_cast<int>(primaryPDG_);
double dimpact = 0.; // valid only if shower core is fixed on the observation plane;
// for skimming showers an offset is needed like in CONEX
- std::array<int, 3> ioseed{static_cast<int>(rng()), static_cast<int>(rng()),
- static_cast<int>(rng())};
+ // SEEDS ARE NOT USED. All random numbers are obtained from
+ // the CORSIKA 8 stream "conex" and "epos"!
+ std::array<int, 3> ioseed{1, 1, 1};
- double xminp = injectionHeight / 1_m;
+ double eprima = primaryEnergy_ / 1_GeV;
+ double xminp = injectionHeight_ / 1_m;
::conex::conexrun_(ipart, eprima, theta, phi, xminp, dimpact, ioseed.data());
}
@@ -225,7 +232,8 @@ namespace corsika {
return true;
}
- inline void CONEXhybrid::solveCE() {
+ template <typename TStack>
+ inline void CONEXhybrid::doCascadeEquations(TStack&) {
::conex::conexcascade_();
diff --git a/corsika/detail/modules/sibyll/Interaction.inl b/corsika/detail/modules/sibyll/Interaction.inl
index e4fdeec751684a0f143fbcc4420a51132bea5cf4..ed42473ec865d339e540e7b7e1904df763136e09 100644
--- a/corsika/detail/modules/sibyll/Interaction.inl
+++ b/corsika/detail/modules/sibyll/Interaction.inl
@@ -206,10 +206,8 @@ namespace corsika::sibyll {
// just for show:
// boost projecticle
[[maybe_unused]] auto const PprojCoM = boost.toCoM(PprojLab);
-
// boost target
[[maybe_unused]] auto const PtargCoM = boost.toCoM(PtargLab);
-
CORSIKA_LOG_DEBUG(
"Interaction: ebeam CoM: {} GeV "
"Interaction: pbeam CoM: {} GeV ",
@@ -328,15 +326,16 @@ namespace corsika::sibyll {
Ecm_final += psib.getEnergy();
}
CORSIKA_LOG_DEBUG(
- "conservation (all GeV):"
- "Ecm_initial(per nucleon)={}, Ecm_final(per nucleon)={}, "
- "Elab_initial={}, Elab_final={}, "
- "diff (%)={}, "
- "E in nucleons={}, "
- "Plab_initial={}, "
- "Plab_final={} ",
+ "conservation (all GeV): "
+ "Ecm_initial(per nucleon)={:.2f}, Ecm_final(per nucleon)={:.2f}, "
+ "Elab_initial={:.2f}, Elab_final={:.2f}, "
+ "Elab-diff (%)={:.2f}, "
+ "m in target nucleons={:.2f}, "
+ "Plab_initial={:.2f}, "
+ "Plab_final={:.2f} ",
Ecm / 1_GeV, Ecm_final * 2. / (get_nwounded() + 1) / 1_GeV, Etot / 1_GeV,
- Elab_final / 1_GeV, (Elab_final / Etot / get_nwounded() - 1) * 100,
+ Elab_final / 1_GeV,
+ (Elab_final / (Etot + get_nwounded() * constants::nucleonMass) - 1) * 100,
constants::nucleonMass * get_nwounded() / 1_GeV,
(pProjectileLab / 1_GeV).getComponents(), (Plab_final / 1_GeV).getComponents());
}
diff --git a/corsika/detail/modules/tracking/Intersect.inl b/corsika/detail/modules/tracking/Intersect.inl
index 6888f09a0607673cb82ae87d58a5ed6c67db11ce..b81cbf5c2eca1202f8ee7bb05aaf9f292424b854 100644
--- a/corsika/detail/modules/tracking/Intersect.inl
+++ b/corsika/detail/modules/tracking/Intersect.inl
@@ -60,7 +60,8 @@ namespace corsika {
for (auto const& node : volumeNode.getChildNodes()) {
Intersections const time_intersections = TDerived::intersect(particle, *node);
- CORSIKA_LOG_DEBUG("intersection times with child volume {}", fmt::ptr(node));
+ CORSIKA_LOG_DEBUG("search intersection times with child volume {}:",
+ fmt::ptr(node));
if (!time_intersections.hasIntersections()) { continue; }
auto const t_entry = time_intersections.getEntry();
auto const t_exit = time_intersections.getExit();
diff --git a/corsika/framework/process/CascadeEquationsProcess.hpp b/corsika/framework/process/CascadeEquationsProcess.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..5e92389f6be3d99e62e705559eabc31e5a50e648
--- /dev/null
+++ b/corsika/framework/process/CascadeEquationsProcess.hpp
@@ -0,0 +1,55 @@
+/*
+ * (c) Copyright 2020 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/process/BaseProcess.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+
+#include <corsika/detail/framework/process/CascadeEquationsProcess.hpp> // for extra traits, method/interface checking
+
+namespace corsika {
+
+ /**
+ * @ingroup Processes
+ * @{
+ *
+ * Processes executing cascade-equations calculations.
+ *
+ * Create a new CascadeEquationsProcess, e.g. for XYModel, via:
+ * @code{.cpp}
+ * class XYModel : public CascadeEquationsProcess<XYModel> {};
+ * @endcode
+ *
+ * and provide the necessary interface method:
+ * @code{.cpp}
+ * template <typename TStack>
+ * void doCascadeEquations(TStack& stack);
+ * @endcode
+ *
+ * Cascade equation processes may generate new particles on the stack. They also
+ * typically generate output.
+ */
+
+ template <typename TDerived>
+ class CascadeEquationsProcess : public BaseProcess<TDerived> {
+ public:
+ };
+
+ /**
+ * ProcessTraits specialization to flag CascadeEquationsProcess objects.
+ */
+ template <typename TProcess>
+ struct is_cascade_equations_process<
+ TProcess, std::enable_if_t<std::is_base_of_v<
+ CascadeEquationsProcess<typename std::decay_t<TProcess>>,
+ typename std::decay_t<TProcess>>>> : std::true_type {};
+
+ //! @}
+
+} // namespace corsika
diff --git a/corsika/framework/process/ProcessSequence.hpp b/corsika/framework/process/ProcessSequence.hpp
index 05191c248d5e1b17ceb5cfe39fc5dfe2fcabe4f6..407f403495c7c4d281880a897366f3b76b22a123 100644
--- a/corsika/framework/process/ProcessSequence.hpp
+++ b/corsika/framework/process/ProcessSequence.hpp
@@ -20,6 +20,7 @@
#include <corsika/framework/process/ContinuousProcessStepLength.hpp>
#include <corsika/framework/process/DecayProcess.hpp>
#include <corsika/framework/process/InteractionProcess.hpp>
+#include <corsika/framework/process/CascadeEquationsProcess.hpp>
#include <corsika/framework/process/ProcessReturn.hpp>
#include <corsika/framework/process/SecondariesProcess.hpp>
#include <corsika/framework/process/StackProcess.hpp>
@@ -209,6 +210,17 @@ namespace corsika {
template <typename TStack>
void doStack(TStack& stack);
+ /**
+ * Execute the CascadeEquationsProcess-es in the ProcessSequence.
+ */
+ template <typename TStack>
+ void doCascadeEquations(TStack& stack);
+
+ /**
+ * Init the CascadeEquationsProcess-es in the ProcessSequence.
+ */
+ void initCascadeEquations();
+
/**
* Calculate the maximum allowed length of the next tracking step, based on all
* ContinuousProcess-es
diff --git a/corsika/framework/process/ProcessTraits.hpp b/corsika/framework/process/ProcessTraits.hpp
index cccd80abaf1a11436befc913eaaf785112c5b123..24910fece0f821a7ac21b4e271ac0bb2c72f8140 100644
--- a/corsika/framework/process/ProcessTraits.hpp
+++ b/corsika/framework/process/ProcessTraits.hpp
@@ -17,7 +17,7 @@
namespace corsika {
/**
- * A traits marker to identify BaseProcess, thus any type of process
+ * A traits marker to identify BaseProcess, thus any type of process.
*/
template <typename TProcess, typename TEnable = void>
struct is_process : std::false_type {};
@@ -26,7 +26,7 @@ namespace corsika {
bool constexpr is_process_v = is_process<TProcess>::value;
/**
- * A traits marker to identify ContinuousProcess
+ * A traits marker to identify ContinuousProcess.
*/
template <typename TProcess, typename Enable = void>
struct is_continuous_process : std::false_type {};
@@ -35,7 +35,7 @@ namespace corsika {
bool constexpr is_continuous_process_v = is_continuous_process<TProcess>::value;
/**
- * A traits marker to identify DecayProcess
+ * A traits marker to identify DecayProcess.
*/
template <typename TProcess, typename Enable = void>
struct is_decay_process : std::false_type {};
@@ -44,7 +44,7 @@ namespace corsika {
bool constexpr is_decay_process_v = is_decay_process<TProcess>::value;
/**
- * A traits marker to identify StackProcess
+ * A traits marker to identify StackProcess.
*/
template <typename TProcess, typename Enable = void>
struct is_stack_process : std::false_type {};
@@ -53,7 +53,17 @@ namespace corsika {
bool constexpr is_stack_process_v = is_stack_process<TProcess>::value;
/**
- * A traits marker to identify SecondariesProcess
+ * A traits marker to identify CascadeEquationsProcess.
+ */
+ template <typename TProcess, typename Enable = void>
+ struct is_cascade_equations_process : std::false_type {};
+
+ template <typename TProcess>
+ bool constexpr is_cascade_equations_process_v =
+ is_cascade_equations_process<TProcess>::value;
+
+ /**
+ * A traits marker to identify SecondariesProcess.
*/
template <typename TProcess, typename Enable = void>
struct is_secondaries_process : std::false_type {};
@@ -62,7 +72,7 @@ namespace corsika {
bool constexpr is_secondaries_process_v = is_secondaries_process<TProcess>::value;
/**
- * A traits marker to identify BoundaryProcess
+ * A traits marker to identify BoundaryProcess.
*/
template <typename TProcess, typename Enable = void>
struct is_boundary_process : std::false_type {};
@@ -71,7 +81,7 @@ namespace corsika {
bool constexpr is_boundary_process_v = is_boundary_process<TProcess>::value;
/**
- * A traits marker to identify InteractionProcess
+ * A traits marker to identify InteractionProcess.
*/
template <typename TProcess, typename Enable = void>
struct is_interaction_process : std::false_type {};
@@ -80,8 +90,8 @@ namespace corsika {
bool constexpr is_interaction_process_v = is_interaction_process<TProcess>::value;
/**
- * A traits marker to identify ProcessSequence that contain a StackProcess
- **/
+ * A traits marker to identify ProcessSequence that contain a StackProcess.
+ */
template <typename TClass>
struct contains_stack_process : std::false_type {};
@@ -89,8 +99,8 @@ namespace corsika {
bool constexpr contains_stack_process_v = contains_stack_process<TClass>::value;
/**
- * traits class to count any type of Process, general version
- **/
+ * 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;
diff --git a/corsika/framework/process/SecondariesProcess.hpp b/corsika/framework/process/SecondariesProcess.hpp
index 975b1e7851174b3be25e43a99ef06c5155a2912c..8fdc6b67d9f5e2c720b59773efcb7131fa57c886 100644
--- a/corsika/framework/process/SecondariesProcess.hpp
+++ b/corsika/framework/process/SecondariesProcess.hpp
@@ -16,24 +16,24 @@
namespace corsika {
/**
- @ingroup Processes
- @{
-
- Processes acting on the secondaries produced by other processes.
-
- Create a new SecondariesProcess, e.g. for XYModel, via
- @code{.cpp}
- class XYModel : public SecondariesProcess<XYModel> {};
- @endcode
-
- and provide the necessary interface method:
- @code{.cpp}
- template <typename TStackView>
- void doSecondaries(TStackView& StackView);
- @endcode
-
- where StackView is an object that can store secondaries on a
- stack and also iterate over these secondaries.
+ * @ingroup Processes
+ * @{
+ *
+ * Processes acting on the secondaries produced by other processes.
+ *
+ * Create a new SecondariesProcess, e.g. for XYModel, via
+ * @code{.cpp}
+ * class XYModel : public SecondariesProcess<XYModel> {};
+ * @endcode
+ *
+ * and provide the necessary interface method:
+ * @code{.cpp}
+ * template <typename TStackView>
+ * void doSecondaries(TStackView& StackView);
+ * @endcode
+ *
+ * where StackView is an object that can store secondaries on a
+ * stack and also iterate over these secondaries.
*/
template <typename TDerived>
@@ -42,8 +42,8 @@ namespace corsika {
};
/**
- * ProcessTraits specialization to flag SecondariesProcess objects
- **/
+ * ProcessTraits specialization to flag SecondariesProcess objects.
+ */
template <typename TProcess>
struct is_secondaries_process<
TProcess, std::enable_if_t<
diff --git a/corsika/framework/process/StackProcess.hpp b/corsika/framework/process/StackProcess.hpp
index 773265f49e102b63ca447191a25040c555179c37..34a9c50bceef5412a6b5f1b19cda29e02b03f8bd 100644
--- a/corsika/framework/process/StackProcess.hpp
+++ b/corsika/framework/process/StackProcess.hpp
@@ -16,32 +16,32 @@
namespace corsika {
/**
- @ingroup Processes
- @{
-
- Process to act on the entire particle stack
-
- Create a new StackProcess, e.g. for XYModel, via
- @code{.cpp}
- class XYModel : public StackProcess<XYModel> {};
- @endcode
-
- and provide the necessary interface method
- @code{.cpp}
- template <typename TStack>
- void XYModel::doStack(TStack&);
- @endcode
-
- Where, of course, Stack is the valid
- class to access particles on the Stack. This methods does
- not need to be templated, they could use the types
- e.g. corsika::setup::Stack directly -- but by the cost of
- loosing all flexibility otherwise provided.
-
- A StackProcess has only one constructor `StackProcess::StackProcess(unsigned int
- const nStep)` where nStep is the number of steps of the cascade stepping after which
- the stack process should be run. Good values are on the order of 1000, which will not
- compromise run time in the end, but provide all the benefits of the StackProcess.
+ * @ingroup Processes
+ * @{
+ *
+ * Process to act on the entire particle stack.
+ *
+ * Create a new StackProcess, e.g. for XYModel, via:
+ * @code{.cpp}
+ * class XYModel : public StackProcess<XYModel> {};
+ * @endcode
+ *
+ * and provide the necessary interface method:
+ * @code{.cpp}
+ * template <typename TStack>
+ * void XYModel::doStack(TStack&);
+ * @endcode
+ *
+ * Where, of course, Stack is the valid
+ * class to access particles on the Stack. This methods does
+ * not need to be templated, they could use the types
+ * e.g. corsika::setup::Stack directly -- but by the cost of
+ * loosing all flexibility otherwise provided.
+ *
+ * A StackProcess has only one constructor `StackProcess::StackProcess(unsigned int
+ * const nStep)` where nStep is the number of steps of the cascade stepping after which
+ * the stack process should be run. Good values are on the order of 1000, which will not
+ * compromise run time in the end, but provide all the benefits of the StackProcess.
*/
template <typename TDerived>
@@ -64,10 +64,10 @@ namespace corsika {
private:
/**
- @name The number of "steps" during the cascade processing after
- which this StackProcess is going to be executed. The logic is
- "iStep_ modulo nStep_"
- @{
+ * @name The number of "steps" during the cascade processing after
+ * which this StackProcess is going to be executed. The logic is
+ * "iStep_ modulo nStep_"
+ * @{
*/
unsigned int nStep_ = 0;
unsigned long int iStep_ = 0;
@@ -75,8 +75,8 @@ namespace corsika {
};
/**
- * ProcessTraits specialization to flag StackProcess objects
- **/
+ * ProcessTraits specialization to flag StackProcess objects.
+ */
template <typename TProcess>
struct is_stack_process<
TProcess,
diff --git a/corsika/modules/conex/CONEXhybrid.hpp b/corsika/modules/conex/CONEXhybrid.hpp
index 2f0b80873dfa4109202f1a01cb11ff571c988cfd..0d3878e474212cab10c6fefaf81ad8b23067d145 100644
--- a/corsika/modules/conex/CONEXhybrid.hpp
+++ b/corsika/modules/conex/CONEXhybrid.hpp
@@ -8,11 +8,13 @@
#pragma once
+#include <corsika/framework/process/SecondariesProcess.hpp>
+#include <corsika/framework/process/CascadeEquationsProcess.hpp>
+
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/core/ParticleProperties.hpp>
#include <corsika/framework/geometry/Point.hpp>
#include <corsika/framework/geometry/Vector.hpp>
-#include <corsika/framework/process/SecondariesProcess.hpp>
#include <corsika/media/ShowerAxis.hpp>
#include <corsika/modules/conex/CONEX_f.hpp>
@@ -23,17 +25,37 @@ namespace corsika {
LengthType constexpr earthRadius{6371315 * meter};
} // namespace conex
- class CONEXhybrid : public SecondariesProcess<CONEXhybrid> {
+ class CONEXhybrid : public CascadeEquationsProcess<CONEXhybrid>,
+ public SecondariesProcess<CONEXhybrid> {
public:
CONEXhybrid(Point center, ShowerAxis const& showerAxis, LengthType groundDist,
LengthType injectionHeight, HEPEnergyType primaryEnergy, PDGCode pdg);
+ /**
+ * Main entry point to pass new particle data towards CONEX. If a
+ * particles is selected for CONEX, it is removed from the CORSIKA
+ * 8 stack.
+ */
template <typename TStackView>
void doSecondaries(TStackView&);
- void solveCE();
+ /**
+ * init currently needs to be called to initializa a new
+ * event. All tables are cleared, etc.
+ */
+ void initCascadeEquations();
+
+ /**
+ * Cascade equations are solved basoned on the data in the tables
+ */
+ template <typename TStack>
+ void doCascadeEquations(TStack& stack);
+ /**
+ * Internal function to fill particle data inside CONEX
+ * tables. Only e.m. particles are selected right now.
+ */
bool addParticle(Code pid, HEPEnergyType energy, HEPEnergyType mass,
Point const& position, Vector<dimensionless_d> const& direction,
TimeType t);
@@ -51,8 +73,11 @@ namespace corsika {
Point const center_; //!< center of CONEX Earth
ShowerAxis const& showerAxis_;
- LengthType groundDist_; //!< length from injection point to shower core
- Point const showerCore_; //!< shower core
+ LengthType groundDist_; //!< length from injection point to shower core
+ LengthType injectionHeight_; //!< starting height of primary particle
+ HEPEnergyType primaryEnergy_; //!< primary particle energy
+ PDGCode primaryPDG_; //!< primary particle PDG
+ Point const showerCore_; //!< shower core
CoordinateSystemPtr const conexObservationCS_; //!< CONEX observation frame
DirectionVector const x_sf_,
y_sf_; //!< unit vectors of CONEX shower frame, z_sf is shower axis direction
diff --git a/corsika/modules/conex/Random.hpp b/corsika/modules/conex/Random.hpp
index a638a330857ebc2c011c82237df2dd4a1a1f08eb..9532eace8ed29d4322def04b5ad5fee64761b366 100644
--- a/corsika/modules/conex/Random.hpp
+++ b/corsika/modules/conex/Random.hpp
@@ -17,6 +17,19 @@
* This file is an integral part of the epos interface. It must be
* linked to the executable linked to epos exactly once
*
+ * Note, that the fortran random numbe interface functions are all
+ * defined in the epos corsika 8 interface:
+ *
+ * ranfst, ranfgt, rmmaqd, ranfini, ranfcv, rmmard, rangen, drangen
+ *
+ * All of them use the epos_random_interface registered as "epos" stream.
+ *
+ * Thus, the fortran part of CONEX will use the "epos" CORSIKA 8 random stream,
+ * only the CONEX c++ part will use the "conex" random stream.
+ *
+ * Since EPOS and CONEX use the same fortran symbols for access to
+ * random numbers this can only be changed by renaming inside the
+ * fortran part.
*/
namespace conex {
diff --git a/examples/hybrid_MC.cpp b/examples/hybrid_MC.cpp
index 04a126366a2d408e074532b447f3357aa4145bed..dabed46e9f4e598b4146eebd673cc9ad07721fd0 100644
--- a/examples/hybrid_MC.cpp
+++ b/examples/hybrid_MC.cpp
@@ -291,8 +291,6 @@ int main(int argc, char** argv) {
cut.reset();
eLoss.reset();
- conex_model.solveCE();
-
auto const hists = sibyllCounted.getHistogram() + sibyllNucCounted.getHistogram() +
urqmdCounted.getHistogram();
diff --git a/tests/modules/testCONEX.cpp b/tests/modules/testCONEX.cpp
index f86dd7b2d642c58278995823b9a250d27b36926c..985043d4dca501fab27359ff8be5eb557615ca8f 100644
--- a/tests/modules/testCONEX.cpp
+++ b/tests/modules/testCONEX.cpp
@@ -46,6 +46,8 @@ const std::string refDataDir = std::string(REFDATADIR); // from cmake
template <typename T>
using MExtraEnvirnoment = MediumPropertyModel<UniformMagneticField<T>>;
+struct DummyStack {};
+
TEST_CASE("CONEXSourceCut") {
logging::set_level(logging::level::info);
@@ -70,7 +72,7 @@ TEST_CASE("CONEXSourceCut") {
builder.setNuclearComposition(
{{Code::Nitrogen, Code::Oxygen},
- {0.7847f, 1.f - 0.7847f}}); // values taken from AIRES manual, Ar removed for now
+ {0.7847, 1. - 0.7847}}); // values taken from AIRES manual, Ar removed for now
builder.addExponentialLayer(1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 4_km);
builder.addExponentialLayer(1144.9069_g / (1_cm * 1_cm), 878153.55_cm, 10_km);
@@ -101,6 +103,7 @@ TEST_CASE("CONEXSourceCut") {
[[maybe_unused]] corsika::sibyll::NuclearInteraction sibyllNuc(sibyll, env);
CONEXhybrid conex(center, showerAxis, t, injectionHeight, E0, get_PDG(Code::Proton));
+ conex.initCascadeEquations();
HEPEnergyType const Eem{1_PeV};
auto const momentum = showerAxis.getDirection() * Eem;
@@ -122,7 +125,8 @@ TEST_CASE("CONEXSourceCut") {
auto const momentumPhoton = showerAxis.getDirection() * 1_TeV;
conex.addParticle(Code::Photon, 1_TeV, 0_eV, emPosition, momentumPhoton.normalized(),
0_s);
- conex.solveCE();
+ DummyStack stack;
+ conex.doCascadeEquations(stack);
}
#include <algorithm>