diff --git a/corsika/detail/modules/writers/EnergyLossWriter.inl b/corsika/detail/modules/writers/EnergyLossWriter.inl
new file mode 100644
index 0000000000000000000000000000000000000000..1fbf3b5ba2c961bb4f61ea39f1647362d5c83d3e
--- /dev/null
+++ b/corsika/detail/modules/writers/EnergyLossWriter.inl
@@ -0,0 +1,218 @@
+/*
+ * (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>
+#include <corsika/framework/utility/FindXmax.hpp>
+
+#include <corsika/media/ShowerAxis.hpp>
+
+#include <exception>
+
+namespace corsika {
+
+ template <typename TOutput>
+ inline EnergyLossWriter<TOutput>::EnergyLossWriter(ShowerAxis const& axis,
+ GrammageType dX,
+ unsigned int const nBins,
+ GrammageType dX_threshold)
+ : showerAxis_(axis)
+ , dX_(dX)
+ , nBins_(nBins)
+ , dX_threshold_(dX_threshold) {}
+
+ template <typename TOutput>
+ inline void EnergyLossWriter<TOutput>::startOfLibrary(
+ boost::filesystem::path const& directory) {
+ TOutput::startOfLibrary(directory);
+ }
+
+ template <typename TOutput>
+ inline void EnergyLossWriter<TOutput>::startOfShower(unsigned int const showerId) {
+ TOutput::startOfShower(showerId);
+ }
+
+ template <typename TOutput>
+ inline void EnergyLossWriter<TOutput>::endOfLibrary() {
+ TOutput::endOfLibrary();
+ }
+
+ template <typename TOutput>
+ inline void EnergyLossWriter<TOutput>::endOfShower(unsigned int const showerId) {
+
+ int iRow{0};
+ for (Profile const& row : profile_) {
+ TOutput::write(showerId, iRow * dX_, row.at(static_cast<int>(ProfileIndex::Total)));
+ iRow++;
+ }
+
+ TOutput::endOfShower(showerId);
+
+ // reset profile
+ profile_.clear();
+ profile_.resize(nBins_);
+ }
+
+ template <typename TOutput>
+ template <typename TTrack>
+ inline void EnergyLossWriter<TOutput>::write(TTrack const& track, Code const PID,
+ HEPEnergyType const dE) {
+
+ GrammageType grammageStart = showerAxis_.getProjectedX(track.getPosition(0));
+ GrammageType grammageEnd = showerAxis_.getProjectedX(track.getPosition(1));
+
+ if (grammageStart > grammageEnd) { // particle going upstream
+ std::swap(grammageStart, grammageEnd);
+ }
+
+ GrammageType const deltaX = grammageEnd - grammageStart;
+
+ if (deltaX < dX_threshold_) {
+ this->write(track.getPosition(0), PID, dE);
+ return;
+ }
+
+ // only register the range that is covered by the profile
+ int const maxBin = int(profile_.size() - 1);
+ int binStart = grammageStart / dX_;
+ if (binStart < 0) binStart = 0;
+ if (binStart > maxBin) binStart = maxBin;
+ int binEnd = grammageEnd / dX_;
+ if (binEnd < 0) binEnd = 0;
+ if (binEnd > maxBin) binEnd = maxBin;
+
+ CORSIKA_LOGGER_TRACE(
+ TOutput::getLogger(), "energy deposit of dE={} GeV between {} and {}", dE / 1_GeV,
+ grammageStart / 1_g * square(1_cm), grammageEnd / 1_g * square(1_cm));
+
+ auto energyCount = HEPEnergyType::zero();
+
+ auto const factor = dE / deltaX;
+ auto fill = [&](int const bin, GrammageType const weight) {
+ auto const increment = factor * weight;
+ profile_[bin][static_cast<int>(ProfileIndex::Total)] += increment;
+ energyCount += increment;
+
+ CORSIKA_LOGGER_TRACE(TOutput::getLogger(),
+ "filling bin={} with weight {} : dE={} GeV ", bin, weight,
+ increment / 1_GeV);
+ };
+
+ // fill longitudinal profile
+ if (binStart == binEnd) {
+ fill(binStart, deltaX);
+ } else {
+ fill(binStart, ((1 + binStart) * dX_ - grammageStart));
+ fill(binEnd, (grammageEnd - binEnd * dX_));
+ for (int bin = binStart + 1; bin < binEnd; ++bin) { fill(bin, dX_); }
+ }
+
+ CORSIKA_LOGGER_TRACE(TOutput::getLogger(), "total energy added to histogram: {} GeV ",
+ energyCount / 1_GeV);
+ }
+
+ template <typename TOutput>
+ inline void EnergyLossWriter<TOutput>::write(Point const& point, Code const PID,
+ HEPEnergyType const dE) {
+ GrammageType grammage = showerAxis_.getProjectedX(point);
+ int const maxBin = int(profile_.size() - 1);
+ int bin = grammage / dX_;
+ if (bin < 0) bin = 0;
+ if (bin > maxBin) bin = maxBin;
+
+ CORSIKA_LOGGER_TRACE(TOutput::getLogger(), "add local energy loss bin={} dE={} GeV ",
+ bin, dE / 1_GeV);
+
+ profile_[bin][static_cast<int>(ProfileIndex::Total)] += dE;
+ }
+
+ template <typename TOutput>
+ inline void EnergyLossWriter<TOutput>::write(GrammageType const Xstart,
+ GrammageType const Xend,
+ HEPEnergyType const dE) {
+ double const bstart = Xstart / dX_;
+ double const bend = Xend / dX_;
+
+ if (abs(bstart - floor(bstart + 0.5)) > 1e-2 ||
+ abs(bend - floor(bend + 0.5)) > 1e-2 || abs(bend - bstart - 1) > 1e-2) {
+ CORSIKA_LOGGER_ERROR(TOutput::getLogger(),
+ "CONEX and Corsika8 dX grammage binning are not the same! "
+ "Xstart={} Xend={} dX={}",
+ Xstart / 1_g * square(1_cm), Xend / 1_g * square(1_cm),
+ dX_ / 1_g * square(1_cm));
+ throw std::runtime_error(
+ "CONEX and Corsika8 dX grammage binning are not the same!");
+ }
+
+ int const bin = int((bend + bstart) / 2);
+ CORSIKA_LOGGER_TRACE(TOutput::getLogger(),
+ "add binned energy loss {} {} bin={} dE={} GeV ", bstart, bend,
+ bin, dE / 1_GeV);
+ profile_[bin][static_cast<int>(ProfileIndex::Total)] += dE;
+ }
+
+ template <typename TOutput>
+ inline HEPEnergyType EnergyLossWriter<TOutput>::getTotal() const {
+ HEPEnergyType tot = HEPEnergyType::zero();
+ for (Profile const& row : profile_)
+ tot += row.at(static_cast<int>(ProfileIndex::Total));
+ return tot;
+ }
+
+ template <typename TOutput>
+ inline YAML::Node EnergyLossWriter<TOutput>::getConfig() const {
+
+ YAML::Node node;
+
+ node["type"] = "EnergyLoss";
+ node["units"]["energy"] = "GeV";
+ node["units"]["grammage"] = "g/cm^2";
+ node["bin-size"] = dX_ / (1_g / square(1_cm));
+ node["nbins"] = nBins_;
+ node["grammage_threshold"] = dX_threshold_ / (1_g / square(1_cm));
+
+ // TODO: add shower axis to config
+
+ return node;
+ }
+
+ template <typename TOutput>
+ inline YAML::Node EnergyLossWriter<TOutput>::getSummary() const {
+
+ // determined Xmax and dEdXmax from quadratic interpolation
+ double maximum = 0;
+ unsigned int iMaximum = 0;
+ for (unsigned int i = 0; i < profile_.size() - 3; ++i) {
+ double value = (profile_[i + 0].at(static_cast<int>(ProfileIndex::Total)) +
+ profile_[i + 1].at(static_cast<int>(ProfileIndex::Total)) +
+ profile_[i + 2].at(static_cast<int>(ProfileIndex::Total))) /
+ 1_GeV;
+ if (value > maximum) {
+ maximum = value;
+ iMaximum = i;
+ }
+ }
+
+ double const dX = dX_ / 1_g * square(1_cm);
+
+ auto [Xmax, dEdXmax] = FindXmax::interpolateProfile(
+ dX * (0.5 + iMaximum), dX * (1.5 + iMaximum), dX * (2.5 + iMaximum),
+ profile_[iMaximum + 0].at(static_cast<int>(ProfileIndex::Total)) / 1_GeV,
+ profile_[iMaximum + 1].at(static_cast<int>(ProfileIndex::Total)) / 1_GeV,
+ profile_[iMaximum + 2].at(static_cast<int>(ProfileIndex::Total)) / 1_GeV);
+
+ YAML::Node summary;
+ summary["total"] = getTotal() / 1_GeV;
+ summary["Xmax"] = Xmax;
+ summary["dEdXmax"] = dEdXmax;
+ return summary;
+ }
+
+} // namespace corsika
diff --git a/corsika/detail/modules/writers/EnergyLossWriterParquet.inl b/corsika/detail/modules/writers/EnergyLossWriterParquet.inl
index 7148bdf8e43b01c77c7ba23503d29ea70a5bdf26..b7e82d669747150dd9ba73b7d1e21f939cef253d 100644
--- a/corsika/detail/modules/writers/EnergyLossWriterParquet.inl
+++ b/corsika/detail/modules/writers/EnergyLossWriterParquet.inl
@@ -18,14 +18,7 @@
namespace corsika {
- inline EnergyLossWriterParquet::EnergyLossWriterParquet(ShowerAxis const& showerAxis,
- GrammageType const dX,
- unsigned int const nBins,
- GrammageType const dX_threshold)
- : showerAxis_(showerAxis)
- , dX_(dX) // profile binning
- , nBins_(nBins)
- , dX_threshold_(dX_threshold) {}
+ inline EnergyLossWriterParquet::EnergyLossWriterParquet() {}
inline void EnergyLossWriterParquet::startOfLibrary(
boost::filesystem::path const& directory) {
@@ -46,160 +39,21 @@ namespace corsika {
output_.buildStreamer();
}
- inline void EnergyLossWriterParquet::startOfShower(unsigned int const) {
+ inline void EnergyLossWriterParquet::write(unsigned int const showerId,
+ GrammageType const grammage,
+ HEPEnergyType const total) {
- // initialize profile
- profile_.clear();
- profile_.resize(nBins_);
- }
-
- inline void EnergyLossWriterParquet::endOfShower(unsigned int const showerId) {
-
- int iRow = 0;
- for (Profile const& row : profile_) {
-
- double const dX = dX_ / 1_g * square(1_cm); // g/cm2
-
- // clang-format off
- *(output_.getWriter()) << showerId << static_cast<float>(iRow * dX)
- << static_cast<float>(row.at(static_cast<int>(ProfileIndex::Total)) / 1_GeV / dX)
- << parquet::EndRow;
- // clang-format on
- iRow++;
- }
- }
-
- inline void EnergyLossWriterParquet::endOfLibrary() { output_.closeStreamer(); }
+ double const dX = grammage / 1_g * square(1_cm); // g/cm2
- template <typename TTrack>
- inline void EnergyLossWriterParquet::write(TTrack const& track, Code const PID,
- HEPEnergyType const dE) {
-
- GrammageType grammageStart = showerAxis_.getProjectedX(track.getPosition(0));
- GrammageType grammageEnd = showerAxis_.getProjectedX(track.getPosition(1));
-
- if (grammageStart > grammageEnd) { // particle going upstream
- std::swap(grammageStart, grammageEnd);
- }
-
- GrammageType const deltaX = grammageEnd - grammageStart;
-
- if (deltaX < dX_threshold_) {
- write(track.getPosition(0), PID, dE);
- return;
- }
-
- // only register the range that is covered by the profile
- int const maxBin = int(profile_.size() - 1);
- int binStart = grammageStart / dX_;
- if (binStart < 0) binStart = 0;
- if (binStart > maxBin) binStart = maxBin;
- int binEnd = grammageEnd / dX_;
- if (binEnd < 0) binEnd = 0;
- if (binEnd > maxBin) binEnd = maxBin;
-
- CORSIKA_LOGGER_TRACE(getLogger(), "energy deposit of dE={} GeV between {} and {}",
- dE / 1_GeV, grammageStart / 1_g * square(1_cm),
- grammageEnd / 1_g * square(1_cm));
-
- auto energyCount = HEPEnergyType::zero();
-
- auto const factor = dE / deltaX;
- auto fill = [&](int const bin, GrammageType const weight) {
- auto const increment = factor * weight;
- profile_[bin][static_cast<int>(ProfileIndex::Total)] += increment;
- energyCount += increment;
-
- CORSIKA_LOGGER_TRACE(getLogger(), "filling bin={} with weight {} : dE={} GeV ", bin,
- weight, increment / 1_GeV);
- };
-
- // fill longitudinal profile
- if (binStart == binEnd) {
- fill(binStart, deltaX);
- } else {
- fill(binStart, ((1 + binStart) * dX_ - grammageStart));
- fill(binEnd, (grammageEnd - binEnd * dX_));
- for (int bin = binStart + 1; bin < binEnd; ++bin) { fill(bin, dX_); }
- }
-
- CORSIKA_LOGGER_TRACE(getLogger(), "total energy added to histogram: {} GeV ",
- energyCount / 1_GeV);
+ // and write the data into the column
+ *(output_.getWriter()) << showerId << static_cast<float>(dX)
+ << static_cast<float>(total / 1_GeV) << parquet::EndRow;
}
- inline void EnergyLossWriterParquet::write(Point const& point, Code const PID,
- HEPEnergyType const dE) {
- GrammageType grammage = showerAxis_.getProjectedX(point);
- int const maxBin = int(profile_.size() - 1);
- int bin = grammage / dX_;
- if (bin < 0) bin = 0;
- if (bin > maxBin) bin = maxBin;
+ inline void EnergyLossWriterParquet::startOfShower(unsigned int const) {}
- CORSIKA_LOGGER_TRACE(getLogger(), "add local energy loss bin={} dE={} GeV ", bin,
- dE / 1_GeV);
+ inline void EnergyLossWriterParquet::endOfShower(unsigned int const showerId) {}
- profile_[bin][static_cast<int>(ProfileIndex::Total)] += dE;
- }
-
- inline void EnergyLossWriterParquet::write(GrammageType const Xstart,
- GrammageType const Xend,
- HEPEnergyType const dE) {
- double const bstart = Xstart / dX_;
- double const bend = Xend / dX_;
-
- if (abs(bstart - floor(bstart + 0.5)) > 1e-2 ||
- abs(bend - floor(bend + 0.5)) > 1e-2 || abs(bend - bstart - 1) > 1e-2) {
- CORSIKA_LOGGER_ERROR(getLogger(),
- "CONEX and Corsika8 dX grammage binning are not the same! "
- "Xstart={} Xend={} dX={}",
- Xstart / 1_g * square(1_cm), Xend / 1_g * square(1_cm),
- dX_ / 1_g * square(1_cm));
- throw std::runtime_error(
- "CONEX and Corsika8 dX grammage binning are not the same!");
- }
-
- int const bin = int((bend + bstart) / 2);
- CORSIKA_LOGGER_TRACE(getLogger(), "add binned energy loss {} {} bin={} dE={} GeV ",
- bstart, bend, bin, dE / 1_GeV);
- profile_[bin][static_cast<int>(ProfileIndex::Total)] += dE;
- }
-
- inline HEPEnergyType EnergyLossWriterParquet::getTotal() const {
- HEPEnergyType tot = HEPEnergyType::zero();
- for (Profile const& row : profile_)
- tot += row.at(static_cast<int>(ProfileIndex::Total));
- return tot;
- }
-
- inline YAML::Node EnergyLossWriterParquet::getSummary() const {
-
- // determined Xmax and dEdXmax from quadratic interpolation
- double maximum = 0;
- unsigned int iMaximum = 0;
- for (unsigned int i = 0; i < profile_.size() - 3; ++i) {
- double value = (profile_[i + 0].at(static_cast<int>(ProfileIndex::Total)) +
- profile_[i + 1].at(static_cast<int>(ProfileIndex::Total)) +
- profile_[i + 2].at(static_cast<int>(ProfileIndex::Total))) /
- 1_GeV;
- if (value > maximum) {
- maximum = value;
- iMaximum = i;
- }
- }
-
- double const dX = dX_ / 1_g * square(1_cm);
-
- auto [Xmax, dEdXmax] = FindXmax::interpolateProfile(
- dX * (0.5 + iMaximum), dX * (1.5 + iMaximum), dX * (2.5 + iMaximum),
- profile_[iMaximum + 0].at(static_cast<int>(ProfileIndex::Total)) / 1_GeV,
- profile_[iMaximum + 1].at(static_cast<int>(ProfileIndex::Total)) / 1_GeV,
- profile_[iMaximum + 2].at(static_cast<int>(ProfileIndex::Total)) / 1_GeV);
-
- YAML::Node summary;
- summary["total"] = getTotal() / 1_GeV;
- summary["Xmax"] = Xmax;
- summary["dEdXmax"] = dEdXmax;
- return summary;
- }
+ inline void EnergyLossWriterParquet::endOfLibrary() { output_.closeStreamer(); }
} // namespace corsika
diff --git a/corsika/modules/writers/EnergyLossWriter.hpp b/corsika/modules/writers/EnergyLossWriter.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..8806acd8c5bc0330fc80243f9d95a06e5c8368d6
--- /dev/null
+++ b/corsika/modules/writers/EnergyLossWriter.hpp
@@ -0,0 +1,92 @@
+/*
+ * (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/media/ShowerAxis.hpp>
+#include <corsika/modules/writers/WriterOff.hpp>
+
+#include <vector>
+#include <array>
+
+namespace corsika {
+
+ template <typename TOutput = WriterOff>
+ class EnergyLossWriter : public TOutput {
+
+ enum class ProfileIndex { Total, Entries };
+ typedef std::array<HEPEnergyType, static_cast<int>(ProfileIndex::Entries)> Profile;
+
+ public:
+ /**
+ * Construct a new writer.
+ */
+ EnergyLossWriter(ShowerAxis const& axis,
+ GrammageType dX = 10_g / square(1_cm), // profile binning
+ unsigned int const nBins = 200, // number of bins
+ GrammageType dX_threshold = 0.0001_g /
+ square(1_cm)); // ignore too short tracks
+
+ void startOfLibrary(boost::filesystem::path const& directory) final override;
+
+ void startOfShower(unsigned int const showerId) final override;
+
+ void endOfShower(unsigned int const showerId) final override;
+
+ void endOfLibrary() final override;
+
+ /**
+ * Add continuous energy loss.
+ */
+ template <typename TTrack>
+ void write(TTrack const& track, Code const PID, HEPEnergyType const dE);
+
+ /**
+ * Add localized energy loss.
+ */
+ void write(Point const& point, Code const PID, HEPEnergyType const dE);
+
+ /**
+ * Add binned energy loss.
+ */
+ void write(GrammageType const Xstart, GrammageType const Xend,
+ HEPEnergyType const dE);
+
+ /**
+ * Get total observed energy loss.
+ *
+ * @return HEPEnergyType The total energy.
+ */
+ HEPEnergyType getTotal() const;
+
+ /**
+ * Return a summary.
+ */
+ YAML::Node getSummary() const;
+
+ /**
+ * Return the configuration of this output.
+ */
+ YAML::Node getConfig() const;
+
+ private:
+ ShowerAxis const& showerAxis_; ///< conversion between geometry and grammage
+ GrammageType dX_; ///< binning of profile.
+ unsigned int nBins_; ///< number of profile bins.
+ GrammageType dX_threshold_; ///< too short tracks are discarded.
+ std::vector<Profile> profile_; // longitudinal profile
+
+ }; // namespace corsika
+
+} // namespace corsika
+
+#include <corsika/detail/modules/writers/EnergyLossWriter.inl>
diff --git a/corsika/modules/writers/EnergyLossWriterParquet.hpp b/corsika/modules/writers/EnergyLossWriterParquet.hpp
index b200c8deb5177ce763539b9f6e8b9b202cfec8ca..5f2cbac28d1f389aa49a7b91dcb81918ad60acc3 100644
--- a/corsika/modules/writers/EnergyLossWriterParquet.hpp
+++ b/corsika/modules/writers/EnergyLossWriterParquet.hpp
@@ -21,33 +21,26 @@ namespace corsika {
class EnergyLossWriterParquet : public BaseOutput {
- enum class ProfileIndex { Total, Entries };
- typedef std::array<HEPEnergyType, static_cast<int>(ProfileIndex::Entries)> Profile;
-
public:
/**
* Construct a new writer.
*/
- EnergyLossWriterParquet(
- ShowerAxis const& axis,
- GrammageType dX = 10_g / square(1_cm), // profile binning
- unsigned int const nBins = 200, // number of bins
- GrammageType dX_threshold = 0.0001_g / square(1_cm)); // ignore too short tracks
+ EnergyLossWriterParquet();
/**
* Called at the start of each library.
*/
- void startOfLibrary(boost::filesystem::path const& directory) final override;
+ void startOfLibrary(boost::filesystem::path const& directory) override;
/**
* Called at the start of each shower.
*/
- void startOfShower(unsigned int const showerId) final override;
+ void startOfShower(unsigned int const showerId) override;
/**
* Called at the end of each shower.
*/
- void endOfShower(unsigned int const showerId) final override;
+ void endOfShower(unsigned int const showerId) override;
/**
* Called at the end of each library.
@@ -55,46 +48,17 @@ namespace corsika {
* This must also increment the run number since we override
* the default behaviour of BaseOutput.
*/
- void endOfLibrary() final override;
-
- /**
- * Add continuous energy loss.
- */
- template <typename TTrack>
- void write(TTrack const& track, Code const PID, HEPEnergyType const dE);
-
- /**
- * Add localized energy loss.
- */
- void write(Point const& point, Code const PID, HEPEnergyType const dE);
-
- /**
- * Add binned energy loss.
- */
- void write(GrammageType const Xstart, GrammageType const Xend,
- HEPEnergyType const dE);
-
- /**
- * Get total observed energy loss.
- *
- * @return HEPEnergyType The total energy.
- */
- HEPEnergyType getTotal() const;
+ void endOfLibrary() override;
/**
- * Returns library-wide summary.
+ * Write energy lost to the file
*/
- YAML::Node getSummary() const;
+ void write(unsigned int const showerId, GrammageType const grammage,
+ HEPEnergyType const total);
- public:
+ private:
ParquetStreamer output_; ///< The primary output file.
- ShowerAxis const& showerAxis_; ///< conversion between geometry and grammage
- GrammageType dX_; ///< binning of profile.
- unsigned int nBins_; ///< number of profile bins.
- GrammageType dX_threshold_; ///< too short tracks are discarded.
- std::vector<Profile> profile_; // longitudinal profile
-
}; // namespace corsika
} // namespace corsika
diff --git a/corsika/modules/writers/SubWriter.hpp b/corsika/modules/writers/SubWriter.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..d75017715114c97f9b845b4caeb63b8fbeae94cd
--- /dev/null
+++ b/corsika/modules/writers/SubWriter.hpp
@@ -0,0 +1,38 @@
+/*
+ * (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/Configurable.hpp>
+
+#include <utility>
+
+namespace corsika {
+
+ template <typename TOutput>
+ class SubWriter : public Configurable {
+
+ TOutput& output_;
+
+ public:
+ SubWriter(TOutput& output)
+ : output_(output) {}
+
+ virtual ~SubWriter() {}
+
+ /*
+ * Write data to the underlying writer.
+ */
+ template <typename... TArgs>
+ void write(TArgs&&... args) {
+ output_.write(std::forward<TArgs>(args)...);
+ }
+
+ }; // class SubWriter
+
+} // namespace corsika
diff --git a/examples/vertical_EAS.cpp b/examples/vertical_EAS.cpp
index 2023d527bd48ffda3af586c6d1cc3721b3b82ced..a8ee6af6f30226424ddf59eb4af3228158d339cd 100644
--- a/examples/vertical_EAS.cpp
+++ b/examples/vertical_EAS.cpp
@@ -41,9 +41,9 @@
#include <corsika/media/ShowerAxis.hpp>
#include <corsika/media/CORSIKA7Atmospheres.hpp>
-// #include <corsika/modules/writers/SubWriter.hpp>
-// #include <corsika/modules/writers/EnergyLossWriter.hpp>
-// #include <corsika/modules/writers/EnergyLossWriterParquet.hpp>
+#include <corsika/modules/writers/SubWriter.hpp>
+#include <corsika/modules/writers/EnergyLossWriter.hpp>
+#include <corsika/modules/writers/EnergyLossWriterParquet.hpp>
#include <corsika/modules/BetheBlochPDG.hpp>
#include <corsika/modules/writers/BetheBlochPDGWriterParquet.hpp>
#include <corsika/modules/LongitudinalProfile.hpp>
@@ -203,13 +203,23 @@ int main(int argc, char** argv) {
// output.add("dEdX", dEdX_output);
// register profile output
+ // construct the overall energy loss writer and register it
+ EnergyLossWriter<EnergyLossWriterParquet> dEdX{showerAxis};
+ output.add("energyloss", dEdX);
+
+ // construct the continuous energy loss model
+ BetheBlochPDG<SubWriter<decltype(dEdX)>> emContinuous{dEdX};
+ // output.add("bethebloch", emContinuous);
+
+ // construct a particle cut
+ ParticleCut<SubWriter<decltype(dEdX)>> cut{60_GeV, 60_GeV, 60_GeV, 60_GeV, true, dEdX};
+ // output.add("cuts", cut);
+
// setup longitudinal profile
- LongitudinalProfile<LongitudinalProfileWriterParquet> profile{showerAxis};
+ LongitudinalProfile<corsika::LongitudinalProfileWriterParquet> profile{showerAxis};
output.add("profile", profile);
- // register ground particle output
- // ParticleWriterParquet particleOutput;
-
+ // create a track writer and register it with the output manager
TrackWriter<TrackWriterParquet> trackWriter;
output.add("tracks", trackWriter);
@@ -246,9 +256,6 @@ int main(int argc, char** argv) {
decaySibyll.printDecayConfig();
- ParticleCut<ParticleCutWriterParquet> cut{60_GeV, 60_GeV, 60_GeV, 60_GeV, true};
- output.add("cuts", cut);
-
corsika::urqmd::UrQMD urqmd;
InteractionCounter urqmdCounted{urqmd};
StackInspector<setup::Stack> stackInspect(50000, false, E0);
@@ -281,16 +288,11 @@ int main(int argc, char** argv) {
obsPlane, DirectionVector(rootCS, {1., 0., 0.}), particleOutput};
output.add("particles", observationLevel);
- // EnergyLossWriter<EnergyLossWriterParquet> dEdX{showerAxis, 10_g / square(1_cm), 200};
- BetheBlochPDG<BetheBlochPDGWriterParquet> emContinuous;
- output.add("bethebloch", emContinuous);
- // dEdX.add(emContinuous);
-
// auto sequence = make_sequence(stackInspect, hadronSequence, decaySequence,
// emContinuous,
// cut, trackWriter, observationLevel, longprof);
auto sequence = make_sequence(stackInspect, hadronSequence, decaySequence, emContinuous,
- cut, trackWriter, observationLevel);
+ cut, trackWriter, observationLevel, profile);
// define air shower object, run simulation
setup::Tracking tracking;