diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 918d02cf4aa7b20b0645be8d79b750667d14c023..e02c515945b6e6f25573ba06b153509301b986de 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -222,7 +222,7 @@ build_test-clang-14:
junit:
- ${CI_PROJECT_DIR}/build/test_outputs/junit*.xml
paths:
- - ${CI_PROJECT_DIR}/build/build_examples/example_outputs/radio_em_shower_outputs #python examples need this
+ - ${CI_PROJECT_DIR}/build/build_examples/example_outputs/radio_em_shower_outputs #python tests/examples need this
- ${CI_PROJECT_DIR}/build/test_outputs/junit*.xml
- ${CI_PROJECT_DIR}/build/CMakeCache.txt #python tests need this
- ${CI_PROJECT_DIR}/build/bin #python tests need this
@@ -411,9 +411,7 @@ python-tests:
allow_failure: false
python-examples:
- stage: python
- tags:
- - corsika
+ extends: .python
image: corsika/devel:u-22.04
needs:
- job: build_test_example-u-22_04
diff --git a/applications/c8_air_shower.cpp b/applications/c8_air_shower.cpp
index 7abd71bbb8240b370296ac4e99ab91e8738b11c9..8912d6a91529770bfb2390a1ccbf4cfd90df9635 100644
--- a/applications/c8_air_shower.cpp
+++ b/applications/c8_air_shower.cpp
@@ -27,6 +27,7 @@
#include <corsika/framework/utility/SaveBoostHistogram.hpp>
#include <corsika/modules/writers/EnergyLossWriter.hpp>
+#include <corsika/modules/writers/InteractionWriter.hpp>
#include <corsika/modules/writers/LongitudinalWriter.hpp>
#include <corsika/modules/writers/PrimaryWriter.hpp>
#include <corsika/modules/writers/SubWriter.hpp>
@@ -218,6 +219,9 @@ 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")
@@ -571,10 +575,15 @@ int main(int argc, char** argv) {
// register ZHS with the output manager
output.add("ZHS", zhs);
+ // make and register the first interaction writer
+ InteractionWriter<setup::Tracking, ParticleWriterParquet> inter_writer(
+ showerAxis, observationLevel);
+ output.add("interactions", inter_writer);
+
// assemble the final process sequence with radio
auto sequence = make_sequence(stackInspect, neutrinoPrimaryPythia, hadronSequence,
decayPythia, emCascade, emContinuous, coreas, zhs,
- longprof, observationLevel, thinning, cut);
+ longprof, observationLevel, inter_writer, thinning, cut);
/* === END: SETUP PROCESS LIST === */
@@ -627,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 bb741627f9cdebe233694a7c55ea5b979a47c46b..63ccdd7dd19b939e75872353618452876523d3da 100644
--- a/corsika/detail/framework/core/Cascade.inl
+++ b/corsika/detail/framework/core/Cascade.inl
@@ -43,7 +43,8 @@ namespace corsika {
, sequence_(pl)
, output_(out)
, stack_(stack)
- , forceInteraction_(false) {
+ , forceInteraction_(false)
+ , forceDecay_(false) {
CORSIKA_LOG_INFO(c8_ascii_);
CORSIKA_LOG_INFO("This is CORSIKA {}.{}.{}.{}", CORSIKA_RELEASE_NUMBER,
CORSIKA_MAJOR_NUMBER, CORSIKA_MINOR_NUMBER, CORSIKA_PATCH_NUMBER);
@@ -97,6 +98,21 @@ 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 when forceDecay is already set");
+ throw std::runtime_error(
+ "Cannot set forceInteraction when forceDecay is already set");
+ }
+ }
+
+ 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 when forceInteraction is already set");
+ throw std::runtime_error(
+ "Cannot set forceDecay when forceInteraction is already set");
+ }
}
template <typename TTracking, typename TProcessList, typename TOutput, typename TStack>
@@ -140,6 +156,22 @@ namespace corsika {
return;
}
+ if (forceDecay_) {
+ CORSIKA_LOG_TRACE("forced decay!");
+ forceDecay_ = false; // just one decay
+ 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/detail/modules/writers/InteractionWriter.inl b/corsika/detail/modules/writers/InteractionWriter.inl
new file mode 100644
index 0000000000000000000000000000000000000000..e4d23bd50f587f25454c73f574175a988f3d07e5
--- /dev/null
+++ b/corsika/detail/modules/writers/InteractionWriter.inl
@@ -0,0 +1,152 @@
+/*
+ * (c) Copyright 2024 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/Logging.hpp>
+
+namespace corsika {
+
+ template <typename TTracking, typename TOutput>
+ inline InteractionWriter<TTracking, TOutput>::InteractionWriter(
+ ShowerAxis const& axis, ObservationPlane<TTracking, TOutput> const& obsPlane)
+ : obsPlane_(obsPlane)
+ , showerAxis_(axis)
+ , interactionCounter_(0)
+ , showerId_(0) {}
+
+ template <typename TTracking, typename TOutput>
+ template <typename TStackView>
+ inline void InteractionWriter<TTracking, TOutput>::doSecondaries(TStackView& vS) {
+ // 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_kin {}", primary.getPID(),
+ primary.getKineticEnergy());
+
+ // get the location of the primary w.r.t. observation plane
+ Vector const displacement = primary.getPosition() - obsPlane_.getPlane().getCenter();
+
+ auto const x = displacement.dot(obsPlane_.getXAxis());
+ auto const y = displacement.dot(obsPlane_.getYAxis());
+ auto const z = displacement.dot(obsPlane_.getPlane().getNormal());
+ auto const nx = primary.getDirection().dot(obsPlane_.getXAxis());
+ auto const ny = primary.getDirection().dot(obsPlane_.getYAxis());
+ auto const nz = primary.getDirection().dot(obsPlane_.getPlane().getNormal());
+
+ auto const px = primary.getMomentum().dot(obsPlane_.getXAxis());
+ auto const py = primary.getMomentum().dot(obsPlane_.getYAxis());
+ auto const pz = primary.getMomentum().dot(obsPlane_.getPlane().getNormal());
+
+ summary_["shower_" + std::to_string(showerId_)]["pdg"] =
+ static_cast<int>(get_PDG(primary.getPID()));
+ summary_["shower_" + std::to_string(showerId_)]["name"] =
+ static_cast<std::string>(get_name(primary.getPID()));
+ summary_["shower_" + std::to_string(showerId_)]["total_energy"] =
+ (primary.getKineticEnergy() + get_mass(primary.getPID())) / 1_GeV;
+ summary_["shower_" + std::to_string(showerId_)]["kinetic_energy"] =
+ primary.getKineticEnergy() / 1_GeV;
+ summary_["shower_" + std::to_string(showerId_)]["x"] = x / 1_m;
+ summary_["shower_" + std::to_string(showerId_)]["y"] = y / 1_m;
+ summary_["shower_" + std::to_string(showerId_)]["z"] = z / 1_m;
+ summary_["shower_" + std::to_string(showerId_)]["nx"] = static_cast<double>(nx);
+ summary_["shower_" + std::to_string(showerId_)]["ny"] = static_cast<double>(ny);
+ summary_["shower_" + std::to_string(showerId_)]["nz"] = static_cast<double>(nz);
+ summary_["shower_" + std::to_string(showerId_)]["px"] =
+ static_cast<double>(px / 1_GeV);
+ summary_["shower_" + std::to_string(showerId_)]["py"] =
+ static_cast<double>(py / 1_GeV);
+ summary_["shower_" + std::to_string(showerId_)]["pz"] =
+ static_cast<double>(pz / 1_GeV);
+ summary_["shower_" + std::to_string(showerId_)]["time"] = primary.getTime() / 1_s;
+ summary_["shower_" + std::to_string(showerId_)]["slant_depth"] =
+ dX / (1_g / 1_cm / 1_cm);
+
+ uint nSecondaries = 0;
+
+ // Loop through secondaries
+ auto particle = vS.begin();
+ while (particle != vS.end()) {
+
+ // Get the momentum of the secondary, write w.r.t. observation plane
+ auto const p_2nd = particle.getMomentum();
+
+ *(output_.getWriter())
+ << showerId_ << static_cast<int>(get_PDG(particle.getPID()))
+ << static_cast<float>(p_2nd.dot(obsPlane_.getXAxis()) / 1_GeV)
+ << static_cast<float>(p_2nd.dot(obsPlane_.getYAxis()) / 1_GeV)
+ << static_cast<float>(p_2nd.dot(obsPlane_.getPlane().getNormal()) / 1_GeV)
+ << parquet::EndRow;
+
+ CORSIKA_LOG_INFO(" 2ndary: {}, E_kin {}", particle.getPID(),
+ particle.getKineticEnergy());
+ ++particle;
+ ++nSecondaries;
+ }
+
+ summary_["shower_" + std::to_string(showerId_)]["n_secondaries"] = nSecondaries;
+ }
+
+ template <typename TTracking, typename TOutput>
+ inline void InteractionWriter<TTracking, TOutput>::startOfLibrary(
+ boost::filesystem::path const& directory) {
+ output_.initStreamer((directory / ("interactions.parquet")).string());
+
+ output_.addField("pdg", parquet::Repetition::REQUIRED, parquet::Type::INT32,
+ parquet::ConvertedType::INT_32);
+ output_.addField("px", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
+ parquet::ConvertedType::NONE);
+ output_.addField("py", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
+ parquet::ConvertedType::NONE);
+ output_.addField("pz", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
+ parquet::ConvertedType::NONE);
+
+ output_.buildStreamer();
+
+ showerId_ = 0;
+ interactionCounter_ = 0;
+ summary_ = YAML::Node();
+ }
+
+ template <typename TTracking, typename TOutput>
+ inline void InteractionWriter<TTracking, TOutput>::startOfShower(
+ unsigned int const showerId) {
+ showerId_ = showerId;
+ interactionCounter_ = 0;
+ }
+
+ template <typename TTracking, typename TOutput>
+ inline void InteractionWriter<TTracking, TOutput>::endOfShower(unsigned int const) {}
+
+ template <typename TTracking, typename TOutput>
+ inline void InteractionWriter<TTracking, TOutput>::endOfLibrary() {
+ output_.closeStreamer();
+ }
+
+ template <typename TTracking, typename TOutput>
+ inline YAML::Node InteractionWriter<TTracking, TOutput>::getConfig() const {
+ YAML::Node node;
+ node["type"] = "Interactions";
+ node["units"]["energy"] = "GeV";
+ node["units"]["length"] = "m";
+ node["units"]["time"] = "ns";
+ node["units"]["grammage"] = "g/cm^2";
+ return node;
+ }
+
+ template <typename TTracking, typename TOutput>
+ inline YAML::Node InteractionWriter<TTracking, TOutput>::getSummary() const {
+ return summary_;
+ }
+
+} // namespace corsika
diff --git a/corsika/framework/core/Cascade.hpp b/corsika/framework/core/Cascade.hpp
index d4b7d9774c7f951dedb33db7df5c36e4daaf9762..87bea31427df998270957159e59d87e4b4a5ccf6 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:
diff --git a/corsika/modules/writers/EnergyLossWriter.hpp b/corsika/modules/writers/EnergyLossWriter.hpp
index 6f42de944526bcc32b6f33f73c83401c8e66df09..d0a201062dffcab7e458219b12b28f6434c5a61d 100644
--- a/corsika/modules/writers/EnergyLossWriter.hpp
+++ b/corsika/modules/writers/EnergyLossWriter.hpp
@@ -161,7 +161,7 @@ namespace corsika {
/**
* Return a summary.
*/
- YAML::Node getSummary() const;
+ YAML::Node getSummary() const override;
/**
* Return the configuration of this output.
diff --git a/corsika/modules/writers/InteractionWriter.hpp b/corsika/modules/writers/InteractionWriter.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..86527001f7a169df3924e200c7aee48ca4843b95
--- /dev/null
+++ b/corsika/modules/writers/InteractionWriter.hpp
@@ -0,0 +1,68 @@
+/*
+ * (c) Copyright 2024 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/SecondariesProcess.hpp>
+
+#include <corsika/media/ShowerAxis.hpp>
+
+#include <corsika/modules/ObservationPlane.hpp>
+
+#include <corsika/output/BaseOutput.hpp>
+#include <corsika/output/ParquetStreamer.hpp>
+
+namespace corsika {
+
+ template <typename TTracking, typename TOutput>
+ class InteractionWriter
+ : public SecondariesProcess<InteractionWriter<TTracking, TOutput>>,
+ public BaseOutput {
+
+ public:
+ /**
+ *
+ * @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);
+
+ /**
+ *
+ * @tparam TStackView
+ */
+ template <typename TStackView>
+ void doSecondaries(TStackView&);
+
+ void startOfLibrary(boost::filesystem::path const&) override;
+ void endOfLibrary() override;
+
+ void startOfShower(unsigned int const) override;
+ void endOfShower(unsigned int const) override;
+
+ YAML::Node getConfig() const override;
+ YAML::Node getSummary() const override;
+
+ auto getInteractionCounter() { return interactionCounter_; }
+
+ private:
+ ObservationPlane<TTracking, TOutput> const obsPlane_;
+ ShowerAxis const& showerAxis_;
+
+ unsigned int interactionCounter_;
+ unsigned int showerId_;
+
+ ParquetStreamer output_;
+ YAML::Node summary_;
+
+ }; // namespace corsika
+
+} // namespace corsika
+
+#include <corsika/detail/modules/writers/InteractionWriter.inl>
diff --git a/corsika/output/OutputManager.hpp b/corsika/output/OutputManager.hpp
index 6befc54ff741c55e63030237bdff2542d0312b81..1359b0d959322460ea9ab3f49621d74e9cfac6a1 100644
--- a/corsika/output/OutputManager.hpp
+++ b/corsika/output/OutputManager.hpp
@@ -35,8 +35,10 @@ namespace corsika {
/**
* Construct an OutputManager instance with a name in a given directory.
*
- * @param name The name of this output collection.
- * @param dir The directory where the output directory will be stored.
+ * @param name The name of this output collection.
+ * @param vseed The seed for the simulation (written to summary file)
+ * @param input_args The command line arguments at runtime (written to summary file)
+ * @param dir The directory where the output directory will be stored.
*/
OutputManager(std::string const& name, const long& vseed,
std::string const& input_args, boost::filesystem::path const& dir);
diff --git a/python/corsika/io/outputs/__init__.py b/python/corsika/io/outputs/__init__.py
index 54731dd42d775073d5d15e3d21e914706b5a7755..7ea762ac7fbcd666d5de961c99e0785bb9acd699 100644
--- a/python/corsika/io/outputs/__init__.py
+++ b/python/corsika/io/outputs/__init__.py
@@ -9,6 +9,7 @@
from .bethe_bloch import BetheBlochPDG
from .energy_loss import EnergyLoss
+from .interaction import Interactions
from .longitudinal_profile import LongitudinalProfile
from .observation_plane import ObservationPlane
from .output import Output
@@ -28,4 +29,5 @@ __all__ = [
"RadioProcess",
"PrimaryParticle",
"Particle",
+ "Interactions",
]
diff --git a/python/corsika/io/outputs/interaction.py b/python/corsika/io/outputs/interaction.py
new file mode 100644
index 0000000000000000000000000000000000000000..21a1b1d54650066d2d3e664238e8dc3f6286d154
--- /dev/null
+++ b/python/corsika/io/outputs/interaction.py
@@ -0,0 +1,133 @@
+"""
+ Read data written by InteractionWriter.
+
+ (c) Copyright 2024 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.
+"""
+
+import logging
+import os.path as op
+from typing import Any
+
+import pyarrow.parquet as pq
+import yaml
+
+from ..converters import arrow_to_numpy
+from .output import Output
+from .primary import Particle
+
+
+class FirstInteraction(Particle):
+ """
+ Simple class that holds the information of the first interaction
+ This is a complete wrapper of the `Particle` class at the moment
+ but is "new" so that it can be distinguished using `type(thing)`
+ """
+
+ def __init__(self, prop_dict: dict):
+ """
+ Parameters
+ ----------
+ prop_dict: dict
+ the entry from the output yaml file for a single shower
+ """
+ Particle.__init__(self, prop_dict)
+
+ # need to define for momentum prop
+ self.px = self.py = self.pz = None
+
+ def __repr__(self) -> str:
+ """
+ Return a string representation of this class.
+ """
+ out_str = "FirstInteraction:\n"
+ for key in self.__dict__.keys():
+ out_str += f"\t{key}: {self.__dict__[key]}\n"
+ return out_str
+
+ @property
+ def momentum(self) -> list:
+ return [self.px, self.py, self.pz]
+
+
+class Interactions(Output):
+ """
+ Reads the interactions and secondaries
+ """
+
+ def __init__(self, path: str):
+ """
+ Load the particle data into a parquet table.
+
+ Parameters
+ ----------
+ path: str
+ The path to the directory containing this output.
+ """
+ super().__init__(path)
+
+ # try and load our data
+ try:
+ self.__data = pq.read_table(op.join(path, "interactions.parquet"))
+ with open(op.join(path, "summary.yaml"), "r") as f:
+ temp = yaml.load(f, Loader=yaml.Loader)
+ self.__projectiles = [FirstInteraction(x) for x in temp.values()]
+ except Exception as e:
+ logging.getLogger("corsika").warn(
+ f"An error occured loading an Interaction: {e}"
+ )
+
+ @property
+ def projectiles(self) -> list:
+ return self.__projectiles
+
+ def is_good(self) -> bool:
+ """
+ Returns true if this output has been read successfully
+ and has the correct files/state/etc.
+
+ Returns
+ -------
+ bool:
+ True if this is a good output.
+ """
+ return self.__data is not None
+
+ def astype(self, dtype: str = "pandas", **kwargs: Any) -> Any:
+ """
+ Load the particle data of the particle interactions.
+
+ All additional keyword arguments are passed to `parquet.read_table`
+
+ Parameters
+ ----------
+ dtype: str
+ The data format to return the data in (i.e. numpy, pandas, etc.)
+
+ Returns
+ -------
+ Any:
+ The return type of this method is determined by `dtype`.
+ """
+ if dtype == "arrow":
+ return self.__data
+ elif dtype == "pandas":
+ return self.__data.to_pandas()
+ elif dtype == "numpy":
+ return arrow_to_numpy.convert_to_numpy(self.__data)
+ else:
+ raise ValueError(
+ (
+ f"Unknown format '{dtype}' for Interaction. "
+ "We currently only support ['arrow', 'pandas', 'numpy']."
+ )
+ )
+
+ def __repr__(self) -> str:
+ """
+ Return a string representation of this class.
+ """
+ return f"Interaction('{self.config['name']}')"
diff --git a/python/examples/first_interactions.py b/python/examples/first_interactions.py
new file mode 100644
index 0000000000000000000000000000000000000000..d86d8d70dd7cbe21adad95ac6eed4d71b8dfa97a
--- /dev/null
+++ b/python/examples/first_interactions.py
@@ -0,0 +1,114 @@
+#!/usr/bin/python3
+
+import argparse
+import os
+
+import matplotlib.pyplot as plt
+import numpy as np
+import particle
+
+import corsika
+
+here = os.path.abspath(os.path.dirname(__file__))
+
+parser = argparse.ArgumentParser()
+parser.add_argument(
+ "--input-dir", required=True, help="output directory of the CORSIKA 8 simulation"
+)
+parser.add_argument(
+ "--output-dir",
+ default=os.path.join(here, "example_plots"),
+ help="output directory for plots",
+)
+args = parser.parse_args()
+
+if not os.path.isdir(args.output_dir):
+ print("Making directory", args.output_dir)
+ os.makedirs(args.output_dir)
+
+# Load the shower simulation output
+lib = corsika.Library(args.input_dir)
+
+# Load the primary particle information from the shower
+interactions = lib.get("interactions").data
+interaction_config = lib.get("interactions").config
+projectiles = lib.get("interactions").projectiles
+
+shower_ids = np.unique(interactions["shower"])
+
+ncols = 2
+nrows = len(shower_ids)
+fig, ax = plt.subplots(
+ ncols=ncols,
+ nrows=nrows,
+ figsize=(ncols * 8, nrows * 8),
+ gridspec_kw={"wspace": 0.2, "hspace": 0.3},
+)
+ax = np.atleast_2d(ax)
+if ax.shape == (1, 2):
+ ax = ax.T
+
+for ish, sh_id in enumerate(shower_ids):
+
+ # Find all of the secondary particles associated with this shower
+ daughters = interactions[interactions["shower"] == sh_id]
+
+ mother = projectiles[ish]
+ norm = np.sqrt(sum(np.array(mother.momentum) ** 2)) # total momentum
+
+ print("Plotting shower", sh_id)
+ print(mother)
+
+ px, py, pz = mother.momentum
+
+ mother_name = f"${particle.Particle.from_pdgid(mother.pdg).latex_name}$"
+ ax[0, sh_id].plot([-px, 0], [-pz, 0], label=mother_name)
+ ax[1, sh_id].plot([-py, 0], [-pz, 0], label=mother_name)
+
+ # Get daughters, filter out zero-momentum particles
+ p_tot = np.sqrt(daughters["px"] ** 2 + daughters["py"] ** 2 + daughters["pz"] ** 2)
+ daughters = daughters[p_tot > 0]
+
+ # keep track of largest values
+ max_pz = abs(pz)
+ max_pxy = max(abs(px), abs(py))
+
+ # plot the daughter particles
+ for i in range(len(daughters)):
+ pz = daughters["pz"].iloc[i]
+ px = daughters["px"].iloc[i]
+ py = daughters["py"].iloc[i]
+
+ pdg = int(daughters["pdg"].iloc[i])
+ try:
+ name = f"${particle.Particle.from_pdgid(pdg).latex_name}$"
+ except Exception:
+ name = "Unknown"
+
+ ax[0, sh_id].plot([px, 0], [pz, 0], label=name)
+ ax[1, sh_id].plot([py, 0], [pz, 0], label=name)
+
+ max_pz = max(max_pz, abs(pz))
+ max_pxy = max(max_pxy, max(abs(px), abs(py)))
+
+ # set plotting style and names
+ for i in range(2):
+ ax[i, sh_id].legend()
+
+ ymin, ymax = ax[i, sh_id].get_ylim()
+ ax[i, sh_id].set_ylim(-max_pz, max_pz)
+ ax[i, sh_id].set_xlim(-max_pxy, max_pxy)
+
+ unit_energy_str = interaction_config["units"]["energy"]
+ title = f"Mother: {mother.kinetic_energy:0.2e} {unit_energy_str}"
+ title += ", " + mother_name
+ title += f"\n{mother.n_secondaries} secondaries"
+ ax[i, sh_id].set_title(title)
+ ax[i, sh_id].set_ylabel(f"$P_z$ ({unit_energy_str})")
+
+ ax[0, sh_id].set_xlabel(f"$P_x$ ({unit_energy_str})")
+ ax[1, sh_id].set_xlabel(f"$P_y$ ({unit_energy_str})")
+
+plot_path = os.path.join(args.output_dir, "interactions.png")
+print("Saving", plot_path)
+fig.savefig(plot_path, bbox_inches="tight")
diff --git a/tests/framework/testCascade.cpp b/tests/framework/testCascade.cpp
index c6aff27cf29a59676d1b448ceb08c2ff8c92ad19..63da8381b2f93dcdf426580d202f7c3a8482e61d 100644
--- a/tests/framework/testCascade.cpp
+++ b/tests/framework/testCascade.cpp
@@ -184,6 +184,32 @@ private:
int calls_ = 0;
};
+class DummyDecay : public DecayProcess<DummyDecay> {
+ // Dummy decay process that puts pre-predefined particles on stack
+public:
+ DummyDecay(Code code, HEPEnergyType ek, DirectionVector dir)
+ : code_(code)
+ , e0_(ek)
+ , dir_(dir) {}
+
+ Code code_;
+ HEPEnergyType e0_;
+ DirectionVector dir_;
+
+ template <typename Particle>
+ TimeType getLifetime(Particle&) const {
+ return 1_s;
+ }
+
+ template <typename TView>
+ void doDecay(TView& view) const {
+ auto projectile = view.getProjectile();
+ auto const dir = projectile.getMomentum().normalized();
+ projectile.addSecondary(std::make_tuple(code_, e0_, dir_));
+ CORSIKA_LOG_INFO("Particle stack size {}", view.getSize());
+ }
+};
+
TEST_CASE("Cascade", "[Cascade]") {
logging::set_level(logging::level::info);
@@ -196,19 +222,24 @@ TEST_CASE("Cascade", "[Cascade]") {
auto env = make_dummy_env();
auto const& rootCS = env.getCoordinateSystem();
+ // Properties of primary particle
+ auto const primCode = Code::Electron;
+ auto const primEk = E0 - get_mass(Code::Electron);
+ auto const primDir = DirectionVector(rootCS, {0, 0, -1});
+
StackInspector<TestCascadeStack> stackInspect(100, true, E0);
NullModel nullModel;
+ DummyDecay decay(primCode, primEk, primDir); // decay product will be primary
HEPEnergyType const Ecrit = 85_MeV;
ProcessSplit split;
ProcessCut cut(Ecrit);
- auto sequence = make_sequence(nullModel, stackInspect, split, cut);
+ auto sequence = make_sequence(nullModel, decay, stackInspect, split, cut);
TestCascadeStack stack;
+
stack.clear();
- stack.addParticle(std::make_tuple(Code::Electron,
- E0 - get_mass(Code::Electron), // Ekin
- DirectionVector(rootCS, {0, 0, -1}),
- Point(rootCS, {0_m, 0_m, 10_km}), 0_ns));
+ stack.addParticle(
+ std::make_tuple(primCode, primEk, primDir, Point(rootCS, {0_m, 0_m, 10_km}), 0_ns));
DummyTracking tracking;
DummyOutputManager output;
@@ -235,6 +266,36 @@ TEST_CASE("Cascade", "[Cascade]") {
CHECK(stack.getSize() == 13);
CHECK(split.getCalls() == 2047);
}
+
+ SECTION("double_forcing_1") {
+ EAS.forceInteraction();
+ REQUIRE_THROWS(EAS.forceDecay());
+ }
+
+ SECTION("double_forcing_2") {
+ EAS.forceDecay();
+ REQUIRE_THROWS(EAS.forceInteraction());
+ }
+
+ SECTION("forced decay") {
+ // Put anything new so that it won't trigger "decays into self" error
+ stack.clear();
+ stack.addParticle(std::make_tuple(Code::PiPlus,
+ 0.5_GeV, // Ekin
+ DirectionVector(rootCS, {0, 0, -1}),
+ Point(rootCS, {0_m, 0_m, 10_km}), 0_ns));
+
+ CHECK(tracking.getCalls() == 0);
+ EAS.forceDecay();
+ CHECK(tracking.getCalls() == 0);
+ CHECK(stack.getEntries() == 1);
+ CHECK(stack.getSize() == 1);
+ EAS.run();
+ CHECK(tracking.getCalls() == 2047);
+ CHECK(stack.getEntries() == 0);
+ CHECK(stack.getSize() == 14);
+ CHECK(split.getCalls() == 2047);
+ }
}
TEST_CASE("Cascade Zero Interaction", "[Cascade]") {
@@ -273,4 +334,4 @@ TEST_CASE("Cascade Zero Interaction", "[Cascade]") {
// expect 100 calls. if it would be less, this means that the particle has been erased
// early by the Cascade.inl algorithm
CHECK(counter.getCalls() == 100);
-}
\ No newline at end of file
+}
diff --git a/tests/output/CMakeLists.txt b/tests/output/CMakeLists.txt
index f85f95aa71a5aeb00e768eb5aa85286ec3f67466..3b954de5f94129ad897fd8785a8aa676b6d75e01 100644
--- a/tests/output/CMakeLists.txt
+++ b/tests/output/CMakeLists.txt
@@ -1,7 +1,8 @@
set (test_output_sources
+ testDummyOutputManager.cpp
+ testInteractionWriter.cpp
TestMain.cpp
testOutputManager.cpp
- testDummyOutputManager.cpp
testParquetStreamer.cpp
testWriterObservationPlane.cpp
testWriterTracks.cpp
diff --git a/tests/output/testInteractionWriter.cpp b/tests/output/testInteractionWriter.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..d9a45e41883643dbf8a841a609dd137f43be022d
--- /dev/null
+++ b/tests/output/testInteractionWriter.cpp
@@ -0,0 +1,124 @@
+/*
+ * (c) Copyright 2024 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.
+ */
+
+#include <boost/filesystem.hpp>
+
+#include <catch2/catch.hpp>
+
+#include <corsika/modules/writers/InteractionWriter.hpp>
+
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/geometry/RootCoordinateSystem.hpp>
+#include <corsika/framework/geometry/Vector.hpp>
+
+#include <corsika/media/Environment.hpp>
+#include <corsika/media/HomogeneousMedium.hpp>
+
+#include <corsika/setup/SetupTrajectory.hpp>
+
+#include <SetupTestStack.hpp>
+#include <SetupTestTrajectory.hpp>
+
+using namespace corsika;
+
+using DummyEnvironmentInterface = IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>;
+using DummyEnvironment = Environment<DummyEnvironmentInterface>;
+
+TEST_CASE("InteractionWriter", "process") {
+ logging::set_level(logging::level::info);
+
+ auto env = std::make_unique<Environment<IMediumModel>>();
+ auto& universe = *(env->getUniverse());
+ CoordinateSystemPtr const& rootCS = env->getCoordinateSystem();
+ auto theMedium = Environment<IMediumModel>::createNode<Sphere>(
+ Point{rootCS, 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity());
+ using MyHomogeneousModel = HomogeneousMedium<IMediumModel>;
+ const auto density = 1_kg / (1_m * 1_m * 1_m);
+ theMedium->setModelProperties<MyHomogeneousModel>(
+ density, NuclearComposition({Code::Nitrogen}, {1.}));
+ universe.addChild(std::move(theMedium));
+
+ std::string const outputDir = "./output_interactions";
+
+ // preparation
+ if (boost::filesystem::exists(outputDir)) { boost::filesystem::remove_all(outputDir); }
+ boost::filesystem::create_directory(outputDir);
+
+ // setup empty particle stack
+ test::Stack stack;
+ stack.clear();
+
+ // Make shower axis
+ auto start = Point(rootCS, -1_m, 0_m, 0_m);
+ auto stop = Point(rootCS, 1_m, 0_m, 0_m);
+ auto length = stop - start;
+ ShowerAxis axis(start, length, *env);
+
+ Plane const plane(Point(rootCS, {0_m, 0_m, 0_m}),
+ DirectionVector(rootCS, {0., 0., 1.}));
+
+ ObservationPlane<setup::Tracking, ParticleWriterParquet> obsPlane{
+ plane, DirectionVector(rootCS, {1., 0., 0.}),
+ true, // plane should "absorb" particles
+ false};
+
+ InteractionWriter<setup::Tracking, ParticleWriterParquet> writer(axis, obsPlane);
+ writer.startOfLibrary(outputDir);
+ writer.startOfShower(0);
+
+ CHECK(!writer.getInteractionCounter());
+
+ // add primary particle to stack
+ auto particle = stack.addParticle(std::make_tuple(Code::Proton, 10_GeV,
+ DirectionVector(rootCS, {1, 0, 0}),
+ Point(rootCS, 0_m, 0_m, 0_m), 0_ns));
+
+ test::StackView view(particle);
+ auto projectile = view.getProjectile();
+
+ std::vector<Code> const particleList = {Code::PiPlus, Code::PiMinus, Code::KPlus,
+ Code::KMinus, Code::K0Long, Code::K0Short,
+ Code::Electron, Code::MuPlus, Code::NuE,
+ Code::Neutron, Code::NuMu};
+
+ for (auto proType : particleList)
+ projectile.addSecondary(
+ std::make_tuple(proType, 1_GeV, DirectionVector(rootCS, {1, 0, 0})));
+
+ // pre-check on stack size
+ CHECK(view.getEntries() == 11);
+ CHECK(stack.getEntries() == 12);
+
+ // correctly initialized counter
+ CHECK(!writer.getInteractionCounter());
+
+ // run the process and check for consistency
+ writer.doSecondaries(view);
+ CHECK(1 == writer.getInteractionCounter());
+ CHECK(view.getEntries() == 11);
+ CHECK(stack.getEntries() == 12);
+
+ // ensure the counter bumps
+ writer.doSecondaries(view);
+ CHECK(2 == writer.getInteractionCounter());
+
+ writer.endOfShower(0);
+ writer.endOfLibrary();
+
+ CHECK(boost::filesystem::exists(outputDir + "/interactions.parquet"));
+
+ auto const config = writer.getConfig();
+ CHECK(config["type"].as<std::string>() == "Interactions");
+
+ auto const summary = writer.getSummary();
+ CHECK(summary["shower_0"]["n_secondaries"].as<int>() == 11);
+
+ // clean up
+ if (boost::filesystem::exists(outputDir)) { boost::filesystem::remove_all(outputDir); }
+}