diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index a52219a52896ab6382d6ab72148da5288d041124..53f2e7b4ce54b88a9e45df329eb9a7348829330c 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -33,6 +33,7 @@ stages:
- build_test
- example
- build_test_example
+ - python
- install
- optional
@@ -542,3 +543,38 @@ sanity:
untracked: true
policy: pull
key: "${CI_COMMIT_REF_SLUG}-gcc"
+
+##########################################################
+# template for all Python jobs
+.python:
+ stage: python
+ tags:
+ - corsika
+ before_script:
+ - apt-get update && apt-get install -y python3-pip
+ script:
+ - cd ${CI_PROJECT_DIR}/python # change into the Python directory
+ - pip3 install --user -e '.[test]' # install the package + test deps
+ - python3 --version
+ - python3 -m mypy corsika
+ - python3 -m isort --atomic --check-only corsika tests
+ - python3 -m black -t py37 corsika tests
+ - python3 -m flake8 corsika tests
+ - python3 -m pytest --cov=corsika tests
+ - cd ${CI_PROJECT_DIR} # reset the directory
+ coverage: '/^TOTAL\s*\d+\s*\d+\s*(.*\%)/'
+
+# the default Python version Ubuntu 18.04 is Python3.8
+python-3.8:
+ extends: .python
+ image: corsika/devel:u-18.04
+ dependencies:
+ - build-u-18_04
+ cache:
+ key: "${CI_COMMIT_REF_SLUG}-gcc"
+ artifacts:
+ when: always
+ expire_in: 1 year
+ paths:
+ - ${CI_PROJECT_DIR}/Python/python-test.log
+ allow_failure: true
diff --git a/.gitmodules b/.gitmodules
index 1bbe5d4d093457574b332f922f45b43b2d2f8b53..5429f2d75f591d52f001615fddde188e8bd012de 100644
--- a/.gitmodules
+++ b/.gitmodules
@@ -12,4 +12,4 @@
path = modules/conex
url = ../../AirShowerPhysics/cxroot.git
branch = master
- shallow = true
+ shallow = true
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 74c59086064789144e328c4d1bc911c410658d5e..d81de38f705fa7f6b5901a5643e74cb2f5f9ab93 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -105,6 +105,7 @@ add_compile_options (
include (conan) # self-provided in 'cmake' directory
#
# download and build all dependencies
+message (STATUS "Installing dependencies from Conan (this may take some time)...")
conan_cmake_run (CONANFILE conanfile.txt
BASIC_SETUP CMAKE_TARGETS
BUILD missing
@@ -210,6 +211,8 @@ target_link_libraries (
CONAN_PKG::eigen
CONAN_PKG::spdlog
CONAN_PKG::boost
+ CONAN_PKG::yaml-cpp
+ CONAN_PKG::arrow
cnpy # for SaveBoostHistogram
)
diff --git a/MAINTAINERS.md b/MAINTAINERS.md
index a57e4bfbbbff4526527819a3f2e4840c58a9b799..2da518cf5b6b9a8e0e3a278bf166fd246137554d 100644
--- a/MAINTAINERS.md
+++ b/MAINTAINERS.md
@@ -23,6 +23,13 @@ Hadron models:
- Felix Riehn <friehn@lip.pt>, Santiago/Lisbon
- Anatoli Fedynitch <anatoli.fedynitch@icecube.wisc.edu> ICRR Tokyo
+Output formats and infrastructure:
+- Remy Prechelt <prechelt@hawaii.edu>, UHM
+- Ralf Ulrich <ralf.ulrich@kit.edu>, KIT
+
+Python library:
+- Remy Prechelt <prechelt@hawaii.edu>, UHM
+
Radio:
- Remy Prechelt <prechelt@hawaii.edu>
- Tim Huege <tim.huege@kit.edu>, KIT
diff --git a/corsika/detail/framework/core/Cascade.inl b/corsika/detail/framework/core/Cascade.inl
index 32127908c714d79dfdd7e22e74615aed449b1b7c..654dfc6d1aad0b603baa76e701f5357fa85ffc13 100644
--- a/corsika/detail/framework/core/Cascade.inl
+++ b/corsika/detail/framework/core/Cascade.inl
@@ -26,11 +26,14 @@
namespace corsika {
- template <typename TTracking, typename TProcessList, typename TStack,
+ template <typename TTracking, typename TProcessList, typename TOutput, typename TStack,
typename TStackView>
- inline void Cascade<TTracking, TProcessList, TStack, TStackView>::run() {
+ inline void Cascade<TTracking, TProcessList, TOutput, TStack, TStackView>::run() {
setNodes(); // put each particle on stack in correct environment volume
+ // start this event (i.e. this shower)
+ output_.startOfShower();
+
while (!stack_.isEmpty()) {
while (!stack_.isEmpty()) {
CORSIKA_LOG_TRACE("Stack: {}", stack_.asString());
@@ -51,11 +54,15 @@ namespace corsika {
// thus, the double loop
// doCascadeEquations();
}
+
+ // end this event (i.e. this shower)
+ output_.endOfShower();
}
- template <typename TTracking, typename TProcessList, typename TStack,
+ template <typename TTracking, typename TProcessList, typename TOutput, typename TStack,
typename TStackView>
- inline void Cascade<TTracking, TProcessList, TStack, TStackView>::forceInteraction() {
+ inline void
+ Cascade<TTracking, TProcessList, TOutput, TStack, TStackView>::forceInteraction() {
CORSIKA_LOG_TRACE("forced interaction!");
setNodes();
auto vParticle = stack_.getNextParticle();
@@ -65,9 +72,9 @@ namespace corsika {
vParticle.erase(); // primary particle is done
}
- template <typename TTracking, typename TProcessList, typename TStack,
+ template <typename TTracking, typename TProcessList, typename TOutput, typename TStack,
typename TStackView>
- inline void Cascade<TTracking, TProcessList, TStack, TStackView>::step(
+ inline void Cascade<TTracking, TProcessList, TOutput, TStack, TStackView>::step(
Particle& vParticle) {
// determine combined total interaction length (inverse)
@@ -258,9 +265,10 @@ namespace corsika {
vParticle.erase();
}
- template <typename TTracking, typename TProcessList, typename TStack,
+ template <typename TTracking, typename TProcessList, typename TOutput, typename TStack,
typename TStackView>
- inline ProcessReturn Cascade<TTracking, TProcessList, TStack, TStackView>::decay(
+ inline ProcessReturn
+ Cascade<TTracking, TProcessList, TOutput, TStack, TStackView>::decay(
TStackView& view, InverseTimeType initial_inv_decay_time) {
CORSIKA_LOG_DEBUG("decay");
@@ -289,11 +297,11 @@ namespace corsika {
return returnCode;
}
- template <typename TTracking, typename TProcessList, typename TStack,
+ template <typename TTracking, typename TProcessList, typename TOutput, typename TStack,
typename TStackView>
- inline ProcessReturn Cascade<TTracking, TProcessList, TStack, TStackView>::interaction(
+ inline ProcessReturn
+ Cascade<TTracking, TProcessList, TOutput, TStack, TStackView>::interaction(
TStackView& view, InverseGrammageType initial_inv_int_length) {
-
CORSIKA_LOG_DEBUG("collide");
#ifdef DEBUG
@@ -322,9 +330,9 @@ namespace corsika {
return returnCode;
}
- template <typename TTracking, typename TProcessList, typename TStack,
+ template <typename TTracking, typename TProcessList, typename TOutput, typename TStack,
typename TStackView>
- inline void Cascade<TTracking, TProcessList, TStack, TStackView>::setNodes() {
+ inline void Cascade<TTracking, TProcessList, TOutput, TStack, TStackView>::setNodes() {
std::for_each(stack_.begin(), stack_.end(), [&](auto& p) {
auto const* numericalNode =
environment_.getUniverse()->getContainingNode(p.getPosition());
@@ -332,9 +340,9 @@ namespace corsika {
});
}
- template <typename TTracking, typename TProcessList, typename TStack,
+ template <typename TTracking, typename TProcessList, typename TOutput, typename TStack,
typename TStackView>
- inline void Cascade<TTracking, TProcessList, TStack, TStackView>::setEventType(
+ inline void Cascade<TTracking, TProcessList, TOutput, TStack, TStackView>::setEventType(
TStackView& view, [[maybe_unused]] history::EventType eventType) {
if constexpr (TStackView::has_event) {
for (auto&& sec : view) { sec.getEvent()->setEventType(eventType); }
diff --git a/corsika/detail/framework/utility/CubicSolver.inl b/corsika/detail/framework/utility/CubicSolver.inl
index 3a68a3d6ebcf39290807c3402c09fcfccfe8ab6e..c15594e02a02bef2683e4d3ae4d4ad02895ef747 100644
--- a/corsika/detail/framework/utility/CubicSolver.inl
+++ b/corsika/detail/framework/utility/CubicSolver.inl
@@ -248,10 +248,6 @@ namespace corsika {
long double x1 = xinfl;
long double f_x1 = cubic_function(xinfl, a, b, c, d);
- CORSIKA_LOG_TRACE("dist={} xinfl={} f_x1={} {} {}", dist, xinfl, f_x1,
- cubic_function(xinfl - 2 / 3 * std::sqrt(dist), a, b, c, d),
- cubic_function(xinfl + 2 / 3 * std::sqrt(dist), a, b, c, d));
-
if (std::abs(f_x1) > epsilon) {
if (std::abs(dist) < epsilon) {
x1 = xinfl - std::cbrt(f_x1);
diff --git a/corsika/detail/modules/ObservationPlane.inl b/corsika/detail/modules/ObservationPlane.inl
index e5be30bedf2f81b65d68faeca314669c58c31ac6..12614573f13888396cc28e3655b9119fa7af2bec 100644
--- a/corsika/detail/modules/ObservationPlane.inl
+++ b/corsika/detail/modules/ObservationPlane.inl
@@ -6,33 +6,26 @@
* the license.
*/
-#include <corsika/modules/ObservationPlane.hpp>
-#include <corsika/framework/core/Logging.hpp>
+#include <corsika/setup/SetupStack.hpp>
#include <corsika/setup/SetupTrajectory.hpp>
-#include <fstream>
-
namespace corsika {
- inline ObservationPlane::ObservationPlane(Plane const& obsPlane,
- DirectionVector const& x_axis,
- std::string const& filename, bool deleteOnHit)
+ template <typename TOutput>
+ ObservationPlane<TOutput>::ObservationPlane(Plane const& obsPlane,
+ DirectionVector const& x_axis,
+ bool deleteOnHit)
: plane_(obsPlane)
- , outputStream_(filename)
, deleteOnHit_(deleteOnHit)
, energy_ground_(0_GeV)
, count_ground_(0)
, xAxis_(x_axis.normalized())
- , yAxis_(obsPlane.getNormal().cross(xAxis_)) {
- CORSIKA_LOG_DEBUG("Plane height: {}", obsPlane.getCenter());
- outputStream_ << "#PDG code, energy / eV, x distance / m, y distance / m"
- << std::endl;
- }
+ , yAxis_(obsPlane.getNormal().cross(xAxis_)) {}
- inline ProcessReturn ObservationPlane::doContinuous(
+ template <typename TOutput>
+ inline ProcessReturn ObservationPlane<TOutput>::doContinuous(
corsika::setup::Stack::particle_type& particle, corsika::setup::Trajectory&,
bool const stepLimit) {
-
/*
The current step did not yet reach the ObservationPlane, do nothing now and wait:
*/
@@ -53,9 +46,9 @@ namespace corsika {
Point const pointOfIntersection = particle.getPosition();
Vector const displacement = pointOfIntersection - plane_.getCenter();
- outputStream_ << static_cast<int>(get_PDG(particle.getPID())) << ' ' << energy / 1_eV
- << ' ' << displacement.dot(xAxis_) / 1_m << ' '
- << displacement.dot(yAxis_) / 1_m << '\n';
+ // add our particles to the output file stream
+ this->write(particle.getPID(), energy, displacement.dot(xAxis_),
+ displacement.dot(yAxis_));
CORSIKA_LOG_TRACE("Particle detected absorbed={}", deleteOnHit_);
@@ -68,7 +61,8 @@ namespace corsika {
}
}
- inline LengthType ObservationPlane::getMaxStepLength(
+ template <typename TOutput>
+ inline LengthType ObservationPlane<TOutput>::getMaxStepLength(
corsika::setup::Stack::particle_type const& particle,
corsika::setup::Trajectory const& trajectory) {
@@ -93,17 +87,62 @@ namespace corsika {
return dist;
}
- inline void ObservationPlane::showResults() const {
+ template <typename TOutput>
+ inline void ObservationPlane<TOutput>::showResults() const {
CORSIKA_LOG_INFO(
" ******************************\n"
" ObservationPlane: \n"
- " energy in ground (GeV) : {}\n"
- " no. of particles in ground : {}\n"
+ " energy an ground (GeV) : {}\n"
+ " no. of particles at ground : {}\n"
" ******************************",
energy_ground_ / 1_GeV, count_ground_);
}
- inline void ObservationPlane::reset() {
+ template <typename TOutput>
+ inline YAML::Node ObservationPlane<TOutput>::getConfig() const {
+ using namespace units::si;
+
+ // construct the top-level node
+ YAML::Node node;
+
+ // basic info
+ node["type"] = "ObservationPlane";
+ node["units"] = "m"; // add default units for values
+
+ // the center of the plane
+ auto const center{plane_.getCenter()};
+
+ // save each component in its native coordinate system
+ auto const center_coords{center.getCoordinates(center.getCoordinateSystem())};
+ node["plane"]["center"].push_back(center_coords.getX() / 1_m);
+ node["plane"]["center"].push_back(center_coords.getY() / 1_m);
+ node["plane"]["center"].push_back(center_coords.getZ() / 1_m);
+
+ // the normal vector of the plane
+ auto const normal{plane_.getNormal().getComponents()};
+ node["plane"]["normal"].push_back(normal.getX().magnitude());
+ node["plane"]["normal"].push_back(normal.getY().magnitude());
+ node["plane"]["normal"].push_back(normal.getZ().magnitude());
+
+ // the x-axis vector
+ auto const xAxis_coords{xAxis_.getComponents(xAxis_.getCoordinateSystem())};
+ node["x-axis"].push_back(xAxis_coords.getX().magnitude());
+ node["x-axis"].push_back(xAxis_coords.getY().magnitude());
+ node["x-axis"].push_back(xAxis_coords.getZ().magnitude());
+
+ // the y-axis vector
+ auto const yAxis_coords{yAxis_.getComponents(yAxis_.getCoordinateSystem())};
+ node["y-axis"].push_back(yAxis_coords.getX().magnitude());
+ node["y-axis"].push_back(yAxis_coords.getY().magnitude());
+ node["y-axis"].push_back(yAxis_coords.getZ().magnitude());
+
+ node["delete_on_hit"] = deleteOnHit_;
+
+ return node;
+ }
+
+ template <typename TOutput>
+ inline void ObservationPlane<TOutput>::reset() {
energy_ground_ = 0_GeV;
count_ground_ = 0;
}
diff --git a/corsika/detail/modules/TrackWriter.inl b/corsika/detail/modules/TrackWriter.inl
index e55d40cb9cfb64d23aaccd21fc11347e10035226..8af3c31dc45b718d8213346e3896e2dffb02812d 100644
--- a/corsika/detail/modules/TrackWriter.inl
+++ b/corsika/detail/modules/TrackWriter.inl
@@ -8,59 +8,48 @@
#pragma once
-#include <corsika/modules/TrackWriter.hpp>
-
#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
-#include <corsika/setup/SetupStack.hpp>
-#include <corsika/setup/SetupTrajectory.hpp>
-
-#include <iomanip>
#include <limits>
namespace corsika {
- inline TrackWriter::TrackWriter(std::string const& filename)
- : filename_(filename) {
- using namespace std::string_literals;
-
- file_.open(filename_);
- file_
- << "# PID, E / eV, start coordinates / m, displacement vector to end / m, steplength / m "s
- << '\n';
- }
-
- inline TrackWriter::~TrackWriter() { file_.close(); }
+ template <typename TOutputWriter>
+ inline TrackWriter<TOutputWriter>::TrackWriter() {}
+ template <typename TOutputWriter>
template <typename TParticle, typename TTrack>
- inline ProcessReturn TrackWriter::doContinuous(TParticle const& vP, TTrack const& vT,
- bool const) {
-
- CORSIKA_LOG_DEBUG("TrackWriter");
+ inline ProcessReturn TrackWriter<TOutputWriter>::doContinuous(TParticle const& vP,
+ TTrack const& vT,
+ bool const) {
auto const start = vT.getPosition(0).getCoordinates();
- auto const delta = vT.getPosition(1).getCoordinates() - start;
- auto const pdg = static_cast<int>(get_PDG(vP.getPID()));
+ auto const end = vT.getPosition(1).getCoordinates();
- // clang-format off
- file_ << std::setw(7) << pdg
- << std::setw(width_) << std::scientific << std::setprecision(precision_) << vP.getEnergy() / 1_eV
- << std::setw(width_) << std::scientific << std::setprecision(precision_) << start[0] / 1_m
- << std::setw(width_) << std::scientific << std::setprecision(precision_) << start[1] / 1_m
- << std::setw(width_) << std::scientific << std::setprecision(precision_) << start[2] / 1_m
- << std::setw(width_) << std::scientific << std::setprecision(precision_) << delta[0] / 1_m
- << std::setw(width_) << std::scientific << std::setprecision(precision_) << delta[1] / 1_m
- << std::setw(width_) << std::scientific << std::setprecision(precision_) << delta[2] / 1_m
- << std::setw(width_) << std::scientific << std::setprecision(precision_) << delta.getNorm() / 1_m
- << '\n';
- // clang-format on
+ // write the track to the file
+ this->write(vP.getPID(), vP.getEnergy(), start, end);
return ProcessReturn::Ok;
}
+ template <typename TOutputWriter>
template <typename TParticle, typename TTrack>
- inline LengthType TrackWriter::getMaxStepLength(TParticle const&, TTrack const&) {
+ inline LengthType TrackWriter<TOutputWriter>::getMaxStepLength(TParticle const&,
+ TTrack const&) {
return meter * std::numeric_limits<double>::infinity();
}
+ template <typename TOutputWriter>
+ YAML::Node TrackWriter<TOutputWriter>::getConfig() const {
+ using namespace units::si;
+
+ YAML::Node node;
+
+ // add default units for values
+ node["type"] = "TrackWriter";
+ node["units"] = "GeV | m";
+
+ return node;
+ }
} // namespace corsika
diff --git a/corsika/detail/modules/writers/ObservationPlaneWriterParquet.inl b/corsika/detail/modules/writers/ObservationPlaneWriterParquet.inl
new file mode 100644
index 0000000000000000000000000000000000000000..60752c56a9fed06b326449530c2d24d5837d233e
--- /dev/null
+++ b/corsika/detail/modules/writers/ObservationPlaneWriterParquet.inl
@@ -0,0 +1,54 @@
+/*
+ * (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
+
+namespace corsika {
+
+ inline ObservationPlaneWriterParquet::ObservationPlaneWriterParquet()
+ : output_() {}
+
+ inline void ObservationPlaneWriterParquet::startOfLibrary(
+ boost::filesystem::path const& directory) {
+
+ // setup the streamer
+ output_.initStreamer((directory / "particles.parquet").string());
+
+ // enable compression with the default level
+ output_.enableCompression();
+
+ // build the schema
+ output_.addField("pdg", parquet::Repetition::REQUIRED, parquet::Type::INT32,
+ parquet::ConvertedType::INT_32);
+ output_.addField("energy", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
+ parquet::ConvertedType::NONE);
+ output_.addField("x", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
+ parquet::ConvertedType::NONE);
+ output_.addField("y", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
+ parquet::ConvertedType::NONE);
+
+ // and build the streamer
+ output_.buildStreamer();
+ }
+
+ inline void ObservationPlaneWriterParquet::endOfShower() { ++shower_; }
+
+ inline void ObservationPlaneWriterParquet::endOfLibrary() { output_.closeStreamer(); }
+
+ inline void ObservationPlaneWriterParquet::write(Code const& pid,
+ HEPEnergyType const& energy,
+ LengthType const& x,
+ LengthType const& y) {
+ // write the next row - we must write `shower_` first.
+ *(output_.getWriter()) << shower_ << static_cast<int>(get_PDG(pid))
+ << static_cast<float>(energy / 1_GeV)
+ << static_cast<float>(x / 1_m) << static_cast<float>(y / 1_m)
+ << parquet::EndRow;
+ }
+
+} // namespace corsika
diff --git a/corsika/detail/modules/writers/TrackWriterParquet.inl b/corsika/detail/modules/writers/TrackWriterParquet.inl
new file mode 100644
index 0000000000000000000000000000000000000000..26c62bdf0e5e8260c431dc114aaed25e4aab0a82
--- /dev/null
+++ b/corsika/detail/modules/writers/TrackWriterParquet.inl
@@ -0,0 +1,68 @@
+/*
+ * (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
+
+namespace corsika {
+
+ inline TrackWriterParquet::TrackWriterParquet()
+ : output_() {}
+
+ inline void TrackWriterParquet::startOfLibrary(
+ boost::filesystem::path const& directory) {
+
+ // setup the streamer
+ output_.initStreamer((directory / "tracks.parquet").string());
+
+ // build the schema
+ output_.addField("pdg", parquet::Repetition::REQUIRED, parquet::Type::INT32,
+ parquet::ConvertedType::INT_32);
+ output_.addField("energy", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
+ parquet::ConvertedType::NONE);
+ output_.addField("start_x", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
+ parquet::ConvertedType::NONE);
+ output_.addField("start_y", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
+ parquet::ConvertedType::NONE);
+ output_.addField("start_z", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
+ parquet::ConvertedType::NONE);
+ output_.addField("end_x", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
+ parquet::ConvertedType::NONE);
+ output_.addField("end_y", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
+ parquet::ConvertedType::NONE);
+ output_.addField("end_z", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
+ parquet::ConvertedType::NONE);
+
+ // and build the streamer
+ output_.buildStreamer();
+ }
+
+ inline void TrackWriterParquet::endOfShower() { ++shower_; }
+
+ inline void TrackWriterParquet::endOfLibrary() { output_.closeStreamer(); }
+
+ inline void TrackWriterParquet::write(Code const& pid, HEPEnergyType const& energy,
+ QuantityVector<length_d> const& start,
+ QuantityVector<length_d> const& end) {
+
+ // write the next row - we must write `shower_` first.
+ // clang-format off
+ *(output_.getWriter())
+ << shower_
+ << static_cast<int>(get_PDG(pid))
+ << static_cast<float>(energy / 1_GeV)
+ << static_cast<float>(start[0] / 1_m)
+ << static_cast<float>(start[1] / 1_m)
+ << static_cast<float>(start[2] / 1_m)
+ << static_cast<float>(end[0] / 1_m)
+ << static_cast<float>(end[1] / 1_m)
+ << static_cast<float>(end[2] / 1_m)
+ << parquet::EndRow;
+ // clang-format on
+ }
+
+} // namespace corsika
diff --git a/corsika/detail/output/BaseOutput.inl b/corsika/detail/output/BaseOutput.inl
new file mode 100644
index 0000000000000000000000000000000000000000..24ec861ac9c2b3083ea834261f7fa0a57fd8559e
--- /dev/null
+++ b/corsika/detail/output/BaseOutput.inl
@@ -0,0 +1,15 @@
+/*
+ * (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
+
+namespace corsika {
+
+ inline BaseOutput::BaseOutput()
+ : shower_(0) {}
+
+} // namespace corsika
diff --git a/corsika/detail/output/DummyOutputManager.inl b/corsika/detail/output/DummyOutputManager.inl
new file mode 100644
index 0000000000000000000000000000000000000000..acf2b69f987dbebeb69c9a77438141afc36a036d
--- /dev/null
+++ b/corsika/detail/output/DummyOutputManager.inl
@@ -0,0 +1,27 @@
+/*
+ * (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
+
+namespace corsika {
+
+ inline DummyOutputManager::DummyOutputManager() {}
+
+ inline DummyOutputManager::~DummyOutputManager() {}
+
+ template <typename TOutput>
+ inline void DummyOutputManager::add(std::string const&, TOutput&) {}
+
+ inline void DummyOutputManager::startOfLibrary() {}
+
+ inline void DummyOutputManager::startOfShower() {}
+
+ inline void DummyOutputManager::endOfShower() {}
+
+ inline void DummyOutputManager::endOfLibrary() {}
+
+} // namespace corsika
diff --git a/corsika/detail/output/OutputManager.inl b/corsika/detail/output/OutputManager.inl
new file mode 100644
index 0000000000000000000000000000000000000000..6a4ec62045acddbf9028f27d9991239c63012553
--- /dev/null
+++ b/corsika/detail/output/OutputManager.inl
@@ -0,0 +1,238 @@
+/*
+ * (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 <algorithm>
+#include <fstream>
+#include <functional>
+
+#include <iomanip>
+#include <ctime>
+#include <sstream>
+
+#include <boost/filesystem.hpp>
+
+#include <fmt/core.h>
+#include <fmt/chrono.h>
+
+namespace corsika {
+
+ inline void OutputManager::writeNode(YAML::Node const& node,
+ boost::filesystem::path const& path) const {
+
+ // construct a YAML emitter for this config file
+ YAML::Emitter out;
+
+ // and write the node to the output
+ out << node;
+
+ // open the output file - this is <output name>.yaml
+ boost::filesystem::ofstream file(path);
+
+ // dump the YAML to the file
+ file << out.c_str() << std::endl;
+ }
+
+ inline void OutputManager::writeTopLevelConfig() const {
+
+ YAML::Node config;
+
+ // some basic info
+ config["name"] = name_; // the simulation name
+ config["creator"] = "CORSIKA8"; // a tag to identify C8 libraries
+ config["version"] = "8.0.0-prealpha"; // the current version
+
+ // write the node to a file
+ writeNode(config, root_ / ("config.yaml"));
+ }
+
+ inline void OutputManager::writeTopLevelSummary() const {
+
+ YAML::Node config;
+
+ // the total number of showers contained in the library
+ config["showers"] = count_;
+
+ // this next section handles writing some time and duration information
+
+ // create a quick lambda function to convert a time-instance to a string
+ auto timeToString = [&](auto const time) -> std::string {
+ // ISO 8601 time format
+ auto format{"%FT%T%z"};
+
+ // convert the clock to a time_t
+ auto time_tc{std::chrono::system_clock::to_time_t(time)};
+
+ // create the string and push the time onto it
+ std::ostringstream oss;
+ oss << std::put_time(std::localtime(&time_tc), format);
+
+ return oss.str();
+ };
+
+ auto end_time{std::chrono::system_clock::now()};
+
+ // now let's construct an estimate of the runtime
+ auto runtime{end_time - start_time};
+
+ // add the time and duration info
+ config["start time"] = timeToString(start_time);
+ config["end time"] = timeToString(end_time);
+ config["runtime"] = fmt::format("{:%H:%M:%S}", runtime);
+
+ // write the node to a file
+ writeNode(config, root_ / ("summary.yaml"));
+ }
+
+ inline void OutputManager::initOutput(std::string const& name) const {
+ // construct the path to this directory
+ auto const path{root_ / name};
+
+ // create the directory for this process.
+ boost::filesystem::create_directory(path);
+
+ // get the config for this output
+ auto config = outputs_.at(name).get().getConfig();
+
+ // and assign the name for this output
+ config["name"] = name;
+
+ // write the config for this output to the file
+ writeNode(config, path / "config.yaml");
+ }
+
+ inline OutputManager::OutputManager(
+ std::string const& name,
+ boost::filesystem::path const& dir = boost::filesystem::current_path())
+ : name_(name)
+ , root_(dir / name) {
+
+ // check if this directory already exists
+ if (boost::filesystem::exists(root_)) {
+ logger->error("Output directory '{}' already exists! Do not overwrite!.",
+ root_.string());
+ throw std::runtime_error("Output directory already exists.");
+ }
+
+ // construct the directory for this library
+ boost::filesystem::create_directories(root_);
+
+ // write the top level config file
+ writeTopLevelConfig();
+ }
+
+ inline OutputManager::~OutputManager() {
+
+ if (state_ == OutputState::ShowerInProgress) {
+ // if this the destructor is called before the shower has been explicitly
+ // ended, print a warning and end the shower before continuing.
+ logger->warn(
+ "OutputManager was destroyed before endOfShower() called."
+ " The last shower in this libray may be incomplete.");
+ endOfShower();
+ }
+
+ // write the top level summary file (summary.yaml)
+ writeTopLevelSummary();
+
+ // if we are being destructed but EndOfLibrary() has not been called,
+ // make sure that we gracefully close all the outputs. This is a supported
+ // method of operation so we don't issue a warning here
+ if (state_ == OutputState::LibraryReady) { endOfLibrary(); }
+ }
+
+ template <typename TOutput>
+ inline void OutputManager::add(std::string const& name, TOutput& output) {
+
+ // check if that name is already in the map
+ if (outputs_.count(name) > 0) {
+ logger->error("'{}' is already registered. All outputs must have unique names.",
+ name);
+ throw std::runtime_error("Output already exists. Do not overwrite!");
+ }
+
+ // if we get here, the name is not already in the map
+ // so we create the output and register it into the map
+ outputs_.insert(std::make_pair(name, std::ref(output)));
+
+ // and initialize this output
+ initOutput(name);
+ }
+
+ inline void OutputManager::startOfLibrary() {
+
+ // this is only valid when we haven't started a library
+ // or have already finished a library
+ if (!(state_ == OutputState::NoInit || state_ == OutputState::LibraryFinished)) {
+
+ throw std::runtime_error("startOfLibrary() called in invalid state.");
+ }
+
+ // we now forward this signal to all of our outputs
+ for (auto& [name, output] : outputs_) {
+
+ // construct the path to this output subdirectory
+ auto const path{root_ / name};
+
+ // and start the library
+ output.get().startOfLibrary(path);
+ }
+
+ // we have now started running
+ state_ = OutputState::LibraryReady;
+ }
+
+ inline void OutputManager::startOfShower() {
+
+ // if this is called and we still in the "no init" state, then
+ // this is the first shower in the library so make sure we start it
+ if (state_ == OutputState::NoInit) { startOfLibrary(); }
+
+ // now start the event for all the outputs
+ for (auto& [name, output] : outputs_) { output.get().startOfShower(); }
+
+ // increment our shower count
+ ++count_;
+
+ // and transition to the in progress state
+ state_ = OutputState::ShowerInProgress;
+ }
+
+ inline void OutputManager::endOfShower() {
+
+ for (auto& [name, output] : outputs_) { output.get().endOfShower(); }
+
+ // switch back to the initialized state
+ state_ = OutputState::LibraryReady;
+ }
+
+ inline void OutputManager::endOfLibrary() {
+
+ // we can only call endOfLibrary when we have already started
+ if (state_ == OutputState::NoInit) {
+ throw std::runtime_error("endOfLibrary() called in invalid state.");
+ }
+
+ // write the summary for each output and forward the endOfLibrary call()
+ for (auto& [name, output] : outputs_) {
+
+ // we get the summary for each output as a YAML node
+ auto summary{outputs_.at(name).get().getSummary()};
+
+ // write the summary for this output to the file
+ writeNode(summary, root_ / name / "summary.yaml");
+
+ // and forward the end of library call
+ output.get().endOfLibrary();
+ }
+
+ // and the library has finished
+ state_ = OutputState::LibraryFinished;
+ }
+
+} // namespace corsika
diff --git a/corsika/detail/output/ParquetStreamer.inl b/corsika/detail/output/ParquetStreamer.inl
new file mode 100644
index 0000000000000000000000000000000000000000..d0bb2f9e01a42fef34c00fd2d41c189e60f252e2
--- /dev/null
+++ b/corsika/detail/output/ParquetStreamer.inl
@@ -0,0 +1,71 @@
+/*
+ * (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
+
+namespace corsika {
+
+ inline ParquetStreamer::ParquetStreamer()
+ : isInit_(false) {}
+
+ inline void ParquetStreamer::initStreamer(std::string const& filepath) {
+
+ // open the file and connect it to our pointer
+ PARQUET_ASSIGN_OR_THROW(outfile_, arrow::io::FileOutputStream::Open(filepath));
+
+ // the default builder settings
+ builder_.created_by("CORSIKA8");
+
+ // add run and event tags to the file
+ addField("shower", parquet::Repetition::REQUIRED, parquet::Type::INT32,
+ parquet::ConvertedType::INT_32);
+ }
+
+ template <typename... TArgs>
+ inline void ParquetStreamer::addField(TArgs&&... args) {
+ fields_.push_back(parquet::schema::PrimitiveNode::Make(args...));
+ }
+
+ inline void ParquetStreamer::enableCompression(int const /*level*/) {
+ // builder_.compression(parquet::Compression::ZSTD);
+ // builder_.compression_level(level);
+ }
+
+ inline void ParquetStreamer::buildStreamer() {
+
+ // build the top level schema
+ schema_ = std::static_pointer_cast<parquet::schema::GroupNode>(
+ parquet::schema::GroupNode::Make("schema", parquet::Repetition::REQUIRED,
+ fields_));
+
+ // and build the writer
+ writer_ = std::make_shared<parquet::StreamWriter>(
+ parquet::ParquetFileWriter::Open(outfile_, schema_, builder_.build()));
+
+ // only now this object is ready to stream
+ isInit_ = true;
+ }
+
+ inline void ParquetStreamer::closeStreamer() {
+ writer_.reset();
+ [[maybe_unused]] auto status = outfile_->Close();
+ isInit_ = false;
+ }
+
+ inline std::shared_ptr<parquet::StreamWriter> ParquetStreamer::getWriter() {
+ if (!isInit()) {
+ throw std::runtime_error(
+ "ParquetStreamer not initialized. Either 1) add the "
+ "corresponding module to "
+ "the OutputManager, or 2) declare the module to write no output using "
+ "NoOutput.");
+ }
+ return writer_;
+ }
+
+} // namespace corsika
diff --git a/corsika/framework/core/Cascade.hpp b/corsika/framework/core/Cascade.hpp
index b9ee5361bdd00301a1ffe2f2c4ecd75563379ec2..1eadbce1160a405f6cc1738d62997a20923aab3a 100644
--- a/corsika/framework/core/Cascade.hpp
+++ b/corsika/framework/core/Cascade.hpp
@@ -57,7 +57,7 @@ namespace corsika {
*
*
*/
- template <typename TTracking, typename TProcessList, typename TStack,
+ template <typename TTracking, typename TProcessList, typename TOutput, typename TStack,
/*
TStackView is needed as explicit template parameter because
of issue 161 and the
@@ -84,10 +84,11 @@ namespace corsika {
~Cascade() = default;
Cascade& operator=(Cascade const&) = default;
Cascade(Environment<MediumInterface> const& env, TTracking& tr, TProcessList& pl,
- TStack& stack)
+ TOutput& out, TStack& stack)
: environment_(env)
, tracking_(tr)
, sequence_(pl)
+ , output_(out)
, stack_(stack) {
CORSIKA_LOG_INFO(c8_ascii_);
CORSIKA_LOG_INFO("Tracking algorithm: {} (version {})", TTracking::getName(),
@@ -139,6 +140,7 @@ namespace corsika {
Environment<MediumInterface> const& environment_;
TTracking& tracking_;
TProcessList& sequence_;
+ TOutput& output_;
TStack& stack_;
default_prng_type& rng_ = RNGManager::getInstance().getRandomStream("cascade");
unsigned int count_ = 0;
diff --git a/corsika/modules/ObservationPlane.hpp b/corsika/modules/ObservationPlane.hpp
index f7ccba999c531333dbe8dedf9b7550c97f09ba19..613708d1090b19f6cb314f09f62e837ee2e55e20 100644
--- a/corsika/modules/ObservationPlane.hpp
+++ b/corsika/modules/ObservationPlane.hpp
@@ -13,9 +13,7 @@
#include <corsika/framework/process/ContinuousProcess.hpp>
#include <corsika/setup/SetupStack.hpp>
#include <corsika/setup/SetupTrajectory.hpp>
-
-#include <fstream>
-#include <string>
+#include <corsika/modules/writers/ObservationPlaneWriterParquet.hpp>
namespace corsika {
@@ -34,11 +32,12 @@ namespace corsika {
small gap in between the two plane in such a scenario, or develop
another more specialized output class.
*/
- class ObservationPlane : public ContinuousProcess<ObservationPlane> {
+ template <typename TOutputWriter = ObservationPlaneWriterParquet>
+ class ObservationPlane : public ContinuousProcess<ObservationPlane<TOutputWriter>>,
+ public TOutputWriter {
public:
- ObservationPlane(Plane const&, DirectionVector const&, std::string const&,
- bool = true);
+ ObservationPlane(Plane const&, DirectionVector const&, bool = true);
ProcessReturn doContinuous(corsika::setup::Stack::particle_type& vParticle,
corsika::setup::Trajectory& vTrajectory,
@@ -50,10 +49,10 @@ namespace corsika {
void showResults() const;
void reset();
HEPEnergyType getEnergyGround() const { return energy_ground_; }
+ YAML::Node getConfig() const;
private:
Plane const plane_;
- std::ofstream outputStream_;
bool const deleteOnHit_;
HEPEnergyType energy_ground_;
unsigned int count_ground_;
diff --git a/corsika/modules/TrackWriter.hpp b/corsika/modules/TrackWriter.hpp
index 8b9b0efe7faca74ade9f9c55ea463186bc111e98..79e1dfc91bbb52d08300699fc3e6ee1e098feb8a 100644
--- a/corsika/modules/TrackWriter.hpp
+++ b/corsika/modules/TrackWriter.hpp
@@ -8,19 +8,17 @@
#pragma once
-#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/process/ContinuousProcess.hpp>
-
-#include <fstream>
-#include <string>
+#include <corsika/modules/writers/TrackWriterParquet.hpp>
namespace corsika {
- class TrackWriter : public ContinuousProcess<TrackWriter> {
+ template <typename TOutputWriter = TrackWriterParquet>
+ class TrackWriter : public ContinuousProcess<TrackWriter<TOutputWriter>>,
+ public TOutputWriter {
public:
- TrackWriter(std::string const& filename);
- ~TrackWriter();
+ TrackWriter();
template <typename TParticle, typename TTrack>
ProcessReturn doContinuous(TParticle const&, TTrack const&, bool const limitFlag);
@@ -28,12 +26,7 @@ namespace corsika {
template <typename TParticle, typename TTrack>
LengthType getMaxStepLength(TParticle const&, TTrack const&);
- private:
- std::string const filename_;
- std::ofstream file_;
-
- int width_ = 14;
- int precision_ = 6;
+ YAML::Node getConfig() const;
};
} // namespace corsika
diff --git a/corsika/modules/writers/ObservationPlaneWriterParquet.hpp b/corsika/modules/writers/ObservationPlaneWriterParquet.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..2a5089c26f129adb6a5281b8474c509096de3ab6
--- /dev/null
+++ b/corsika/modules/writers/ObservationPlaneWriterParquet.hpp
@@ -0,0 +1,59 @@
+/*
+ * (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/output/BaseOutput.hpp>
+#include <corsika/output/ParquetStreamer.hpp>
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+
+namespace corsika {
+
+ class ObservationPlaneWriterParquet : public BaseOutput {
+
+ ParquetStreamer output_; ///< The primary output file.
+
+ public:
+ /**
+ * Construct an ObservationPlane.
+ *
+ * @param name The name of this output.
+ */
+ ObservationPlaneWriterParquet();
+
+ /**
+ * Called at the start of each library.
+ */
+ void startOfLibrary(boost::filesystem::path const& directory) final override;
+
+ /**
+ * Called at the end of each shower.
+ */
+ void endOfShower() final override;
+
+ /**
+ * Called at the end of each library.
+ *
+ * This must also increment the run number since we override
+ * the default behaviour of BaseOutput.
+ */
+ void endOfLibrary() final override;
+
+ protected:
+ /**
+ * Write a particle to the file.
+ */
+ void write(Code const& pid, units::si::HEPEnergyType const& energy,
+ units::si::LengthType const& x, units::si::LengthType const& y);
+
+ }; // class ObservationPlaneWriterParquet
+
+} // namespace corsika
+
+#include <corsika/detail/modules/writers/ObservationPlaneWriterParquet.inl>
diff --git a/corsika/modules/writers/TrackWriterParquet.hpp b/corsika/modules/writers/TrackWriterParquet.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..7caf90558b2ae0fc5c1cc13ca38b7d009e70b335
--- /dev/null
+++ b/corsika/modules/writers/TrackWriterParquet.hpp
@@ -0,0 +1,62 @@
+/*
+ * (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/output/BaseOutput.hpp>
+#include <corsika/output/ParquetStreamer.hpp>
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/QuantityVector.hpp>
+
+namespace corsika {
+
+ class TrackWriterParquet : public BaseOutput {
+
+ public:
+ /**
+ * Construct a new writer.
+ *
+ * @param name The name of this output.
+ */
+ TrackWriterParquet();
+
+ /**
+ * Called at the start of each library.
+ */
+ void startOfLibrary(boost::filesystem::path const& directory) final override;
+
+ /**
+ * Called at the end of each shower.
+ */
+ void endOfShower() final override;
+
+ /**
+ * Called at the end of each library.
+ *
+ * This must also increment the run number since we override
+ * the default behaviour of BaseOutput.
+ */
+ void endOfLibrary() final override;
+
+ protected:
+ /**
+ * Write a track to the file.
+ */
+ void write(Code const& pid, units::si::HEPEnergyType const& energy,
+ QuantityVector<length_d> const& start,
+ QuantityVector<length_d> const& end);
+
+ private:
+ ParquetStreamer output_; ///< The primary output file.
+
+ }; // class TrackWriterParquet
+
+} // namespace corsika
+
+#include <corsika/detail/modules/writers/TrackWriterParquet.inl>
diff --git a/corsika/output/BaseOutput.hpp b/corsika/output/BaseOutput.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..0024c5b7a918e4705fbbadd3bdcdd59dfcab0449
--- /dev/null
+++ b/corsika/output/BaseOutput.hpp
@@ -0,0 +1,75 @@
+/*
+ * (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 <boost/filesystem.hpp>
+
+#include <yaml-cpp/yaml.h>
+
+namespace corsika {
+
+ /**
+ * This is the base class for all outputs so that they
+ * can be stored in homogeneous containers.
+ */
+ class BaseOutput {
+
+ protected:
+ BaseOutput();
+
+ public:
+ /**
+ * Called at the start of each run.
+ */
+ virtual void startOfLibrary(boost::filesystem::path const& directory) = 0;
+
+ /**
+ * Called at the start of each event/shower.
+ */
+ virtual void startOfShower() {}
+
+ /**
+ * Called at the end of each event/shower.
+ */
+ virtual void endOfShower() = 0;
+
+ /**
+ * Called at the end of each run.
+ */
+ virtual void endOfLibrary() = 0;
+
+ /**
+ * Get the configuration of this output.
+ */
+ virtual YAML::Node getConfig() const = 0;
+
+ /**
+ * Get any summary information for the entire library.
+ */
+ virtual YAML::Node getSummary() { return YAML::Node(); }
+
+ /**
+ * Flag to indicate readiness.
+ */
+ bool isInit() const { return is_init_; }
+
+ protected:
+ /**
+ * Set init flag.
+ */
+ void setInit(bool const v) { is_init_ = v; }
+
+ int shower_{0}; ///< The current event number.
+
+ private:
+ bool is_init_{false}; ///< flag to indicate readiness
+ };
+
+} // namespace corsika
+
+#include <corsika/detail/output/BaseOutput.inl>
diff --git a/corsika/output/DummyOutputManager.hpp b/corsika/output/DummyOutputManager.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..fc5cf4647d5daa227fceceda5961bb051a840058
--- /dev/null
+++ b/corsika/output/DummyOutputManager.hpp
@@ -0,0 +1,70 @@
+/*
+ * (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
+
+namespace corsika {
+
+ /*!
+ * An output manager that does nothing.
+ */
+ class DummyOutputManager final {
+
+ public:
+ /**
+ * 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.
+ */
+ DummyOutputManager();
+
+ /**
+ * Handle graceful closure of the outputs upon destruction.
+ */
+ ~DummyOutputManager();
+
+ /**
+ * Register an existing output to this manager.
+ *
+ * @param name The unique name of this output.
+ * @param args... These are perfect forwarded to the
+ * constructor of the output.
+ */
+ template <typename TOutput>
+ void add(std::string const& name, TOutput& output);
+
+ /**
+ * Called at the start of each library.
+ *
+ * This iteratively calls startOfLibrary on each registered output.
+ */
+ void startOfLibrary();
+
+ /**
+ * Called at the start of each event/shower.
+ * This iteratively calls startOfEvent on each registered output.
+ */
+ void startOfShower();
+
+ /**
+ * Called at the end of each event/shower.
+ * This iteratively calls endOfEvent on each registered output.
+ */
+ void endOfShower();
+
+ /**
+ * Called at the end of each library.
+ * This iteratively calls endOfLibrary on each registered output.
+ */
+ void endOfLibrary();
+
+ }; // class DummyOutputManager
+
+} // namespace corsika
+
+#include <corsika/detail/output/DummyOutputManager.inl>
diff --git a/corsika/output/NoOutput.hpp b/corsika/output/NoOutput.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..289c19ad720f9f03d621d3b50ac031c40cc794cf
--- /dev/null
+++ b/corsika/output/NoOutput.hpp
@@ -0,0 +1,62 @@
+/*
+ * (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/core/ParticleProperties.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+
+namespace corsika {
+
+ /**
+ * This is the base class for all outputs so that they
+ * can be stored in homogeneous containers.
+ */
+ class NoOutput : public BaseOutput {
+
+ protected:
+ NoOutput() {}
+
+ public:
+ /**
+ * Called at the start of each run.
+ */
+ void startOfLibrary(boost::filesystem::path const&) final override {}
+
+ /**
+ * Called at the start of each event/shower.
+ */
+ void startOfShower() final override {}
+
+ /**
+ * Called at the end of each event/shower.
+ */
+ void endOfShower() final override {}
+
+ /**
+ * Called at the end of each run.
+ */
+ void endOfLibrary() final override {}
+
+ /**
+ * Get the configuration of this output.
+ */
+ YAML::Node getConfig() const { return YAML::Node(); };
+
+ /**
+ * Get any summary information for the entire library.
+ */
+ YAML::Node getSummary() final override { return YAML::Node(); };
+
+ protected:
+ void write(Code const&, units::si::HEPEnergyType const&, units::si::LengthType const&,
+ units::si::LengthType const&) {}
+ };
+
+} // namespace corsika
+
+#include <corsika/detail/output/BaseOutput.inl>
diff --git a/corsika/output/OutputManager.hpp b/corsika/output/OutputManager.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..2582220520454d028b660ca06f5ac93ed8a18b73
--- /dev/null
+++ b/corsika/output/OutputManager.hpp
@@ -0,0 +1,118 @@
+/*
+ * (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 <chrono>
+#include <string>
+#include <boost/filesystem.hpp>
+#include <corsika/output/BaseOutput.hpp>
+#include <corsika/framework/core/Logging.hpp>
+
+namespace corsika {
+
+ /*!
+ * Manages CORSIKA 8 output streams.
+ */
+ class OutputManager final {
+
+ /**
+ * Indicates the current state of this manager.
+ */
+ enum class OutputState {
+ NoInit,
+ LibraryReady,
+ ShowerInProgress,
+ LibraryFinished,
+ };
+
+ OutputState state_{OutputState::NoInit}; ///< The current state of this manager.
+ std::string const name_; ///< The name of this simulation file.
+ boost::filesystem::path const root_; ///< The top-level directory for the output.
+ int count_{0}; ///< The current ID of this shower.
+ std::chrono::time_point<std::chrono::system_clock> const start_time{
+ std::chrono::system_clock::now()}; ///< The time the manager is created.
+ inline static auto logger{get_logger("output")}; ///< A custom logger.
+ /**
+ * The outputs that have been registered with us.
+ */
+ std::map<std::string, std::reference_wrapper<BaseOutput>> outputs_;
+
+ /**
+ * Write a YAML-node to a file.
+ */
+ void writeNode(YAML::Node const& node, boost::filesystem::path const& path) const;
+
+ /**
+ * Write the top-level config of this simulation.
+ */
+ void writeTopLevelConfig() const;
+
+ /**
+ * Initialize the "registered" output with a given name.
+ */
+ void initOutput(std::string const& name) const;
+
+ /**
+ * Write the top-level summary of this library.
+ */
+ void writeTopLevelSummary() const;
+
+ public:
+ /**
+ * 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.
+ */
+ OutputManager(std::string const& name, boost::filesystem::path const& dir);
+
+ /**
+ * Handle graceful closure of the outputs upon destruction.
+ */
+ ~OutputManager();
+
+ /**
+ * Register an existing output to this manager.
+ *
+ * @param name The unique name of this output.
+ * @param args... These are perfect forwarded to the
+ * constructor of the output.
+ */
+ template <typename TOutput>
+ void add(std::string const& name, TOutput& output);
+
+ /**
+ * Called at the start of each library.
+ *
+ * This iteratively calls startOfLibrary on each registered output.
+ */
+ void startOfLibrary();
+
+ /**
+ * Called at the start of each event/shower.
+ * This iteratively calls startOfEvent on each registered output.
+ */
+ void startOfShower();
+
+ /**
+ * Called at the end of each event/shower.
+ * This iteratively calls endOfEvent on each registered output.
+ */
+ void endOfShower();
+
+ /**
+ * Called at the end of each library.
+ * This iteratively calls endOfLibrary on each registered output.
+ */
+ void endOfLibrary();
+
+ }; // class OutputManager
+
+} // namespace corsika
+
+#include <corsika/detail/output/OutputManager.inl>
diff --git a/corsika/output/ParquetStreamer.hpp b/corsika/output/ParquetStreamer.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..dcf1f547f72338b45459583a591150487de937c8
--- /dev/null
+++ b/corsika/output/ParquetStreamer.hpp
@@ -0,0 +1,83 @@
+
+/*
+ * (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 <string>
+
+// clang-format-off, NOTE: the order of these includes is *important*
+// you will get unhelpful compiler errors about unknown
+// operator definitions if these are reordered
+#include <parquet/stream_writer.h>
+#include <parquet/arrow/schema.h>
+#include <arrow/io/file.h>
+// clang-format-on
+
+namespace corsika {
+
+ /**
+ * This class automates the construction of simple tabular
+ * Parquet files using the parquet::StreamWriter.
+ */
+ class ParquetStreamer {
+
+ public:
+ /**
+ * ParquetStreamer's take no constructor arguments.
+ */
+ ParquetStreamer();
+
+ /**
+ * Initialize the streamer to write to a given file.
+ */
+ void initStreamer(std::string const& filepath);
+
+ /**
+ * Add a field to this streamer.
+ */
+ template <typename... TArgs>
+ void addField(TArgs&&... args);
+
+ /**
+ * Enable compression for this streamer.
+ */
+ void enableCompression(int const level = 3);
+
+ /**
+ * Finalize the streamer construction.
+ */
+ void buildStreamer();
+
+ /**
+ * Finish writing this stream.
+ */
+ void closeStreamer();
+
+ /**
+ * Return a reference to the underlying writer.
+ */
+ std::shared_ptr<parquet::StreamWriter> getWriter();
+
+ /**
+ * @return status of streamer
+ */
+ bool isInit() const { return isInit_; }
+
+ private:
+ bool isInit_ = false; ///< flag to handle state of writer
+ parquet::WriterProperties::Builder builder_; ///< The writer properties builder.
+ parquet::schema::NodeVector fields_; ///< The fields in this file.
+ std::shared_ptr<parquet::schema::GroupNode> schema_; ///< The schema for this file.
+ std::shared_ptr<arrow::io::FileOutputStream> outfile_; ///< The output file.
+ std::shared_ptr<parquet::StreamWriter> writer_; ///< The stream writer to 'outfile'
+
+ }; // class ParquetStreamer
+} // namespace corsika
+
+#include <corsika/detail/output/ParquetStreamer.inl>
diff --git a/examples/boundary_example.cpp b/examples/boundary_example.cpp
index bca45481a3841f690e80ea5d936063a4ffe91acf..ebb202597bc8f7d09e88b815e89e5f55cf5a3225 100644
--- a/examples/boundary_example.cpp
+++ b/examples/boundary_example.cpp
@@ -14,6 +14,8 @@
#include <corsika/framework/utility/CorsikaFenv.hpp>
#include <corsika/framework/core/Logging.hpp>
+#include <corsika/output/OutputManager.hpp>
+
#include <corsika/setup/SetupEnvironment.hpp>
#include <corsika/setup/SetupStack.hpp>
#include <corsika/setup/SetupTrajectory.hpp>
@@ -42,6 +44,7 @@
#include <iostream>
#include <limits>
#include <typeinfo>
+#include <fstream>
using namespace corsika;
using namespace std;
@@ -112,6 +115,8 @@ int main() {
world->addChild(std::move(target));
universe.addChild(std::move(world));
+ OutputManager output("boundary_outputs");
+
// setup processes, decays and interactions
setup::Tracking tracking;
@@ -121,7 +126,9 @@ int main() {
ParticleCut cut(50_GeV, true, true);
- TrackWriter trackWriter("boundary_tracks.dat");
+ TrackWriter trackWriter;
+ output.add("tracks", trackWriter); // register TrackWriter
+
MyBoundaryCrossingProcess<true> boundaryCrossing("crossings.dat");
// assemble all processes into an ordered process list
@@ -164,7 +171,7 @@ int main() {
}
// define air shower object, run simulation
- Cascade EAS(env, tracking, sequence, stack);
+ Cascade EAS(env, tracking, sequence, output, stack);
EAS.run();
@@ -174,4 +181,6 @@ int main() {
(cut.getCutEnergy() + cut.getInvEnergy() + cut.getEmEnergy());
CORSIKA_LOG_INFO("Total energy (GeV): {} relative difference (%): {}", Efinal / 1_GeV,
(Efinal / E0 - 1.) * 100);
+
+ output.endOfLibrary();
}
diff --git a/examples/cascade_example.cpp b/examples/cascade_example.cpp
index b9195515c189d992adaedba6f59cf99496ab44e5..f7758db0a3163521274ee5b0b33e21cc478c4405 100644
--- a/examples/cascade_example.cpp
+++ b/examples/cascade_example.cpp
@@ -11,10 +11,11 @@
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/random/RNGManager.hpp>
#include <corsika/framework/geometry/Sphere.hpp>
-
#include <corsika/framework/utility/CorsikaFenv.hpp>
#include <corsika/framework/core/Logging.hpp>
+#include <corsika/output/OutputManager.hpp>
+
#include <corsika/media/Environment.hpp>
#include <corsika/media/HomogeneousMedium.hpp>
#include <corsika/media/NuclearComposition.hpp>
@@ -54,7 +55,7 @@ using namespace std;
//
int main() {
- logging::set_level(logging::level::trace);
+ logging::set_level(logging::level::info);
std::cout << "cascade_example" << std::endl;
@@ -106,6 +107,8 @@ int main() {
rootCS, 0_m, 0_m,
height_atmosphere); // this is the CORSIKA 7 start of atmosphere/universe
+ OutputManager output("cascade_outputs");
+
ShowerAxis const showerAxis{injectionPos, Vector{rootCS, 0_m, 0_m, -100_km}, env};
{
@@ -139,7 +142,9 @@ int main() {
// cascade with only HE model ==> HE cut
ParticleCut cut(80_GeV, true, true);
- TrackWriter trackWriter("tracks.dat");
+ TrackWriter trackWriter;
+ output.add("tracks", trackWriter); // register TrackWriter
+
BetheBlochPDG eLoss{showerAxis};
// assemble all processes into an ordered process list
@@ -147,7 +152,7 @@ int main() {
make_sequence(stackInspect, sibyll, sibyllNuc, decay, eLoss, cut, trackWriter);
// define air shower object, run simulation
- Cascade EAS(env, tracking, sequence, stack);
+ Cascade EAS(env, tracking, sequence, output, stack);
EAS.run();
@@ -161,4 +166,6 @@ int main() {
cout << "total dEdX energy (GeV): " << eLoss.getTotal() / 1_GeV << endl
<< "relative difference (%): " << eLoss.getTotal() / E0 * 100 << endl;
cut.reset();
+
+ output.endOfLibrary();
}
diff --git a/examples/cascade_proton_example.cpp b/examples/cascade_proton_example.cpp
index 4b24f3b0148045a2db97546c37c6a2916e3c1d1e..a02c914e4ad1d554160d36745e5d00a2ee6d981e 100644
--- a/examples/cascade_proton_example.cpp
+++ b/examples/cascade_proton_example.cpp
@@ -11,10 +11,11 @@
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/random/RNGManager.hpp>
#include <corsika/framework/geometry/Sphere.hpp>
-
#include <corsika/framework/utility/CorsikaFenv.hpp>
#include <corsika/framework/core/Logging.hpp>
+#include <corsika/output/OutputManager.hpp>
+
#include <corsika/media/Environment.hpp>
#include <corsika/media/HomogeneousMedium.hpp>
#include <corsika/media/NuclearComposition.hpp>
@@ -66,6 +67,8 @@ int main() {
// initialize random number sequence(s)
RNGManager::getInstance().registerRandomStream("cascade");
+ OutputManager output("cascade_proton_outputs");
+
// setup environment, geometry
using EnvType = setup::Environment;
EnvType env;
@@ -115,7 +118,7 @@ int main() {
// setup processes, decays and interactions
setup::Tracking tracking;
- StackInspector<setup::Stack> stackInspect(1, true, E0);
+ StackInspector<setup::Stack> stackInspect(1000, true, E0);
RNGManager::getInstance().registerRandomStream("sibyll");
RNGManager::getInstance().registerRandomStream("pythia");
@@ -130,7 +133,9 @@ int main() {
// HadronicElasticModel::HadronicElasticInteraction
// hadronicElastic(env);
- TrackWriter trackWriter("tracks.dat");
+ TrackWriter trackWriter;
+ output.add("tracks", trackWriter); // register TrackWriter
+
ShowerAxis const showerAxis{injectionPos, Vector{rootCS, 0_m, 0_m, -100_km}, env};
BetheBlochPDG eLoss{showerAxis};
@@ -141,7 +146,7 @@ int main() {
auto sequence = make_sequence(pythia, decay, eLoss, cut, trackWriter, stackInspect);
// define air shower object, run simulation
- Cascade EAS(env, tracking, sequence, stack);
+ Cascade EAS(env, tracking, sequence, output, stack);
EAS.run();
cout << "Result: E0=" << E0 / 1_GeV << endl;
@@ -150,4 +155,6 @@ int main() {
cut.getCutEnergy() + cut.getInvEnergy() + cut.getEmEnergy();
cout << "total energy (GeV): " << Efinal / 1_GeV << endl
<< "relative difference (%): " << (Efinal / E0 - 1.) * 100 << endl;
+
+ output.endOfLibrary();
}
diff --git a/examples/em_shower.cpp b/examples/em_shower.cpp
index e4d8b3cc432a1852f96dbababca6a0c024fc037b..9634b18703905b104b7c882139972d36125cc0f4 100644
--- a/examples/em_shower.cpp
+++ b/examples/em_shower.cpp
@@ -18,6 +18,8 @@
#include <corsika/framework/process/InteractionCounter.hpp>
#include <corsika/framework/core/Logging.hpp>
+#include <corsika/output/OutputManager.hpp>
+
#include <corsika/media/Environment.hpp>
#include <corsika/media/LayeredSphericalAtmosphereBuilder.hpp>
#include <corsika/media/NuclearComposition.hpp>
@@ -66,7 +68,6 @@ using MyExtraEnv = MediumPropertyModel<UniformMagneticField<T>>;
int main(int argc, char** argv) {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
if (argc != 2) {
std::cerr << "usage: em_shower <energy/GeV>" << std::endl;
@@ -137,6 +138,7 @@ int main(int argc, char** argv) {
std::cout << "shower axis length: " << (showerCore - injectionPos).getNorm() * 1.02
<< std::endl;
+ OutputManager output("em_shower_outputs");
ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.02, env};
// setup processes, decays and interactions
@@ -146,7 +148,8 @@ int main(int argc, char** argv) {
corsika::proposal::ContinuousProcess emContinuous(env);
InteractionCounter emCascadeCounted(emCascade);
- TrackWriter trackWriter("tracks.dat");
+ TrackWriter trackWriter;
+ output.add("tracks", trackWriter); // register TrackWriter
// long. profile; columns for gamma, e+, e- still need to be added
LongitudinalProfile longprof{showerAxis};
@@ -154,12 +157,13 @@ int main(int argc, char** argv) {
Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.}));
ObservationPlane observationLevel(obsPlane, DirectionVector(rootCS, {1., 0., 0.}),
"particles.dat");
+ output.add("obsplane", observationLevel);
auto sequence = make_sequence(emCascadeCounted, emContinuous, longprof, cut,
observationLevel, trackWriter);
// define air shower object, run simulation
setup::Tracking tracking;
- Cascade EAS(env, tracking, sequence, stack);
+ Cascade EAS(env, tracking, sequence, output, stack);
// to fix the point of first interaction, uncomment the following two lines:
// EAS.setNodes();
@@ -183,4 +187,6 @@ int main(int argc, char** argv) {
save_hist(hists.labHist(), "inthist_lab_emShower.npz", true);
save_hist(hists.CMSHist(), "inthist_cms_emShower.npz", true);
longprof.save("longprof_emShower.txt");
+
+ output.endOfLibrary();
}
diff --git a/examples/geometry_example.cpp b/examples/geometry_example.cpp
index 7a25708a8d434b8b35e89a8cd23b1e4ad35e0f1f..596e6c1842011d56d8db27dc6e02fe81ab8bdbb9 100644
--- a/examples/geometry_example.cpp
+++ b/examples/geometry_example.cpp
@@ -21,7 +21,6 @@ using namespace corsika;
int main() {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
CORSIKA_LOG_INFO("geometry_example");
diff --git a/examples/helix_example.cpp b/examples/helix_example.cpp
index cf58531b4170aa5b74a3693033b71c700cde102e..e980931b74ace074bf6d0459cf69e6d158f32b91 100644
--- a/examples/helix_example.cpp
+++ b/examples/helix_example.cpp
@@ -21,7 +21,6 @@ using namespace corsika;
int main() {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
CORSIKA_LOG_INFO("helix_example");
diff --git a/examples/hybrid_MC.cpp b/examples/hybrid_MC.cpp
index af379ca89bc2a14ccc706e653ca6e3fa25e7d49c..f1602f5ce24a7359d8e1de500369f70d32177b5f 100644
--- a/examples/hybrid_MC.cpp
+++ b/examples/hybrid_MC.cpp
@@ -25,6 +25,8 @@
#include <corsika/framework/core/Cascade.hpp>
#include <corsika/framework/geometry/PhysicalGeometry.hpp>
+#include <corsika/output/OutputManager.hpp>
+
#include <corsika/media/Environment.hpp>
#include <corsika/media/FlatExponential.hpp>
#include <corsika/media/LayeredSphericalAtmosphereBuilder.hpp>
@@ -86,7 +88,6 @@ using MyExtraEnv = MediumPropertyModel<UniformMagneticField<T>>;
int main(int argc, char** argv) {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
CORSIKA_LOG_INFO("hybrid_MC");
@@ -171,6 +172,7 @@ int main(int argc, char** argv) {
std::cout << "shower axis length: " << (showerCore - injectionPos).getNorm() * 1.02
<< std::endl;
+ OutputManager output("hybrid_MC_outputs");
ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.02, env};
// setup processes, decays and interactions
@@ -219,6 +221,7 @@ int main(int argc, char** argv) {
Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.}));
ObservationPlane observationLevel(obsPlane, DirectionVector(rootCS, {1., 0., 0.}),
"particles.dat");
+ output.add("obsplane", observationLevel);
corsika::urqmd::UrQMD urqmd_model;
InteractionCounter urqmdCounted{urqmd_model};
@@ -240,7 +243,7 @@ int main(int argc, char** argv) {
// define air shower object, run simulation
setup::Tracking tracking;
- Cascade EAS(env, tracking, sequence, stack);
+ Cascade EAS(env, tracking, sequence, output, stack);
// to fix the point of first interaction, uncomment the following two lines:
// EAS.SetNodes();
@@ -272,4 +275,6 @@ int main(int argc, char** argv) {
std::ofstream finish("finished");
finish << "run completed without error" << std::endl;
+
+ output.endOfLibrary();
}
diff --git a/examples/vertical_EAS.cpp b/examples/vertical_EAS.cpp
index e5f6455d3f6782cef10d69c128e38534c50a2d8d..928197970e8c27f08b93d561f33aab1defea5796 100644
--- a/examples/vertical_EAS.cpp
+++ b/examples/vertical_EAS.cpp
@@ -25,6 +25,8 @@
#include <corsika/framework/core/Cascade.hpp>
#include <corsika/framework/geometry/PhysicalGeometry.hpp>
+#include <corsika/output/OutputManager.hpp>
+
#include <corsika/media/Environment.hpp>
#include <corsika/media/FlatExponential.hpp>
#include <corsika/media/HomogeneousMedium.hpp>
@@ -97,7 +99,6 @@ using MyExtraEnv = MediumPropertyModel<UniformMagneticField<T>>;
int main(int argc, char** argv) {
- corsika_logger->set_pattern("[%n:%^%-8l%$] %s:%#: %v");
logging::set_level(logging::level::info);
CORSIKA_LOG_INFO("vertical_EAS");
@@ -134,6 +135,7 @@ int main(int argc, char** argv) {
{{Code::Nitrogen, Code::Oxygen},
{0.7847f, 1.f - 0.7847f}}); // values taken from AIRES manual, Ar removed for now
+ builder.addExponentialLayer(1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 2_km);
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);
builder.addExponentialLayer(1305.5948_g / (1_cm * 1_cm), 636143.04_cm, 40_km);
@@ -141,22 +143,6 @@ int main(int argc, char** argv) {
builder.addLinearLayer(1e9_cm, 112.8_km + constants::EarthRadius::Mean);
builder.assemble(env);
- CORSIKA_LOG_DEBUG(
- "environment setup: universe={}, layer1={}, layer2={}, layer3={}, layer4={}, "
- "layer5={}",
- fmt::ptr(env.getUniverse()->getContainingNode(
- Point(rootCS, {constants::EarthRadius::Mean + 130_km, 0_m, 0_m}))),
- fmt::ptr(env.getUniverse()->getContainingNode(
- Point(rootCS, {constants::EarthRadius::Mean + 110_km, 0_m, 0_m}))),
- fmt::ptr(env.getUniverse()->getContainingNode(
- Point(rootCS, {constants::EarthRadius::Mean + 50_km, 0_m, 0_m}))),
- fmt::ptr(env.getUniverse()->getContainingNode(
- Point(rootCS, {constants::EarthRadius::Mean + 20_km, 0_m, 0_m}))),
- fmt::ptr(env.getUniverse()->getContainingNode(
- Point(rootCS, {constants::EarthRadius::Mean + 5_km, 0_m, 0_m}))),
- fmt::ptr(env.getUniverse()->getContainingNode(
- Point(rootCS, {constants::EarthRadius::Mean + 2_km, 0_m, 0_m}))));
-
// pre-setup particle stack
unsigned short const A = std::stoi(std::string(argv[1]));
Code beamCode;
@@ -214,6 +200,9 @@ int main(int argc, char** argv) {
ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.5, env};
+ // create the output manager that we then register outputs with
+ OutputManager output("vertical_EAS_outputs");
+
// setup processes, decays and interactions
// corsika::qgsjetII::Interaction qgsjet;
@@ -280,8 +269,6 @@ int main(int argc, char** argv) {
string const labHist_file = "inthist_lab_verticalEAS_" + to_string(i_shower) + ".npz";
string const cMSHist_file = "inthist_cms_verticalEAS_" + to_string(i_shower) + ".npz";
string const longprof_file = "longprof_verticalEAS_" + to_string(i_shower) + ".txt";
- string const tracks_file = "tracks_" + to_string(i_shower) + ".dat";
- string const particles_file = "particles_" + to_string(i_shower) + ".dat";
std::cout << std::endl;
std::cout << "Shower " << i_shower << "/" << number_showers << std::endl;
@@ -289,6 +276,46 @@ int main(int argc, char** argv) {
// setup particle stack, and add primary particle
setup::Stack stack;
stack.clear();
+ unsigned short const A = std::stoi(std::string(argv[1]));
+ Code beamCode;
+ HEPEnergyType mass;
+ unsigned short Z = 0;
+ if (A > 0) {
+ beamCode = Code::Nucleus;
+ Z = std::stoi(std::string(argv[2]));
+ mass = get_nucleus_mass(A, Z);
+ } else {
+ int pdg = std::stoi(std::string(argv[2]));
+ beamCode = convert_from_PDG(PDGCode(pdg));
+ mass = get_mass(beamCode);
+ }
+ HEPEnergyType const E0 = 1_GeV * std::stof(std::string(argv[3]));
+ double theta = 0.;
+ auto const thetaRad = theta / 180. * M_PI;
+
+ auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
+ return sqrt((Elab - m) * (Elab + m));
+ };
+ HEPMomentumType P0 = elab2plab(E0, mass);
+ auto momentumComponents = [](double thetaRad, HEPMomentumType ptot) {
+ return std::make_tuple(ptot * sin(thetaRad), 0_eV, -ptot * cos(thetaRad));
+ };
+
+ auto const [px, py, pz] = momentumComponents(thetaRad, P0);
+ auto plab = MomentumVector(rootCS, {px, py, pz});
+ cout << "input particle: " << beamCode << endl;
+ cout << "input angles: theta=" << theta << endl;
+ cout << "input momentum: " << plab.getComponents() / 1_GeV
+ << ", norm = " << plab.getNorm() << endl;
+
+ auto const observationHeight = 0_km + builder.getEarthRadius();
+ auto const injectionHeight = 111.75_km + builder.getEarthRadius();
+ auto const t = (injectionHeight - observationHeight) / cos(thetaRad);
+ Point const showerCore{rootCS, 0_m, 0_m, observationHeight};
+ Point const injectionPos =
+ showerCore + DirectionVector{rootCS, {-sin(thetaRad), 0, cos(thetaRad)}} * t;
+
+ std::cout << "point of injection: " << injectionPos.getCoordinates() << std::endl;
if (A > 1) {
stack.addParticle(std::make_tuple(beamCode, E0, plab, injectionPos, 0_ns, A, Z));
@@ -308,9 +335,64 @@ int main(int argc, char** argv) {
}
}
- TrackWriter trackWriter(tracks_file);
- ObservationPlane observationLevel(obsPlane, DirectionVector(rootCS, {1., 0., 0.}),
- particles_file);
+ // we make the axis much longer than the inj-core distance since the
+ // profile will go beyond the core, depending on zenith angle
+ std::cout << "shower axis length: " << (showerCore - injectionPos).getNorm() * 1.5
+ << std::endl;
+
+ ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.5, env,
+ false};
+
+ // setup processes, decays and interactions
+
+ // corsika::qgsjetII::Interaction qgsjet;
+ corsika::sibyll::Interaction sibyll;
+ InteractionCounter sibyllCounted(sibyll);
+
+ corsika::sibyll::NuclearInteraction sibyllNuc(sibyll, env);
+ InteractionCounter sibyllNucCounted(sibyllNuc);
+
+ corsika::pythia8::Decay decayPythia;
+
+ // use sibyll decay routine for decays of particles unknown to pythia
+ corsika::sibyll::Decay decaySibyll{{
+ Code::N1440Plus,
+ Code::N1440MinusBar,
+ Code::N1440_0,
+ Code::N1440_0Bar,
+ Code::N1710Plus,
+ Code::N1710MinusBar,
+ Code::N1710_0,
+ Code::N1710_0Bar,
+
+ Code::Pi1300Plus,
+ Code::Pi1300Minus,
+ Code::Pi1300_0,
+
+ Code::KStar0_1430_0,
+ Code::KStar0_1430_0Bar,
+ Code::KStar0_1430_Plus,
+ Code::KStar0_1430_MinusBar,
+ }};
+
+ decaySibyll.printDecayConfig();
+
+ ParticleCut cut{60_GeV, true, true};
+ // corsika::proposal::Interaction emCascade(env);
+ // corsika::proposal::ContinuousProcess emContinuous(env);
+ // InteractionCounter emCascadeCounted(emCascade);
+ BetheBlochPDG emContinuous(showerAxis);
+
+ OnShellCheck reset_particle_mass(1.e-3, 1.e-1, false);
+ TrackWriter trackWriter;
+ output.add("tracks", trackWriter); // register TrackWriter
+
+ LongitudinalProfile longprof{showerAxis};
+
+ Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.}));
+ ObservationPlane observationLevel(obsPlane, DirectionVector(rootCS, {1., 0., 0.}));
+ // register the observation plane with the output
+ output.add("obsplane", observationLevel);
auto sequence =
make_sequence(stackInspect, hadronSequence, reset_particle_mass, decaySequence,
@@ -318,7 +400,7 @@ int main(int argc, char** argv) {
// define air shower object, run simulation
setup::Tracking tracking;
- Cascade EAS(env, tracking, sequence, stack);
+ Cascade EAS(env, tracking, sequence, output, stack);
// to fix the point of first interaction, uncomment the following two lines:
// EAS.forceInteraction();
@@ -326,16 +408,16 @@ int main(int argc, char** argv) {
EAS.run();
cut.showResults();
- emContinuous.showResults();
+ // emContinuous.showResults();
observationLevel.showResults();
const HEPEnergyType Efinal = cut.getCutEnergy() + cut.getInvEnergy() +
- cut.getEmEnergy() + emContinuous.getEnergyLost() +
+ cut.getEmEnergy() + // emContinuous.getEnergyLost() +
observationLevel.getEnergyGround();
cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl
<< "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
observationLevel.reset();
cut.reset();
- emContinuous.reset();
+ // emContinuous.reset();
auto const hists = sibyllCounted.getHistogram() + sibyllNucCounted.getHistogram() +
urqmdCounted.getHistogram();
@@ -343,5 +425,7 @@ int main(int argc, char** argv) {
save_hist(hists.labHist(), labHist_file, true);
save_hist(hists.CMSHist(), cMSHist_file, true);
longprof.save(longprof_file);
+
+ output.endOfLibrary();
}
}
diff --git a/python/.gitignore b/python/.gitignore
new file mode 100644
index 0000000000000000000000000000000000000000..7b9e44cc96089e3d033d74eb46965d968bd28d52
--- /dev/null
+++ b/python/.gitignore
@@ -0,0 +1,147 @@
+# Byte-compiled / optimized / DLL files
+__pycache__/
+*.py[cod]
+*$py.class
+
+# ignore any generated output files
+*.npz
+*.dat
+
+# C extensions
+*.so
+
+# Distribution / packaging
+.Python
+build/
+develop-eggs/
+dist/
+downloads/
+eggs/
+.eggs/
+lib/
+lib64/
+parts/
+sdist/
+var/
+wheels/
+pip-wheel-metadata/
+share/python-wheels/
+*.egg-info/
+.installed.cfg
+*.egg
+MANIFEST
+
+# PyInstaller
+# Usually these files are written by a python script from a template
+# before PyInstaller builds the exe, so as to inject date/other infos into it.
+*.manifest
+*.spec
+
+# Installer logs
+pip-log.txt
+pip-delete-this-directory.txt
+
+# Unit test / coverage reports
+htmlcov/
+.tox/
+.nox/
+.coverage
+.coverage.*
+.cache
+nosetests.xml
+coverage.xml
+*.cover
+*.py,cover
+.hypothesis/
+.pytest_cache/
+cover/
+
+# Translations
+*.mo
+*.pot
+
+# Django stuff:
+*.log
+local_settings.py
+db.sqlite3
+db.sqlite3-journal
+
+# Flask stuff:
+instance/
+.webassets-cache
+
+# Scrapy stuff:
+.scrapy
+
+# Sphinx documentation
+docs/_build/
+
+# PyBuilder
+.pybuilder/
+target/
+
+# Jupyter Notebook
+.ipynb_checkpoints
+
+# IPython
+profile_default/
+ipython_config.py
+
+# pyenv
+# For a library or package, you might want to ignore these files since the code is
+# intended to run in multiple environments; otherwise, check them in:
+# .python-version
+
+# pipenv
+# According to pypa/pipenv#598, it is recommended to include Pipfile.lock in version control.
+# However, in case of collaboration, if having platform-specific dependencies or dependencies
+# having no cross-platform support, pipenv may install dependencies that don't work, or not
+# install all needed dependencies.
+#Pipfile.lock
+
+# PEP 582; used by e.g. github.com/David-OConnor/pyflow
+__pypackages__/
+
+# Celery stuff
+celerybeat-schedule
+celerybeat.pid
+
+# SageMath parsed files
+*.sage.py
+
+# Environments
+.env
+.venv
+env/
+venv/
+ENV/
+env.bak/
+venv.bak/
+
+# Spyder project settings
+.spyderproject
+.spyproject
+
+# Rope project settings
+.ropeproject
+
+# mkdocs documentation
+/site
+
+# mypy
+.mypy_cache/
+.dmypy.json
+dmypy.json
+
+# Pyre type checker
+.pyre/
+
+# pytype static type analyzer
+.pytype/
+
+# Cython debug symbols
+cython_debug/
+
+# static files generated from Django application using `collectstatic`
+media
+static
diff --git a/python/Makefile b/python/Makefile
new file mode 100644
index 0000000000000000000000000000000000000000..ccc13257fe572b456857657929c41d8199ad6529
--- /dev/null
+++ b/python/Makefile
@@ -0,0 +1,30 @@
+##
+# ##
+# corsika
+#
+# @file
+
+# find python3
+PYTHON=`which python3`
+
+# our testing targets
+.PHONY: tests flake black isort all
+
+all: mypy isort black flake tests
+
+tests:
+ ${PYTHON} -m pytest --cov=corsika tests
+
+flake:
+ ${PYTHON} -m flake8 corsika tests
+
+black:
+ ${PYTHON} -m black -t py37 corsika tests
+
+isort:
+ ${PYTHON} -m isort --atomic corsika tests
+
+mypy:
+ ${PYTHON} -m mypy corsika
+
+# end
diff --git a/python/README.md b/python/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..7e35dfe389dde4d372b1e38301c41c42a45c6fa8
--- /dev/null
+++ b/python/README.md
@@ -0,0 +1,8 @@
+# CORSIKA 8 - Python Library
+
+To install this into your global environment using `pip` (not recommended), run
+
+``` shell
+pip install --user git+https://gitlab.ikp.kit.edu/AirShowerPhysics/corsika/-/tree/master/python
+```
+
diff --git a/python/corsika/__init__.py b/python/corsika/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..078d688098095f06e2a14d2db4c203370993d578
--- /dev/null
+++ b/python/corsika/__init__.py
@@ -0,0 +1,17 @@
+"""
+ A Python interface to CORSIKA 8.
+
+ (c) Copyright 2018 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.
+"""
+
+from . import io
+from .io.library import Library
+
+# all imported objects
+__all__ = ["io", "Library"]
+
+__version__: str = "8.0.0-alpha"
diff --git a/python/corsika/io/__init__.py b/python/corsika/io/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..296572d6dee46a4b95a61b8ef18c549f120a1205
--- /dev/null
+++ b/python/corsika/io/__init__.py
@@ -0,0 +1,15 @@
+"""
+ The 'io' module provides for reading CORSIKA8 output files.
+
+ (c) Copyright 2018 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.
+"""
+
+from .hist import read_hist
+from .library import Library
+
+# all exported objects
+__all__ = ["read_hist", "Library"]
diff --git a/python/corsika/io/hist.py b/python/corsika/io/hist.py
new file mode 100644
index 0000000000000000000000000000000000000000..615db0faec0ec63cff8fe00fc5ebe1d8ce08f2ff
--- /dev/null
+++ b/python/corsika/io/hist.py
@@ -0,0 +1,70 @@
+"""
+ This file supports reading boost_histograms from CORSIKA8.
+
+ (c) Copyright 2018 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 boost_histogram as bh
+import numpy as np
+
+
+def read_hist(filename: str) -> bh.Histogram:
+ """
+ Read a histogram produced with CORSIKA8's `SaveBoostHistogram()` function.
+
+ Parameters
+ ----------
+ filename: str
+ The filename of the .npy file containing the histogram.
+
+ Returns
+ -------
+ hist: bh.Histogram
+ An initialized bh.Histogram instance.
+
+ Throws
+ ------
+ ValueError:
+ If the histogram type is not supported.
+ """
+
+ # load the filenames
+ d = np.load(filename)
+
+ # extract the axis and overflows information
+ axistypes = d["axistypes"].view("c")
+ overflow = d["overflow"]
+ underflow = d["underflow"]
+
+ # this is where we store the axes that we extract from the file.
+ axes = []
+
+ # we now loop over the axes
+ for i, (at, has_overflow, has_underflow) in enumerate(
+ zip(axistypes, overflow, underflow)
+ ):
+
+ # continuous histogram
+ if at == b"c":
+ axes.append(
+ bh.axis.Variable(
+ d[f"binedges_{i}"], overflow=has_overflow, underflow=has_underflow
+ )
+ )
+ # discrete histogram
+ elif at == b"d":
+ axes.append(bh.axis.IntCategory(d[f"bins_{i}"], growth=(not has_overflow)))
+
+ # and unknown histogram type
+ else:
+ raise ValueError(f"'{at}' is not a valid C8 histogram axistype.")
+
+ # create the histogram and fill it in
+ h = bh.Histogram(*axes)
+ h.view(flow=True)[:] = d["data"]
+
+ return h
diff --git a/python/corsika/io/library.py b/python/corsika/io/library.py
new file mode 100644
index 0000000000000000000000000000000000000000..236fc97ac972a11c0fb366876bd92cae20d6ab18
--- /dev/null
+++ b/python/corsika/io/library.py
@@ -0,0 +1,231 @@
+"""
+ This file allows for reading/working with C8 libraries.
+
+ (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.
+"""
+import logging
+import os
+import os.path as op
+import re
+from typing import Any, Dict, Optional, List
+
+import yaml
+
+from . import outputs
+
+
+class Library(object):
+ """
+ Represents a library ("run") of showers produced by C8.
+ """
+
+ def __init__(self, path: str):
+ """
+
+ Parameters
+ ----------
+ path: str
+ The path to the directory containing the library.
+
+
+ Raises
+ ------
+ ValueError
+ If `path` does not contain a valid CORSIKA8 library.
+ """
+
+ # check that this is a valid library
+ if not self.__valid_library(path):
+ raise ValueError(f"'{path}' does not contain a valid CORSIKA8 library.")
+
+ # store the top-level path
+ self.path = path
+
+ # load the config and summary files
+ self.config = self.load_config(path)
+ self.summary = self.load_summary(path)
+
+ # build the list of outputs
+ self.__outputs = self.__build_outputs(path)
+
+ @property
+ def names(self) -> List[str]:
+ """
+ Return the list of registered outputs.
+ """
+ return list(self.__outputs.keys())
+
+ @property
+ def modules(self) -> Dict[str, str]:
+ """
+ Return the list of registered outputs.
+ """
+ pass
+
+ def get(self, name: str) -> Optional[outputs.Output]:
+ """
+ Return the output with a given name.
+ """
+ if name in self.__outputs:
+ return self.__outputs[name]
+ else:
+ msg = f"Output with name '{name}' not available in this library."
+ logging.getLogger("corsika").warn(msg)
+ return None
+
+ @staticmethod
+ def load_config(path: str) -> Dict[str, Any]:
+ """
+ Load the top-level config from a given library path.
+
+
+ Parameters
+ ----------
+ path: str
+ The path to the directory containing the library.
+
+ Returns
+ -------
+ dict:
+ The config as a python dictionary.
+
+ Raises
+ ------
+ FileNotFoundError
+ If the config file cannot be found
+
+ """
+ with open(op.join(path, "config.yaml"), "r") as f:
+ return yaml.load(f, Loader=yaml.Loader)
+
+ @staticmethod
+ def load_summary(path: str) -> Dict[str, Any]:
+ """
+ Load the top-level summary from a given library path.
+
+
+ Parameters
+ ----------
+ path: str
+ The path to the directory containing the library.
+
+ Returns
+ -------
+ dict:
+ The summary as a python dictionary.
+
+ Raises
+ ------
+ FileNotFoundError
+ If the summary file cannot be found
+
+ """
+ with open(op.join(path, "summary.yaml"), "r") as f:
+ return yaml.load(f, Loader=yaml.Loader)
+
+ @staticmethod
+ def __valid_library(path: str) -> bool:
+ """
+ Check if the library pointed to by 'path' is a valid C8 library.
+
+ Parameters
+ ----------
+ path: str
+ The path to the directory containing the library.
+
+ Returns
+ -------
+ bool:
+ True if this is a valid C8 library.
+
+ """
+
+ # check that the config file exists
+ if not op.exists(op.join(path, "config.yaml")):
+ return False
+
+ # the config file exists, we load it
+ config = Library.load_config(path)
+
+ # and check that the config's "writer" key is correct
+ return config["creator"] == "CORSIKA8"
+
+ @staticmethod
+ def __build_outputs(path: str) -> Dict[str, outputs.Output]:
+ """
+ Build the outputs contained in this library.
+
+ This will print a warning message if a particular
+ output is invalid but will continue to load additional
+ outputs afterwards.
+
+ Parameters
+ ----------
+ path: str
+ The path to the directory containing this library.
+
+ Returns
+ -------
+ Dict[str, Output]:
+ A dictionary mapping names to initialized outputs.
+ """
+
+ # get a list of the subdirectories in the library
+ _, dirs, _ = next(os.walk(path))
+
+ # this is the dictionary where we store our components
+ components: Dict[str, Any] = {}
+
+ # loop over the subdirectories
+ for subdir in dirs:
+
+ # read the config file for this output
+ config = Library.load_config(op.join(path, subdir))
+
+ # the name keyword is our unique identifier
+ name = config.get("name")
+
+ # get the "type" keyword to identify this component
+ out_type = config.get("type")
+
+ # if `out_type` was None, this is an invalid output
+ if out_type is None or name is None:
+ msg = (
+ f"'{subdir}' does not contain a valid config."
+ "Missing 'type' or 'name' keyword."
+ )
+ logging.getLogger("corsika").warn(msg)
+ continue # skip to the next output, don't error
+
+ # we now have a valid component type, get the corresponding
+ # type from the proccesses subdirectory
+ try:
+
+ # instantiate the output and store it in our dict
+ component = getattr(outputs, out_type)(op.join(path, subdir))
+
+ # check if the read failed
+ if not component.is_good():
+ msg = (
+ f"'{name}' encountered an error while reading. "
+ "This process will be not be loaded."
+ )
+ logging.getLogger("corsika").warn(msg)
+ else:
+ components[name] = component
+
+ except AttributeError as e:
+ msg = (
+ f"Unable to instantiate an instance of '{out_type}' "
+ f"for a process called '{name}'"
+ )
+ logging.getLogger("corsika").warn(msg)
+ logging.getLogger("corsika").warn(e)
+ continue # skip to the next output, don't error
+
+ # and we are done building - return the constructed outputs
+ return components
diff --git a/python/corsika/io/outputs/__init__.py b/python/corsika/io/outputs/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..fa38acdde3c2dc8dd3cc351c16d27159cecd1efb
--- /dev/null
+++ b/python/corsika/io/outputs/__init__.py
@@ -0,0 +1,14 @@
+"""
+
+ (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.
+"""
+
+from .observation_plane import ObservationPlane
+from .track_writer import TrackWriter
+from .output import Output
+
+__all__ = ["Output", "ObservationPlane", "TrackWriter"]
diff --git a/python/corsika/io/outputs/observation_plane.py b/python/corsika/io/outputs/observation_plane.py
new file mode 100644
index 0000000000000000000000000000000000000000..271d145b4f076f1c7d747d6cc5e05b4a183d3655
--- /dev/null
+++ b/python/corsika/io/outputs/observation_plane.py
@@ -0,0 +1,87 @@
+"""
+ Read data written by ObservationPlane.
+
+ (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.
+"""
+import logging
+import os.path as op
+from typing import Any
+
+import pyarrow.parquet as pq
+
+from .output import Output
+
+
+class ObservationPlane(Output):
+ """
+ Read particle data from an ObservationPlane.
+ """
+
+ 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, "particles.parquet"))
+ except Exception as e:
+ logging.getLogger("corsika").warn(
+ f"An error occured loading an ObservationPlane: {e}"
+ )
+
+ 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 from this observation plane.
+
+ 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()
+ else:
+ raise ValueError(
+ (
+ f"Unknown format '{dtype}' for ObservationPlane. "
+ "We currently only support ['arrow', 'pandas']."
+ )
+ )
+
+ def __repr__(self) -> str:
+ """
+ Return a string representation of this class.
+ """
+ return f"ObservationPlane('{self.config['name']}')"
diff --git a/python/corsika/io/outputs/output.py b/python/corsika/io/outputs/output.py
new file mode 100644
index 0000000000000000000000000000000000000000..83a4c7928b786b9bad4301bec39b7e03aaf567f8
--- /dev/null
+++ b/python/corsika/io/outputs/output.py
@@ -0,0 +1,165 @@
+"""
+ This file defines the API for all output readers.
+
+ (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.
+"""
+import os.path as op
+from abc import ABC, abstractmethod
+from typing import Any, Dict
+
+import yaml
+
+
+class Output(ABC):
+ """
+ This class defines the abstract interface for all classes
+ that wish to provide reading support for CORSIKA8 outputs.
+ """
+
+ def __init__(self, path: str):
+ """
+ __init__ must load the output files and check
+ that it is valid.
+
+ Parameters
+ ----------
+ path: str
+ The path to the directory containing this output.
+ """
+ # load and store our path and config
+ self.path = path
+ self.__config = self.load_config(path)
+ self.__summary = self.load_summary(path)
+
+ @abstractmethod
+ 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.
+ """
+ pass
+
+ @abstractmethod
+ def astype(self, dtype: str, **kwargs: Any) -> Any:
+ """
+ Return the data for this output in the data format given by 'dtype'
+
+ Parameters
+ ----------
+ dtype: str
+ The data format to return the data in (i.e. numpy, pandas, etc.)
+ *args: Any
+ Additional arguments can be accepted by the output.
+ **kwargs: Any
+ Additional keyword arguments can be accepted by the output.
+
+ Returns
+ -------
+ Any:
+ The return type of this method is determined by `dtype`.
+ """
+ pass
+
+ @property
+ def config(self) -> Dict[str, Any]:
+ """
+ Return the config file for this output.
+
+ Parameters
+ ----------
+
+ Returns
+ -------
+ Dict[str, any]
+ The configuration file for this output.
+ """
+ return self.__config
+
+ @property
+ def summary(self) -> Dict[str, Any]:
+ """
+ Return the summary file for this output.
+
+ Parameters
+ ----------
+
+ Returns
+ -------
+ Dict[str, any]
+ The summary file for this output.
+ """
+ return self.__summary
+
+ @property
+ def data(self) -> Any:
+ """
+ Return the data in its default format.
+
+ We try to use Pandas as the default format for most data.
+
+ Parameters
+ ----------
+
+ Returns
+ -------
+ Any:
+ The data in its default format.
+ """
+ return self.astype()
+
+ @staticmethod
+ def load_config(path: str) -> Dict[str, Any]:
+ """
+ Load the top-level config from a given library path.
+
+ Parameters
+ ----------
+ path: str
+ The path to the directory containing the library.
+
+ Returns
+ -------
+ dict:
+ The config as a python dictionary.
+
+ Raises
+ ------
+ FileNotFoundError
+ If the config file cannot be found
+
+ """
+ with open(op.join(path, "config.yaml"), "r") as f:
+ return yaml.load(f, Loader=yaml.Loader)
+
+ @staticmethod
+ def load_summary(path: str) -> Dict[str, Any]:
+ """
+ Load the top-level summary from a given library path.
+
+ Parameters
+ ----------
+ path: str
+ The path to the directory containing the library.
+
+ Returns
+ -------
+ dict:
+ The summary as a python dictionary.
+
+ Raises
+ ------
+ FileNotFoundError
+ If the summary file cannot be found
+
+ """
+ with open(op.join(path, "summary.yaml"), "r") as f:
+ return yaml.load(f, Loader=yaml.Loader)
diff --git a/python/corsika/io/outputs/track_writer.py b/python/corsika/io/outputs/track_writer.py
new file mode 100644
index 0000000000000000000000000000000000000000..0859660d590182578e08ad4beaef5bb9a5c53498
--- /dev/null
+++ b/python/corsika/io/outputs/track_writer.py
@@ -0,0 +1,87 @@
+"""
+ Read data written by TrackWriter
+
+ (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.
+"""
+import logging
+import os.path as op
+from typing import Any
+
+import pyarrow.parquet as pq
+
+from .output import Output
+
+
+class TrackWriter(Output):
+ """
+ Read particle data from a TrackWriter
+ """
+
+ 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, "tracks.parquet"))
+ except Exception as e:
+ logging.getLogger("corsika").warn(
+ f"An error occured loading a TrackWriter: {e}"
+ )
+
+ 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 from this track writer.
+
+ 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()
+ else:
+ raise ValueError(
+ (
+ f"Unknown format '{dtype}' for TrackWriter. "
+ "We currently only support ['arrow', 'pandas']."
+ )
+ )
+
+ def __repr__(self) -> str:
+ """
+ Return a string representation of this class.
+ """
+ return f"TrackWriter('{self.config['name']}')"
diff --git a/python/setup.cfg b/python/setup.cfg
new file mode 100644
index 0000000000000000000000000000000000000000..f58d04eeb061f000b2d890aa26f542f70f3aa0c0
--- /dev/null
+++ b/python/setup.cfg
@@ -0,0 +1,69 @@
+[flake8]
+# use a slightly longer line to be consistent with black
+max-line-length = 88
+
+# E231 is missing whitespace after a comma
+# this contradicts the black formatting rules
+# and therefore creates spurious errors
+ignore = E231
+
+# we set various directories we want to exclude
+exclude =
+ # Don't bother checking in cache directories
+ __pycache__
+
+
+[isort]
+# use parenthesis for multi-line imports
+use_parentheses = true
+
+
+[mypy]
+# the primary Python version
+python_version = 3.7
+
+# allow returning Any
+# this creates excessive errors when using libraries
+# that don't have MyPy typing support
+warn_return_any = False
+
+# don't allow untyped functions
+disallow_untyped_defs = True
+
+# warn if any part of this config is mispelled
+warn_unused_configs = True
+
+# warn for missing type information
+warn_incomplete_stub = True
+
+# warn us if we don't return from a function explicitly
+warn_no_return = True
+
+# use incremental typing to speed things up
+incremental = True
+
+# show error contexts
+show_error_context = True
+
+# and show the column numbers for errors
+show_column_numbers = True
+
+# ignore missing types for setuptools
+[mypy-setuptools.*]
+ignore_missing_imports = True
+
+# ignore missing types for numpy
+[mypy-numpy.*]
+ignore_missing_imports = True
+
+# ignore missing types for matplotlib
+[mypy-matplotlib.*]
+ignore_missing_imports = True
+
+# ignore missing types for boost_histogram
+[mypy-boost_histogram.*]
+ignore_missing_imports = True
+
+# ignore missing types for pyarow
+[mypy-pyarrow.*]
+ignore_missing_imports = True
diff --git a/python/setup.py b/python/setup.py
new file mode 100644
index 0000000000000000000000000000000000000000..196c2ecbbccbc7e5caafee93be0af1c41746c94f
--- /dev/null
+++ b/python/setup.py
@@ -0,0 +1,51 @@
+from os import path
+from setuptools import setup
+
+# the stereo version
+__version__ = "8.0.0-alpha"
+
+# get the absolute path of this project
+here = path.abspath(path.dirname(__file__))
+
+# Get the long description from the README file
+with open(path.join(here, "README.md"), encoding="utf-8") as f:
+ long_description = f.read()
+
+# the standard setup info
+setup(
+ name="corsika",
+ version=__version__,
+ description="A Python package for working with CORSIKA 8.",
+ long_description=long_description,
+ long_description_content_type="text/markdown",
+ url="https://gitlab.ikp.kit.edu/AirShowerPhysics/corsika",
+ author="CORSIKA 8 Collaboration",
+ author_email="corsika-devel@lists.kit.edu",
+ classifiers=[
+ "Development Status :: 3 - Alpha",
+ "Intended Audience :: Science/Research",
+ "Topic :: Scientific/Engineering :: Physics",
+ "Programming Language :: Python :: 3.6",
+ "Programming Language :: Python :: 3.7",
+ "Programming Language :: Python :: 3.8",
+ ],
+ keywords=["cosmic ray", "physics", "air shower", "simulation"],
+ packages=["corsika"],
+ python_requires=">=3.6*, <4",
+ install_requires=["numpy", "pyyaml", "pyarrow", "boost_histogram"],
+ extras_require={
+ "test": [
+ "pytest",
+ "black",
+ "mypy",
+ "isort",
+ "coverage",
+ "pytest-cov",
+ "flake8",
+ ],
+ "pandas": ["pandas"],
+ },
+ scripts=[],
+ project_urls={"code": "https://gitlab.ikp.kit.edu/AirShowerPhysics/corsika"},
+ include_package_data=False,
+)
diff --git a/python/tests/__init__.py b/python/tests/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..c7d9b672b4c1bde7cc7725c619c0390a79e682c7
--- /dev/null
+++ b/python/tests/__init__.py
@@ -0,0 +1,56 @@
+"""
+ Tests for `corsika`
+
+ (c) Copyright 2018 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 os
+import os.path as op
+
+__all__ = ["build_directory"]
+
+
+def find_build_directory() -> str:
+ """
+ Return the absolute path to the CORSIKA8 build directory.
+
+ Returns
+ -------
+ str
+ The absolute path to the build directory.
+
+ Raises
+ ------
+ RuntimeError
+ If the build directory cannot be found or it is empty.
+ """
+
+ # check if we are running on Gitlab
+ gitlab_build_dir = os.getenv("CI_BUILDS_DIR")
+
+ # if we are running on Gitlab
+ if gitlab_build_dir is not None:
+ build_dir = op.abspath(gitlab_build_dir)
+ else: # otherwise, we are running locally
+ build_dir = op.abspath(
+ op.join(
+ op.dirname(__file__),
+ op.pardir,
+ op.pardir,
+ "build",
+ )
+ )
+
+ # check that the build directory contains 'CMakeCache.txt'
+ if not op.exists(op.join(build_dir, "CMakeCache.txt")):
+ raise RuntimeError("Python tests cannot find C8 build directory.")
+
+ return build_dir
+
+
+# find the build_directory once
+build_directory = find_build_directory()
diff --git a/python/tests/io/__init__.py b/python/tests/io/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..0d08be944644953fd08405a4502456e76c1203b9
--- /dev/null
+++ b/python/tests/io/__init__.py
@@ -0,0 +1,9 @@
+"""
+ Tests for `corsika.io`
+
+ (c) Copyright 2018 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.
+"""
diff --git a/python/tests/io/test_hist.py b/python/tests/io/test_hist.py
new file mode 100644
index 0000000000000000000000000000000000000000..7cf5b38fb7798c81680cf14d108e3b31dd3dcaf1
--- /dev/null
+++ b/python/tests/io/test_hist.py
@@ -0,0 +1,70 @@
+"""
+ Tests for `corsika.io.hist`
+
+ (c) Copyright 2018 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 os
+import os.path as op
+import subprocess
+
+import corsika
+
+from .. import build_directory
+
+# the directory containing 'testSaveBoostHistogram'
+bindir = op.join(build_directory, "Framework", "Utilities")
+
+
+def generate_hist() -> str:
+ """
+ Generate a test with `testSaveBoostHistogram`.
+
+ Returns
+ -------
+ str
+ The path to the generated histogram.
+ bool
+ If True, this file was regenerated.
+ """
+
+ # we construct the name of the bin
+ bin = op.join(bindir, "testSaveBoostHistogram")
+
+ # check if a histogram already exists
+ if op.exists(op.join(bin, "hist.npz")):
+ return op.join(bin, "hist.npz"), False
+ else:
+ # run the program - this generates "hist.npz" in the CWD
+ subprocess.call(bin)
+
+ return op.join(os.getcwd(), "hist.npz"), True
+
+
+def test_corsika_io() -> None:
+ """
+ Test I can corsika.io without a further import.
+ """
+ corsika.io.read_hist
+
+
+def test_corsika_read_hist() -> None:
+ """
+ Check that I can read in the test histograms with `read_hist`.
+ """
+
+ # generate a test histogram
+ filename, delete = generate_hist()
+
+ # and try and read in a histogram
+ h = corsika.io.read_hist(filename)
+
+ # and delete the generated histogram
+ if delete:
+ os.remove(filename)
+
+ # and check that it isn't empty
+ assert not h.empty()
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 3a12d7cc7401177052d3374be42c583073b3ae20..b8f1e709815fa0d41f9628482981cc4d879caec9 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -3,3 +3,4 @@ add_subdirectory (framework)
add_subdirectory (media)
add_subdirectory (stack)
add_subdirectory (modules)
+add_subdirectory (output)
diff --git a/tests/framework/testCascade.cpp b/tests/framework/testCascade.cpp
index 3e930911b095699bd09a10813c38c167c56b5dcc..44d277bae1bf03a7fdfe9fa12eef3ae2c2197f95 100644
--- a/tests/framework/testCascade.cpp
+++ b/tests/framework/testCascade.cpp
@@ -24,6 +24,8 @@
#include <corsika/media/HomogeneousMedium.hpp>
#include <corsika/media/NuclearComposition.hpp>
+#include <corsika/output/DummyOutputManager.hpp>
+
#include <SetupTestTrajectory.hpp>
#include <catch2/catch.hpp>
@@ -163,8 +165,10 @@ TEST_CASE("Cascade", "[Cascade]") {
Point(rootCS, {0_m, 0_m, 10_km}), 0_ns));
DummyTracking tracking;
- Cascade<DummyTracking, decltype(sequence), TestCascadeStack, TestCascadeStackView> EAS(
- env, tracking, sequence, stack);
+ DummyOutputManager output;
+ Cascade<DummyTracking, decltype(sequence), DummyOutputManager, TestCascadeStack,
+ TestCascadeStackView>
+ EAS(env, tracking, sequence, output, stack);
SECTION("full cascade") {
EAS.run();
diff --git a/tests/modules/testObservationPlane.cpp b/tests/modules/testObservationPlane.cpp
index 90a12288bd0b24c4c3021adaf66630fd1265d688..62983a2e08f49eef546e017b104cb104ead3de41 100644
--- a/tests/modules/testObservationPlane.cpp
+++ b/tests/modules/testObservationPlane.cpp
@@ -17,16 +17,17 @@
#include <corsika/framework/core/ParticleProperties.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/output/NoOutput.hpp>
+
#include <SetupTestEnvironment.hpp>
#include <SetupTestStack.hpp>
#include <SetupTestTrajectory.hpp>
using namespace corsika;
-TEST_CASE("ObservationPlane", "[proccesses][observation_plane]") {
+TEST_CASE("ObservationPlane", "interface") {
logging::set_level(logging::level::trace);
- corsika_logger->set_pattern("[%n:%^%-8l%$]: %v");
auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Oxygen);
auto const& cs = *csPtr;
@@ -54,23 +55,38 @@ TEST_CASE("ObservationPlane", "[proccesses][observation_plane]") {
SECTION("horizontal plane") {
Plane const obsPlane(Point(cs, {10_m, 0_m, 0_m}), DirectionVector(cs, {1., 0., 0.}));
- ObservationPlane obs(obsPlane, DirectionVector(cs, {0., 1., 0.}), "particles.dat",
- true);
+ ObservationPlane<NoOutput> obs(obsPlane, DirectionVector(cs, {0., 1., 0.}));
LengthType const length = obs.getMaxStepLength(particle, no_used_track);
ProcessReturn const ret = obs.doContinuous(particle, no_used_track, true);
CHECK(length / 10_m == Approx(1).margin(1e-4));
CHECK(ret == ProcessReturn::ParticleAbsorbed);
+
+ // particle does not reach plane:
+ {
+ setup::Trajectory no_hit_track =
+ setup::testing::make_track<setup::Trajectory>(line, 1_nm / constants::c);
+ LengthType const no_hit = obs.getMaxStepLength(particle, no_hit_track);
+ CHECK(no_hit == std::numeric_limits<double>::infinity() * 1_m);
+ }
+
+ // particle past plane:
+ {
+ particle.setPosition({cs, {0_m, 0_m, -1_m}});
+ setup::Trajectory no_hit_track =
+ setup::testing::make_track<setup::Trajectory>(line, 1_nm / constants::c);
+ LengthType const no_hit = obs.getMaxStepLength(particle, no_hit_track);
+ CHECK(no_hit == std::numeric_limits<double>::infinity() * 1_m);
+ }
}
SECTION("transparent plane") {
Plane const obsPlane(Point(cs, {1_m, 0_m, 0_m}), DirectionVector(cs, {1., 0., 0.}));
- ObservationPlane obs(obsPlane, DirectionVector(cs, {0., 0., 1.}), "particles.dat",
- false);
+ ObservationPlane<NoOutput> obs(obsPlane, DirectionVector(cs, {0., 0., 1.}), false);
LengthType const length = obs.getMaxStepLength(particle, no_used_track);
- ProcessReturn const ret = obs.doContinuous(particle, no_used_track, true);
+ ProcessReturn const ret = obs.doContinuous(particle, no_used_track, false);
CHECK(length / 1_m == Approx(1).margin(1e-4));
CHECK(ret == ProcessReturn::Ok);
@@ -85,8 +101,7 @@ TEST_CASE("ObservationPlane", "[proccesses][observation_plane]") {
Plane const obsPlane(Point(cs, {10_m, 5_m, 5_m}),
DirectionVector(cs, {1, 0.1, -0.05}));
- ObservationPlane obs(obsPlane, DirectionVector(cs, {0., 1., 0.}), "particles.dat",
- true);
+ ObservationPlane<NoOutput> obs(obsPlane, DirectionVector(cs, {0., 1., 0.}));
LengthType const length = obs.getMaxStepLength(particle, no_used_track);
ProcessReturn const ret = obs.doContinuous(particle, no_used_track, true);
@@ -94,4 +109,11 @@ TEST_CASE("ObservationPlane", "[proccesses][observation_plane]") {
CHECK(length / 10_m == Approx(1.1375).margin(1e-4));
CHECK(ret == ProcessReturn::ParticleAbsorbed);
}
+
+ SECTION("output") {
+ Plane const obsPlane(Point(cs, {1_m, 0_m, 0_m}), DirectionVector(cs, {1., 0., 0.}));
+ ObservationPlane<NoOutput> obs(obsPlane, DirectionVector(cs, {0., 0., 1.}), false);
+ auto const cfg = obs.getConfig();
+ CHECK(cfg["type"]);
+ }
}
diff --git a/tests/output/CMakeLists.txt b/tests/output/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..98a43b4c218bf8d9e21bb7b9e56ac1275cbc9b0e
--- /dev/null
+++ b/tests/output/CMakeLists.txt
@@ -0,0 +1,16 @@
+set (test_output_sources
+ TestMain.cpp
+ testOutputManager.cpp
+ testDummyOutputManager.cpp
+ testParquetStreamer.cpp
+ testWriterObservationPlane.cpp
+ testWriterTracks.cpp
+ )
+
+CORSIKA_ADD_TEST (testOutput SOURCES ${test_output_sources})
+
+target_compile_definitions (
+ testOutput
+ PRIVATE
+ REFDATADIR="${CMAKE_CURRENT_SOURCE_DIR}"
+)
diff --git a/tests/output/TestMain.cpp b/tests/output/TestMain.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..51532584b8e03b35d79301543ac8f80b598ba544
--- /dev/null
+++ b/tests/output/TestMain.cpp
@@ -0,0 +1,11 @@
+/*
+ * (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.
+ */
+
+#define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one
+ // cpp file
+#include <catch2/catch.hpp>
diff --git a/tests/output/testDummyOutputManager.cpp b/tests/output/testDummyOutputManager.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..cc571832f601c1e9814ca25e1de9d4db152d1afc
--- /dev/null
+++ b/tests/output/testDummyOutputManager.cpp
@@ -0,0 +1,32 @@
+/*
+ * (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.
+ */
+
+#include <catch2/catch.hpp>
+
+#include <corsika/output/DummyOutputManager.hpp>
+#include <corsika/framework/core/Logging.hpp>
+
+using namespace corsika;
+
+class DummyOutput {};
+
+TEST_CASE("DummyOutputManager", "interface") {
+
+ logging::set_level(logging::level::info);
+
+ // output manager performs nothing, no action, just interface
+ DummyOutputManager output;
+
+ DummyOutput test;
+ output.add("test", test);
+
+ output.startOfLibrary();
+ output.startOfShower();
+ output.endOfShower();
+ output.endOfLibrary();
+}
diff --git a/tests/output/testOutputManager.cpp b/tests/output/testOutputManager.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..58f90589c2206ed1e01607cd216a6b76c53a20c9
--- /dev/null
+++ b/tests/output/testOutputManager.cpp
@@ -0,0 +1,184 @@
+/*
+ * (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.
+ */
+
+#include <catch2/catch.hpp>
+
+#include <boost/filesystem.hpp>
+
+#include <corsika/framework/core/Logging.hpp>
+
+#include <corsika/output/OutputManager.hpp>
+#include <corsika/output/NoOutput.hpp>
+
+using namespace corsika;
+
+struct DummyNoOutput : public NoOutput {
+ void check() {
+ NoOutput::startOfLibrary("./");
+ NoOutput::startOfShower();
+ NoOutput::endOfShower();
+ NoOutput::endOfLibrary();
+ NoOutput::getConfig();
+ NoOutput::getSummary();
+ }
+ void checkWrite() { NoOutput::write(Code::Unknown, 1_eV, 1_m, 1_m); }
+};
+
+struct DummyOutput : public BaseOutput {
+
+ mutable bool isConfig_ = false;
+ mutable bool isSummary_ = false;
+ bool startLibrary_ = false;
+ bool startShower_ = false;
+ bool endLibrary_ = false;
+ bool endShower_ = false;
+
+ void startOfLibrary(boost::filesystem::path const&) { startLibrary_ = true; }
+
+ void startOfShower() {
+ BaseOutput::startOfShower();
+ startShower_ = true;
+ }
+
+ void endOfShower() { endShower_ = true; }
+
+ void endOfLibrary() { endLibrary_ = true; }
+
+ YAML::Node getConfig() const {
+ isConfig_ = true;
+ return YAML::Node();
+ }
+
+ YAML::Node getSummary() {
+ isSummary_ = true;
+ return BaseOutput::getSummary();
+ }
+};
+
+TEST_CASE("OutputManager") {
+
+ logging::set_level(logging::level::info);
+
+ SECTION("standard") {
+
+ // preparation
+ if (boost::filesystem::exists("./out_test")) {
+ boost::filesystem::remove_all("./out_test");
+ }
+
+ // output manager performs nothing, no action, just interface
+ OutputManager output("check", "./out_test");
+
+ CHECK(boost::filesystem::is_directory("./out_test/check"));
+
+ DummyOutput test;
+ output.add("test", test);
+
+ CHECK_THROWS(output.add(
+ "test",
+ test)); // should emit warning which cannot be catched, but no action or failure
+
+ CHECK(test.isConfig_);
+ test.isConfig_ = false;
+
+ output.startOfLibrary();
+ CHECK(test.startLibrary_);
+ test.startLibrary_ = false;
+
+ output.startOfShower();
+ CHECK(test.startShower_);
+ test.startShower_ = false;
+
+ output.endOfShower();
+ CHECK(test.endShower_);
+ test.endShower_ = false;
+
+ output.endOfLibrary();
+ CHECK(test.endLibrary_);
+ CHECK(test.isSummary_);
+ test.isSummary_ = false;
+ test.endLibrary_ = false;
+
+ CHECK(boost::filesystem::exists("./out_test/check/test/summary.yaml"));
+ }
+
+ SECTION("auto-write") {
+
+ // preparation
+ if (boost::filesystem::exists("./out_test")) {
+ boost::filesystem::remove_all("./out_test");
+ }
+
+ // output manager performs nothing, no action, just interface
+ OutputManager* output = new OutputManager("check", "./out_test");
+
+ CHECK(boost::filesystem::is_directory("./out_test/check"));
+
+ DummyOutput test;
+ output->add("test", test);
+ output->startOfLibrary();
+ output->startOfShower();
+
+ // check support for closing automatically
+ delete output;
+ output = 0;
+
+ CHECK(boost::filesystem::exists("./out_test/check/test/summary.yaml"));
+ }
+
+ SECTION("failures") {
+
+ logging::set_level(logging::level::info);
+
+ // preparation
+ if (boost::filesystem::exists("./out_test")) {
+ boost::filesystem::remove_all("./out_test");
+ }
+
+ // output manager performs nothing, no action, just interface
+ OutputManager output("check", "./out_test");
+ CHECK_THROWS(new OutputManager("check", "./out_test"));
+
+ // CHECK_THROWS(output.startOfShower());
+ // CHECK_THROWS(output.endOfShower());
+ CHECK_THROWS(output.endOfLibrary());
+
+ output.startOfLibrary();
+
+ CHECK_THROWS(output.startOfLibrary());
+ // CHECK_THROWS(output.endOfShower());
+ // CHECK_THROWS(output.endOfLibrary());
+
+ output.startOfShower();
+
+ CHECK_THROWS(output.startOfLibrary());
+ // CHECK_THROWS(output.startOfShower());
+ // CHECK_THROWS(output.endOfLibrary());
+
+ output.endOfShower();
+
+ CHECK_THROWS(output.startOfLibrary());
+ // CHECK_THROWS(output.startOfShower());
+ // CHECK_THROWS(output.endOfShower());
+
+ output.endOfLibrary();
+
+ // CHECK_THROWS(output.endOfShower());
+ // CHECK_THROWS(output.startOfShower());
+ // CHECK_THROWS(output.endOfLibrary());
+ }
+
+ SECTION("NoOutput") {
+ // this is one of the classes where testing is a bit useless, but we can at least make
+ // sure the interface exists.
+ DummyNoOutput nothing;
+
+ nothing.check();
+ nothing.checkWrite();
+ }
+}
\ No newline at end of file
diff --git a/tests/output/testParquetStreamer.cpp b/tests/output/testParquetStreamer.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..387e5cade02d9e9cf1f1c1532cc1594bffb65d4a
--- /dev/null
+++ b/tests/output/testParquetStreamer.cpp
@@ -0,0 +1,56 @@
+/*
+ * (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.
+ */
+
+#include <catch2/catch.hpp>
+
+#include <boost/filesystem.hpp>
+
+#include <corsika/output/ParquetStreamer.hpp>
+#include <corsika/framework/core/Logging.hpp>
+
+using namespace corsika;
+
+TEST_CASE("ParquetStreamer") {
+
+ logging::set_level(logging::level::info);
+
+ SECTION("standard") {
+
+ // preparation
+ if (boost::filesystem::exists("./parquet_test.parquet")) {
+ boost::filesystem::remove_all("./parquet_test.parquet");
+ }
+
+ ParquetStreamer test;
+ CHECK_FALSE(test.isInit());
+ CHECK_THROWS(test.getWriter());
+ test.initStreamer("./parquet_test.parquet");
+
+ test.addField("testint", parquet::Repetition::REQUIRED, parquet::Type::INT32,
+ parquet::ConvertedType::INT_32);
+ test.addField("testfloat", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
+ parquet::ConvertedType::NONE);
+
+ test.enableCompression(1);
+
+ test.buildStreamer();
+ CHECK(test.isInit());
+
+ int testint = 1;
+ double testfloat = 2.0;
+
+ std::shared_ptr<parquet::StreamWriter> writer = test.getWriter();
+ (*writer) << static_cast<int>(testint) << static_cast<int>(testint)
+ << static_cast<float>(testfloat) << parquet::EndRow;
+
+ test.closeStreamer();
+ CHECK_THROWS(test.getWriter());
+ CHECK_FALSE(test.isInit());
+ CHECK(boost::filesystem::exists("./parquet_test.parquet"));
+ }
+}
\ No newline at end of file
diff --git a/tests/output/testWriterObservationPlane.cpp b/tests/output/testWriterObservationPlane.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..4b780a2b4bbcde53d531f4a6a113d4b2a7087f3a
--- /dev/null
+++ b/tests/output/testWriterObservationPlane.cpp
@@ -0,0 +1,50 @@
+/*
+ * (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.
+ */
+
+#include <catch2/catch.hpp>
+
+#include <boost/filesystem.hpp>
+
+#include <corsika/modules/writers/ObservationPlaneWriterParquet.hpp>
+
+#include <corsika/framework/core/Logging.hpp>
+#include <corsika/framework/geometry/QuantityVector.hpp>
+
+using namespace corsika;
+
+struct TestWriterPlane : public ObservationPlaneWriterParquet {
+
+ YAML::Node getConfig() const { return YAML::Node(); }
+
+ void checkWrite() {
+ ObservationPlaneWriterParquet::write(Code::Unknown, 1_eV, 2_m, 3_m);
+ }
+};
+
+TEST_CASE("ObservationPlaneWriterParquet") {
+
+ logging::set_level(logging::level::info);
+
+ SECTION("standard") {
+
+ // preparation
+ if (boost::filesystem::exists("./output_dir")) {
+ boost::filesystem::remove_all("./output_dir");
+ }
+ boost::filesystem::create_directory("./output_dir");
+
+ TestWriterPlane test;
+ test.startOfLibrary("./output_dir");
+ test.startOfShower();
+ test.checkWrite();
+ test.endOfShower();
+ test.endOfLibrary();
+
+ CHECK(boost::filesystem::exists("./output_dir/particles.parquet"));
+ }
+}
diff --git a/tests/output/testWriterTracks.cpp b/tests/output/testWriterTracks.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..38cf0cf5985a956aaef03b589b93a1b6b9f9b260
--- /dev/null
+++ b/tests/output/testWriterTracks.cpp
@@ -0,0 +1,49 @@
+/*
+ * (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.
+ */
+
+#include <catch2/catch.hpp>
+
+#include <boost/filesystem.hpp>
+
+#include <corsika/modules/writers/TrackWriterParquet.hpp>
+
+#include <corsika/framework/core/Logging.hpp>
+
+using namespace corsika;
+
+struct TestWriterTrack : public TrackWriterParquet {
+
+ YAML::Node getConfig() const { return YAML::Node(); }
+
+ void checkWrite() {
+ TrackWriterParquet::write(Code::Unknown, 1_eV, {2_m, 3_m, 4_m}, {5_m, 6_m, 7_m});
+ }
+};
+
+TEST_CASE("TrackWriterParquet") {
+
+ logging::set_level(logging::level::info);
+
+ SECTION("standard") {
+
+ // preparation
+ if (boost::filesystem::exists("./output_dir_tracks")) {
+ boost::filesystem::remove_all("./output_dir_tracks");
+ }
+ boost::filesystem::create_directory("./output_dir_tracks");
+
+ TestWriterTrack test;
+ test.startOfLibrary("./output_dir_tracks");
+ test.startOfShower();
+ test.checkWrite();
+ test.endOfShower();
+ test.endOfLibrary();
+
+ CHECK(boost::filesystem::exists("./output_dir_tracks/tracks.parquet"));
+ }
+}