From 1c170b3b3f520d1e6c5a58cf010419184eb76139 Mon Sep 17 00:00:00 2001
From: Alan Coleman <alanc@udel.edu>
Date: Fri, 8 Mar 2024 10:38:33 +0100
Subject: [PATCH] first draft of interaction writer
---
.../modules/writers/InteractionWriter.inl | 119 ++++++++++++++++++
corsika/modules/writers/InteractionWriter.hpp | 61 +++++++++
tests/output/CMakeLists.txt | 3 +-
tests/output/testInteractionWriter.cpp | 119 ++++++++++++++++++
4 files changed, 301 insertions(+), 1 deletion(-)
create mode 100644 corsika/detail/modules/writers/InteractionWriter.inl
create mode 100644 corsika/modules/writers/InteractionWriter.hpp
create mode 100644 tests/output/testInteractionWriter.cpp
diff --git a/corsika/detail/modules/writers/InteractionWriter.inl b/corsika/detail/modules/writers/InteractionWriter.inl
new file mode 100644
index 000000000..61784ae1c
--- /dev/null
+++ b/corsika/detail/modules/writers/InteractionWriter.inl
@@ -0,0 +1,119 @@
+/*
+ * (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 {
+
+ inline InteractionWriter::InteractionWriter(ShowerAxis const& axis)
+ : showerAxis_(axis)
+ , interactionCounter_(0)
+ , showerId_(0)
+ , nSecondaries_(0)
+ {}
+
+ template <typename TStackView>
+ inline void InteractionWriter::doSecondaries(TStackView& vS) {
+ if (interactionCounter_++) {
+ return;
+ }
+
+ 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());
+
+ *(output_.getWriter()) << showerId_ << true
+ << static_cast<int>(primary.getPID())
+ << static_cast<float>(primary.getKineticEnergy() / 1_GeV)
+ << static_cast<float>(1.0)
+ << static_cast<float>(0.0)
+ << static_cast<float>(0.0)
+ << static_cast<double>(primary.getTime() / 1_ns)
+ << static_cast<float>(dX / (1_g / 1_cm / 1_cm))
+ << parquet::EndRow;
+
+ // Loop through secondaries
+ auto particle = vS.begin();
+ while (particle != vS.end()) {
+
+ *(output_.getWriter()) << showerId_ << false
+ << static_cast<int>(particle.getPID())
+ << static_cast<float>(particle.getKineticEnergy() / 1_GeV)
+ << static_cast<float>(1.0)
+ << static_cast<float>(0.0)
+ << static_cast<float>(0.0)
+ << 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(), particle.getKineticEnergy());
+ ++particle;
+ ++nSecondaries_;
+ }
+ }
+
+ inline void InteractionWriter::startOfLibrary(boost::filesystem::path const& directory) {
+ output_.initStreamer((directory / ("interactions.parquet")).string());
+
+ output_.addField("primary", parquet::Repetition::REQUIRED, parquet::Type::BOOLEAN,
+ parquet::ConvertedType::NONE);
+ output_.addField("pdg", parquet::Repetition::REQUIRED, parquet::Type::INT32,
+ parquet::ConvertedType::INT_32);
+ output_.addField("kinetic_energy", parquet::Repetition::REQUIRED,
+ parquet::Type::FLOAT, parquet::ConvertedType::NONE);
+ output_.addField("nx", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
+ parquet::ConvertedType::NONE);
+ output_.addField("ny", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
+ parquet::ConvertedType::NONE);
+ output_.addField("nz", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
+ parquet::ConvertedType::NONE);
+ output_.addField("time", parquet::Repetition::REQUIRED, parquet::Type::DOUBLE,
+ parquet::ConvertedType::NONE);
+ output_.addField("dX", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
+ parquet::ConvertedType::NONE);
+
+ output_.buildStreamer();
+
+ showerId_ = 0;
+ interactionCounter_ = 0;
+ nSecondaries_ = 0;
+ summary_ = YAML::Node();
+ }
+
+ inline void InteractionWriter::startOfShower(unsigned int const showerId) {
+ showerId_ = showerId;
+ interactionCounter_ = 0;
+ nSecondaries_ = 0;
+ }
+
+ inline void InteractionWriter::endOfShower(unsigned int const) {
+ summary_["shower_" + std::to_string(showerId_)]["n_secondaries"] = nSecondaries_;
+ }
+
+ inline void InteractionWriter::endOfLibrary() {
+ output_.closeStreamer();
+ }
+
+ inline YAML::Node InteractionWriter::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;
+ }
+
+ inline YAML::Node InteractionWriter::getSummary() const {
+ return summary_;
+ }
+
+} // namespace corsika
diff --git a/corsika/modules/writers/InteractionWriter.hpp b/corsika/modules/writers/InteractionWriter.hpp
new file mode 100644
index 000000000..6ba3fee30
--- /dev/null
+++ b/corsika/modules/writers/InteractionWriter.hpp
@@ -0,0 +1,61 @@
+/*
+ * (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/output/BaseOutput.hpp>
+#include <corsika/output/ParquetStreamer.hpp>
+
+namespace corsika {
+
+ class InteractionWriter : public SecondariesProcess<InteractionWriter>, BaseOutput{
+
+ public:
+ /**
+ *
+ * @param axis - axis used for calculating grammage
+ */
+ InteractionWriter(ShowerAxis const& axis);
+
+ /**
+ *
+ * @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:
+ ShowerAxis const& showerAxis_;
+ unsigned int interactionCounter_;
+ unsigned int showerId_;
+ unsigned int nSecondaries_;
+
+ ParquetStreamer output_;
+ YAML::Node summary_;
+
+ }; // namespace corsika
+
+} // namespace corsika
+
+#include <corsika/detail/modules/writers/InteractionWriter.inl>
diff --git a/tests/output/CMakeLists.txt b/tests/output/CMakeLists.txt
index f85f95aa7..3b954de5f 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 000000000..302c0ccd6
--- /dev/null
+++ b/tests/output/testInteractionWriter.cpp
@@ -0,0 +1,119 @@
+/*
+ * (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 <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);
+
+
+ InteractionWriter writer(axis);
+ 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); }
+}
--
GitLab