diff --git a/corsika/detail/modules/ParticleCut.inl b/corsika/detail/modules/ParticleCut.inl
index b5853df3bd27b00cf8af56a37b82dca0d82613e9..789cf5b5972a1a00ccbb6b02117afad12bb05c33 100644
--- a/corsika/detail/modules/ParticleCut.inl
+++ b/corsika/detail/modules/ParticleCut.inl
@@ -13,11 +13,13 @@
namespace corsika {
template <typename TOutput>
- inline ParticleCut<TOutput>::ParticleCut(TOutput& output, HEPEnergyType const eEleCut,
+ template <typename... TArgs>
+ inline ParticleCut<TOutput>::ParticleCut(HEPEnergyType const eEleCut,
HEPEnergyType const ePhoCut,
HEPEnergyType const eHadCut,
- HEPEnergyType const eMuCut, bool const inv)
- : output_(output)
+ HEPEnergyType const eMuCut, bool const inv,
+ TArgs&&... args)
+ : TOutput(std::forward<TArgs>(args)...)
, doCutInv_(inv)
, energy_cut_(0_GeV)
, energy_timecut_(0_GeV)
@@ -45,10 +47,8 @@ namespace corsika {
template <typename TOutput>
inline ParticleCut<TOutput>::ParticleCut(
- TOutput& output, std::unordered_map<Code const, HEPEnergyType const> const& eCuts,
- bool const inv)
- : output_(output)
- , doCutInv_(inv)
+ std::unordered_map<Code const, HEPEnergyType const> const& eCuts, bool const inv)
+ : doCutInv_(inv)
, energy_cut_(0_GeV)
, energy_timecut_(0_GeV)
, energy_invcut_(0_GeV)
@@ -127,8 +127,8 @@ namespace corsika {
auto particle = vS.begin();
while (particle != vS.end()) {
if (checkCutParticle(particle)) {
- output_.write(particle.getPosition(), particle.getPID(),
- particle.getKineticEnergy());
+ this->write(particle.getPosition(), particle.getPID(),
+ particle.getKineticEnergy());
particle.erase();
}
++particle; // next entry in SecondaryView
@@ -141,8 +141,7 @@ namespace corsika {
inline ProcessReturn ParticleCut<TOutput>::doContinuous(TParticle& particle, TTrajectory const&,
bool const) {
if (checkCutParticle(particle)) {
- output_.write(particle.getPosition(), particle.getPID(),
- particle.getKineticEnergy());
+ this->write(particle.getPosition(), particle.getPID(), particle.getKineticEnergy());
CORSIKA_LOG_TRACE("removing during continuous");
// signal to upstream code that this particle was deleted
return ProcessReturn::ParticleAbsorbed;
@@ -181,4 +180,18 @@ namespace corsika {
energy_timecut_ = 0_GeV;
}
+ template <typename TOutput>
+ inline YAML::Node ParticleCut<TOutput>::getConfig() const {
+
+ YAML::Node node;
+ node["type"] = "ParticleCut";
+ node["units"]["energy"] = "GeV";
+ node["energy_invcut"] = energy_invcut_ / 1_GeV;
+ node["inv_count"] = inv_count_;
+ node["energy_cut"] = energy_cut_ / 1_GeV;
+ node["energy_timecut_"] = energy_timecut_ / 1_GeV;
+
+ return node;
+ }
+
} // namespace corsika
diff --git a/corsika/detail/modules/writers/ParticleCutWriterParquet.inl b/corsika/detail/modules/writers/ParticleCutWriterParquet.inl
new file mode 100644
index 0000000000000000000000000000000000000000..1405e9416a8ff916b6775b5b18ee1b3e4441ffc5
--- /dev/null
+++ b/corsika/detail/modules/writers/ParticleCutWriterParquet.inl
@@ -0,0 +1,64 @@
+/*
+ * (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 ParticleCutWriterParquet::ParticleCutWriterParquet()
+ : output_() {}
+
+ inline void ParticleCutWriterParquet::startOfLibrary(
+ boost::filesystem::path const& directory) {
+
+ // setup the streamer
+ output_.initStreamer((directory / "energyloss.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("x", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
+ parquet::ConvertedType::NONE);
+ output_.addField("y", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
+ parquet::ConvertedType::NONE);
+ output_.addField("z", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
+ parquet::ConvertedType::NONE);
+
+ // and build the streamer
+ output_.buildStreamer();
+ }
+
+ inline void ParticleCutWriterParquet::startOfShower(unsigned int const showerId) {
+ showerId_ = showerId;
+ }
+
+ inline void ParticleCutWriterParquet::endOfShower(unsigned int const) {}
+
+ inline void ParticleCutWriterParquet::endOfLibrary() { output_.closeStreamer(); }
+
+ inline void ParticleCutWriterParquet::write(Point const& point, Code const pid,
+ HEPEnergyType energy) {
+
+ auto location{point.getCoordinates()};
+
+ // write the next row - we must write `shower_` first.
+ // clang-format off
+ *(output_.getWriter())
+ << showerId_
+ << static_cast<int>(get_PDG(pid))
+ << static_cast<float>(energy / 1_GeV)
+ << static_cast<float>(location[0] / 1_m)
+ << static_cast<float>(location[1] / 1_m)
+ << static_cast<float>(location[2] / 1_m)
+ << parquet::EndRow;
+ // clang-format on
+ }
+
+} // namespace corsika
diff --git a/corsika/modules/ParticleCut.hpp b/corsika/modules/ParticleCut.hpp
index 3db5d0d30ec7819d3800de359e70b892d8c8c7fc..7ce7139ed4ee6c7cd013a299a5b98b992b235647 100644
--- a/corsika/modules/ParticleCut.hpp
+++ b/corsika/modules/ParticleCut.hpp
@@ -15,6 +15,8 @@
#include <corsika/framework/process/SecondariesProcess.hpp>
#include <corsika/framework/process/ContinuousProcess.hpp>
+#include <corsika/modules/writers/WriterOff.hpp>
+
namespace corsika {
/**
simple ParticleCut process. Goes through the secondaries of an interaction and
@@ -24,9 +26,10 @@ namespace corsika {
for each particle. Special constructors for cuts by the following groups are
implemented: (electrons,positrons), photons, hadrons and muons.
**/
- template <typename TOutput = EnergyLossWriterOff>
+ template <typename TOutput = WriterOff>
class ParticleCut : public SecondariesProcess<ParticleCut<TOutput>>,
- public ContinuousProcess<ParticleCut<TOutput>> {
+ public ContinuousProcess<ParticleCut<TOutput>>,
+ public TOutput {
public:
/**
@@ -34,30 +37,17 @@ namespace corsika {
* hadrons (including nuclei with energy per nucleon) and muons
* invisible particles (neutrinos) can be cut or not
*/
- ParticleCut(TOutput& output, HEPEnergyType const eEleCut, HEPEnergyType const ePhoCut,
- HEPEnergyType const eHadCut, HEPEnergyType const eMuCut, bool const inv);
-
- /**
- * Same, but with no output.
- */
+ template <typename... TArgs>
ParticleCut(HEPEnergyType const eEleCut, HEPEnergyType const ePhoCut,
- HEPEnergyType const eHadCut, HEPEnergyType const eMuCut, bool const inv)
- : ParticleCut(*(new EnergyLossWriterOff()), eEleCut, ePhoCut, eHadCut, eMuCut,
- inv) {}
+ HEPEnergyType const eHadCut, HEPEnergyType const eMuCut, bool const inv,
+ TArgs&&... args);
/**
* Threshold for specific particles redefined. EM and invisible particles can be set
* to be discarded altogether.
*/
- ParticleCut(TOutput& output,
- std::unordered_map<Code const, HEPEnergyType const> const& eCuts,
- bool const inv);
- /**
- * Same, but with no output.
- */
ParticleCut(std::unordered_map<Code const, HEPEnergyType const> const& eCuts,
- bool const inv)
- : ParticleCut(*(new EnergyLossWriterOff()), eCuts, inv) {}
+ bool const inv);
template <typename TStackView>
void doSecondaries(TStackView&);
@@ -97,6 +87,9 @@ namespace corsika {
//! returns number of invisible particles
unsigned int getNumberInvParticles() const { return inv_count_; }
+ // get configuration of this node
+ YAML::Node getConfig() const override;
+
private:
template <typename TParticle>
bool checkCutParticle(TParticle const& p);
@@ -108,7 +101,6 @@ namespace corsika {
bool isInvisible(Code const&) const;
private:
- TOutput& output_;
bool doCutInv_;
HEPEnergyType energy_cut_ = 0 * electronvolt;
HEPEnergyType energy_timecut_ = 0 * electronvolt;
diff --git a/corsika/modules/writers/ParticleCutWriterParquet.hpp b/corsika/modules/writers/ParticleCutWriterParquet.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..466bf55719e044a2749b26a39a8388e0d43cfbf3
--- /dev/null
+++ b/corsika/modules/writers/ParticleCutWriterParquet.hpp
@@ -0,0 +1,61 @@
+/*
+ * (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 ParticleCutWriterParquet : public BaseOutput {
+
+ public:
+ /**
+ * Construct a new writer.
+ */
+ ParticleCutWriterParquet();
+
+ /**
+ * Called at the start of each library.
+ */
+ void startOfLibrary(boost::filesystem::path const& directory) final override;
+
+ /**
+ * Called at the start of each shower.
+ */
+ void startOfShower(unsigned int const showerId) final override;
+
+ /**
+ * Called at the end of each shower.
+ */
+ void endOfShower(unsigned int const showerId) 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;
+
+ /**
+ * Write a particle cut.
+ */
+ void write(Point const& point, Code const PID, HEPEnergyType const energy);
+
+ private:
+ ParquetStreamer output_; ///< The parquet streamer for this process.
+ unsigned int showerId_{0};
+ }; // namespace corsika
+
+} // namespace corsika
+
+#include <corsika/detail/modules/writers/ParticleCutWriterParquet.inl>
diff --git a/examples/vertical_EAS.cpp b/examples/vertical_EAS.cpp
index 8f9f984dc5da4d35b761cd92187e885149b22a02..2023d527bd48ffda3af586c6d1cc3721b3b82ced 100644
--- a/examples/vertical_EAS.cpp
+++ b/examples/vertical_EAS.cpp
@@ -51,6 +51,7 @@
#include <corsika/modules/StackInspector.hpp>
#include <corsika/modules/TrackWriter.hpp>
#include <corsika/modules/ParticleCut.hpp>
+#include <corsika/modules/writers/ParticleCutWriterParquet.hpp>
#include <corsika/modules/Pythia8.hpp>
#include <corsika/modules/Sibyll.hpp>
#include <corsika/modules/UrQMD.hpp>
@@ -245,7 +246,8 @@ int main(int argc, char** argv) {
decaySibyll.printDecayConfig();
- ParticleCut cut{60_GeV, 60_GeV, 60_GeV, 60_GeV, true};
+ ParticleCut<ParticleCutWriterParquet> cut{60_GeV, 60_GeV, 60_GeV, 60_GeV, true};
+ output.add("cuts", cut);
corsika::urqmd::UrQMD urqmd;
InteractionCounter urqmdCounted{urqmd};
diff --git a/python/corsika/io/outputs/__init__.py b/python/corsika/io/outputs/__init__.py
index b86fb6e70e83d06d4570ba65dc2c9a627177d5e7..d39d9a8b4bd68465ab1d799f10c15bb20032cbc9 100644
--- a/python/corsika/io/outputs/__init__.py
+++ b/python/corsika/io/outputs/__init__.py
@@ -11,6 +11,14 @@ from .observation_plane import ObservationPlane
from .track_writer import TrackWriter
from .longitudinal_profile import LongitudinalProfile
from .bethe_bloch import BetheBlochPDG
+from .particle_cut import ParticleCut
from .output import Output
-__all__ = ["Output", "ObservationPlane", "TrackWriter", "LongitudinalProfile"]
+__all__ = [
+ "Output",
+ "ObservationPlane",
+ "TrackWriter",
+ "LongitudinalProfile",
+ "BetheBlochPDG",
+ "ParticleCut",
+]
diff --git a/python/corsika/io/outputs/particle_cut.py b/python/corsika/io/outputs/particle_cut.py
new file mode 100644
index 0000000000000000000000000000000000000000..a55b5fe718cf788fdba051dfa891e77b2d75ede7
--- /dev/null
+++ b/python/corsika/io/outputs/particle_cut.py
@@ -0,0 +1,87 @@
+"""
+ Read data written by ParticleCut.
+
+ (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 ParticleCut(Output):
+ """
+ Read particle data from an ParticleCut.
+ """
+
+ 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, "energyloss.parquet"))
+ except Exception as e:
+ logging.getLogger("corsika").warn(
+ f"An error occured loading a ParticleCut: {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 particle cut instance.
+
+ 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 ParticleCut. "
+ "We currently only support ['arrow', 'pandas']."
+ )
+ )
+
+ def __repr__(self) -> str:
+ """
+ Return a string representation of this class.
+ """
+ return f"ParticleCut('{self.config['name']}')"