From f954a7973006088b3add05bef8305c8414d545e1 Mon Sep 17 00:00:00 2001
From: Alan Coleman <alanc@udel.edu>
Date: Mon, 11 Mar 2024 18:48:09 +0100
Subject: [PATCH] addressing comments by Felix
---
applications/c8_air_shower.cpp | 9 +++++++++
corsika/detail/framework/core/Cascade.inl | 10 +++++-----
corsika/detail/modules/writers/InteractionWriter.inl | 9 +++++++--
corsika/modules/writers/InteractionWriter.hpp | 3 ++-
python/corsika/io/outputs/interaction.py | 2 +-
5 files changed, 24 insertions(+), 9 deletions(-)
diff --git a/applications/c8_air_shower.cpp b/applications/c8_air_shower.cpp
index 473e05b86..7a1e19dbf 100644
--- a/applications/c8_air_shower.cpp
+++ b/applications/c8_air_shower.cpp
@@ -219,6 +219,10 @@ int main(int argc, char** argv) {
app.add_flag("--force-interaction", force_interaction,
"Force the location of the first interaction.")
->group("Misc.");
+ bool force_decay = false;
+ app.add_flag("--force-decay", force_decay,
+ "Force the primary to immediately decay")
+ ->group("Misc.");
bool disable_interaction_hists = false;
app.add_flag("--disable-interaction-histograms", disable_interaction_hists,
"Store interaction histograms")
@@ -632,6 +636,11 @@ int main(int argc, char** argv) {
EAS.forceInteraction();
}
+ if (force_decay) {
+ CORSIKA_LOG_INFO("Forcing the primary to decay");
+ EAS.forceDecay();
+ }
+
primaryWriter.recordPrimary(primaryProperties);
// run the shower
diff --git a/corsika/detail/framework/core/Cascade.inl b/corsika/detail/framework/core/Cascade.inl
index 5819be9e6..6498cc520 100644
--- a/corsika/detail/framework/core/Cascade.inl
+++ b/corsika/detail/framework/core/Cascade.inl
@@ -99,8 +99,8 @@ namespace corsika {
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");
+ CORSIKA_LOG_ERROR("Cannot set forceInteraction when forceDecay is already set");
+ throw std::runtime_error("Cannot set forceInteraction when forceDecay is already set");
}
}
@@ -108,8 +108,8 @@ namespace corsika {
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");
+ CORSIKA_LOG_ERROR("Cannot set forceDecay when forceInteraction is already set");
+ throw std::runtime_error("Cannot set forceDecay when forceInteraction is already set");
}
}
@@ -156,7 +156,7 @@ namespace corsika {
if (forceDecay_) {
CORSIKA_LOG_TRACE("forced decay!");
- forceDecay_ = false; // just one (first) interaction
+ forceDecay_ = false; // just one decay
stack_view_type secondaries(particle);
decay(secondaries, sequence_.getInverseLifetime(particle));
if (secondaries.getSize() == 1 && secondaries.getProjectile().getPID() ==
diff --git a/corsika/detail/modules/writers/InteractionWriter.inl b/corsika/detail/modules/writers/InteractionWriter.inl
index 1e0108dbe..b1cfd4d4e 100644
--- a/corsika/detail/modules/writers/InteractionWriter.inl
+++ b/corsika/detail/modules/writers/InteractionWriter.inl
@@ -24,13 +24,17 @@ namespace corsika {
template <typename TTracking, typename TOutput>
template <typename TStackView>
inline void InteractionWriter<TTracking, TOutput>::doSecondaries(TStackView& vS) {
- if (interactionCounter_++) { return; }
+ // Dumps the stack to the output parquet stream
+ // The primary and secondaries are all written along with the slant depth
+
+ if (interactionCounter_++) { return; } // Only run on the first interaction
auto primary = vS.getProjectile();
auto dX = showerAxis_.getProjectedX(primary.getPosition());
CORSIKA_LOG_INFO("First interaction at dX {}", dX);
CORSIKA_LOG_INFO("Primary: {}, E {}", primary.getPID(), primary.getKineticEnergy());
+ // Get the displacement of the primary particle w.r.t. observation plane
Vector const dr_prim = primary.getPosition() - obsPlane_.getPlane().getCenter();
auto const p_prim = primary.getMomentum();
@@ -49,6 +53,7 @@ namespace corsika {
auto particle = vS.begin();
while (particle != vS.end()) {
+ // Get the displacement of the secondary w.r.t. observation plane
Vector const dr_2nd = particle.getPosition() - obsPlane_.getPlane().getCenter();
auto const p_2nd = particle.getMomentum();
@@ -63,7 +68,7 @@ namespace corsika {
<< static_cast<double>(particle.getTime() / 1_ns)
<< static_cast<float>(dX / (1_g / 1_cm / 1_cm)) << parquet::EndRow;
- CORSIKA_LOG_INFO(" 2ndary: {}, E {}", particle.getPID(),
+ CORSIKA_LOG_INFO(" 2ndary: {}, E_kin {}", particle.getPID(),
particle.getKineticEnergy());
++particle;
++nSecondaries_;
diff --git a/corsika/modules/writers/InteractionWriter.hpp b/corsika/modules/writers/InteractionWriter.hpp
index 6515ad195..47179405d 100644
--- a/corsika/modules/writers/InteractionWriter.hpp
+++ b/corsika/modules/writers/InteractionWriter.hpp
@@ -27,7 +27,8 @@ namespace corsika {
public:
/**
*
- * @param axis - axis used for calculating grammage
+ * @param axis - axis used for calculating slant depth of interaction
+ * @param obsPlane - used to define the location of the particles
*/
InteractionWriter(ShowerAxis const& axis,
ObservationPlane<TTracking, TOutput> const& obsPlane);
diff --git a/python/corsika/io/outputs/interaction.py b/python/corsika/io/outputs/interaction.py
index 64a010319..3a3d83da3 100644
--- a/python/corsika/io/outputs/interaction.py
+++ b/python/corsika/io/outputs/interaction.py
@@ -1,5 +1,5 @@
"""
- Read data written by EnergyLoss.
+ Read data written by InteractionWriter.
(c) Copyright 2024 CORSIKA Project, corsika-project@lists.kit.edu
--
GitLab