From 56ffcb680b3001015d004ec8ad93695fd6d03527 Mon Sep 17 00:00:00 2001
From: Alan Coleman <alanc@udel.edu>
Date: Fri, 8 Mar 2024 10:37:52 +0100
Subject: [PATCH] add forceDecay to Cascade
---
corsika/detail/framework/core/Cascade.inl | 29 +++++++++++++++++++++++
corsika/framework/core/Cascade.hpp | 10 ++++++++
2 files changed, 39 insertions(+)
diff --git a/corsika/detail/framework/core/Cascade.inl b/corsika/detail/framework/core/Cascade.inl
index bb741627f..2a8b32489 100644
--- a/corsika/detail/framework/core/Cascade.inl
+++ b/corsika/detail/framework/core/Cascade.inl
@@ -97,6 +97,19 @@ namespace corsika {
template <typename TTracking, typename TProcessList, typename TOutput, typename TStack>
inline void Cascade<TTracking, TProcessList, TOutput, TStack>::forceInteraction() {
forceInteraction_ = true;
+ if (forceDecay_) {
+ CORSIKA_LOG_ERROR("Cannot set forceInteraction and forceDecay");
+ throw std::runtime_error("Cannot set forceInteraction and forceDecay");
+ }
+ }
+
+ template <typename TTracking, typename TProcessList, typename TOutput, typename TStack>
+ inline void Cascade<TTracking, TProcessList, TOutput, TStack>::forceDecay() {
+ forceDecay_ = true;
+ if (forceInteraction_) {
+ CORSIKA_LOG_ERROR("Cannot set forceDecay and forceInteraction");
+ throw std::runtime_error("Cannot set forceDecay and forceInteraction");
+ }
}
template <typename TTracking, typename TProcessList, typename TOutput, typename TStack>
@@ -140,6 +153,22 @@ namespace corsika {
return;
}
+ if (forceDecay_) {
+ CORSIKA_LOG_TRACE("forced decay!");
+ forceDecay_ = false; // just one (first) interaction
+ stack_view_type secondaries(particle);
+ decay(secondaries, sequence_.getInverseLifetime(particle));
+ if (secondaries.getSize() == 1 && secondaries.getProjectile().getPID() ==
+ secondaries.getNextParticle().getPID()) {
+ throw std::runtime_error(
+ fmt::format("Particle {} decays into itself!",
+ get_name(secondaries.getProjectile().getPID())));
+ }
+ sequence_.doSecondaries(secondaries);
+ particle.erase(); // primary particle is done
+ return;
+ }
+
// calculate interaction length in medium
GrammageType const total_lambda =
(composition.getAverageMassNumber() * constants::u) / total_cx_pre;
diff --git a/corsika/framework/core/Cascade.hpp b/corsika/framework/core/Cascade.hpp
index d4b7d9774..87bea3142 100644
--- a/corsika/framework/core/Cascade.hpp
+++ b/corsika/framework/core/Cascade.hpp
@@ -92,9 +92,18 @@ namespace corsika {
* Force an interaction of the top particle of the stack at its current position.
* Note that setNodes() or an equivalent procedure needs to be called first if you
* want to call forceInteraction() for the primary interaction.
+ * Incompatible with forceDecay()
*/
void forceInteraction();
+ /**
+ * Force an decay of the top particle of the stack at its current position.
+ * Note that setNodes() or an equivalent procedure needs to be called first if you
+ * want to call forceDecay() for the primary interaction.
+ * Incompatible with forceInteraction()
+ */
+ void forceDecay();
+
private:
/**
* The Step function is executed for each particle from the
@@ -122,6 +131,7 @@ namespace corsika {
TStack& stack_;
default_prng_type& rng_ = RNGManager<>::getInstance().getRandomStream("cascade");
bool forceInteraction_;
+ bool forceDecay_;
unsigned int count_ = 0;
// but this here temporarily. Should go into dedicated file later:
--
GitLab