diff --git a/corsika/detail/modules/ObservationPlane.inl b/corsika/detail/modules/ObservationPlane.inl
index eba312b29975cf5ba031e74533ea2457dad8aa5b..8f4a658b504d896f674943f4729cc2a30336329f 100644
--- a/corsika/detail/modules/ObservationPlane.inl
+++ b/corsika/detail/modules/ObservationPlane.inl
@@ -18,9 +18,7 @@ namespace corsika {
, plane_(obsPlane)
, xAxis_(x_axis.normalized())
, yAxis_(obsPlane.getNormal().cross(xAxis_))
- , deleteOnHit_(deleteOnHit)
- , energy_ground_(0_GeV)
- , count_ground_(0) {}
+ , deleteOnHit_(deleteOnHit) {}
template <typename TTracking, typename TOutput>
template <typename TParticle, typename TTrajectory>
@@ -60,8 +58,6 @@ namespace corsika {
CORSIKA_LOG_TRACE("Particle detected absorbed={}", deleteOnHit_);
if (deleteOnHit_) {
- count_ground_++;
- energy_ground_ += energy;
return ProcessReturn::ParticleAbsorbed;
} else {
return ProcessReturn::Ok;
@@ -94,17 +90,6 @@ namespace corsika {
return trajectory.getLength(fractionOfIntersection);
}
- template <typename TTracking, typename TOutput>
- inline void ObservationPlane<TTracking, TOutput>::showResults() const {
- CORSIKA_LOG_INFO(
- "\n ******************************\n"
- " ObservationPlane: \n"
- " energy an ground (GeV) : {}\n"
- " no. of particles at ground : {}\n"
- " ******************************",
- energy_ground_ / 1_GeV, count_ground_);
- }
-
template <typename TTracking, typename TOutput>
inline YAML::Node ObservationPlane<TTracking, TOutput>::getConfig() const {
using namespace units::si;
@@ -148,10 +133,4 @@ namespace corsika {
return node;
}
- template <typename TTracking, typename TOutput>
- inline void ObservationPlane<TTracking, TOutput>::reset() {
- energy_ground_ = 0_GeV;
- count_ground_ = 0;
- }
-
} // namespace corsika
diff --git a/corsika/detail/modules/StackInspector.inl b/corsika/detail/modules/StackInspector.inl
index 06705ac8f6be0892750e8c3942e4710d273b6dc8..f801abf624705b5bb5891bd34304c3ed0570936b 100644
--- a/corsika/detail/modules/StackInspector.inl
+++ b/corsika/detail/modules/StackInspector.inl
@@ -54,20 +54,24 @@ namespace corsika {
}
}
- auto const now = std::chrono::system_clock::now();
- std::chrono::duration<double> const elapsed_seconds = now - StartTime_;
- std::time_t const now_time = std::chrono::system_clock::to_time_t(now);
+ std::chrono::system_clock::time_point const now = std::chrono::system_clock::now();
+ std::chrono::duration<double> const elapsed_seconds = now - StartTime_; // seconds
auto const dE = E0_ - Etot;
if (dE < dE_threshold_) return;
double const progress = dE / E0_;
+ // for printout
+ std::time_t const now_time = std::chrono::system_clock::to_time_t(now);
std::time_t const start_time = std::chrono::system_clock::to_time_t(StartTime_);
if (progress > 0) {
double const eta_seconds = elapsed_seconds.count() / progress;
- std::time_t const eta_time = std::chrono::system_clock::to_time_t(
- StartTime_ + std::chrono::seconds((int)eta_seconds));
+ std::chrono::system_clock::time_point const eta =
+ StartTime_ + std::chrono::seconds((int)eta_seconds);
+
+ // for printout
+ std::time_t const eta_time = std::chrono::system_clock::to_time_t(eta);
int const yday0 = std::localtime(&start_time)->tm_yday;
int const yday1 = std::localtime(&eta_time)->tm_yday;
diff --git a/corsika/detail/modules/urqmd/UrQMD.inl b/corsika/detail/modules/urqmd/UrQMD.inl
index 9e070f1c62df8860f9bbb4d89df26f9c87728e09..1772df106dcb04d8744eac83959e833f5fa27f9a 100644
--- a/corsika/detail/modules/urqmd/UrQMD.inl
+++ b/corsika/detail/modules/urqmd/UrQMD.inl
@@ -13,6 +13,7 @@
#include <corsika/framework/core/ParticleProperties.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/core/EnergyMomentumOperations.hpp>
#include <corsika/framework/geometry/QuantityVector.hpp>
#include <corsika/framework/geometry/Vector.hpp>
#include <corsika/framework/utility/COMBoost.hpp>
@@ -280,9 +281,11 @@ namespace corsika::urqmd {
CORSIKA_LOG_DEBUG(" {} {} {} ", i, code, momentum.getComponents());
HEPEnergyType const mass = get_mass(code);
- HEPEnergyType Ekin = sqrt(momentum.getSquaredNorm() + mass * mass) - mass;
+ HEPEnergyType const Ekin = calculate_kinetic_energy(momentum.getNorm(), mass);
if (Ekin <= 0_GeV) {
- CORSIKA_LOG_WARN("Negative kinetic energy {} {}. Skipping.", code, Ekin);
+ if (Ekin < 0_GeV) {
+ CORSIKA_LOG_WARN("Negative kinetic energy {} {}. Skipping.", code, Ekin);
+ }
view.addSecondary(
std::make_tuple(code, 0_eV, DirectionVector{originalCS, {0, 0, 0}}));
} else {
diff --git a/corsika/detail/modules/writers/EnergyLossWriter.inl b/corsika/detail/modules/writers/EnergyLossWriter.inl
index f4ef75adab355276b34e0e5754cf6c913c4c5eed..e4101afeb2145f70e9f3069bd43b1a7c380b4a6e 100644
--- a/corsika/detail/modules/writers/EnergyLossWriter.inl
+++ b/corsika/detail/modules/writers/EnergyLossWriter.inl
@@ -23,7 +23,8 @@ namespace corsika {
GrammageType dX,
unsigned int const nBins,
GrammageType dX_threshold)
- : showerAxis_(axis)
+ : TOutput(dEdX_output::ProfileIndexNames)
+ , showerAxis_(axis)
, dX_(dX)
, nBins_(nBins)
, dX_threshold_(dX_threshold) {}
@@ -51,9 +52,9 @@ namespace corsika {
inline void EnergyLossWriter<TOutput>::endOfShower(unsigned int const showerId) {
int iRow{0};
- for (Profile const& row : profile_) {
+ for (dEdX_output::Profile const& row : profile_) {
// here: write to underlying writer (e.g. parquet)
- TOutput::write(showerId, iRow * dX_, row.at(static_cast<int>(ProfileIndex::Total)));
+ TOutput::write(showerId, iRow * dX_, row);
iRow++;
}
@@ -100,7 +101,7 @@ namespace corsika {
CORSIKA_LOGGER_TRACE(TOutput::getLogger(),
"filling bin={} with weight {} : dE={} GeV ", bin, weight,
increment / 1_GeV);
- profile_[bin][static_cast<int>(ProfileIndex::Total)] += increment;
+ profile_[bin][static_cast<int>(dEdX_output::ProfileIndex::Total)] += increment;
energyCount += increment;
};
@@ -129,7 +130,7 @@ namespace corsika {
CORSIKA_LOGGER_TRACE(TOutput::getLogger(), "add local energy loss bin={} dE={} GeV ",
bin, dE / 1_GeV);
- profile_[bin][static_cast<int>(ProfileIndex::Total)] += dE;
+ profile_[bin][static_cast<int>(dEdX_output::ProfileIndex::Total)] += dE;
}
template <typename TOutput>
@@ -141,11 +142,12 @@ namespace corsika {
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={} g/cm2",
- Xstart / 1_g * square(1_cm), Xend / 1_g * square(1_cm),
- dX_ / 1_g * square(1_cm));
+ CORSIKA_LOGGER_ERROR(
+ TOutput::getLogger(),
+ "CascadeEquation (CONEX) and Corsika8 dX grammage binning are not the same! "
+ "Xstart={} Xend={} dX={} g/cm2",
+ 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!");
}
@@ -154,14 +156,14 @@ namespace corsika {
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;
+ profile_[bin][static_cast<int>(dEdX_output::ProfileIndex::Total)] += dE;
}
template <typename TOutput>
- inline HEPEnergyType EnergyLossWriter<TOutput>::getTotal() const {
+ inline HEPEnergyType EnergyLossWriter<TOutput>::getEnergyLost() const {
HEPEnergyType tot = HEPEnergyType::zero();
- for (Profile const& row : profile_)
- tot += row.at(static_cast<int>(ProfileIndex::Total));
+ for (dEdX_output::Profile const& row : profile_)
+ tot += row.at(static_cast<int>(dEdX_output::ProfileIndex::Total));
return tot;
}
@@ -177,8 +179,6 @@ namespace corsika {
node["nbins"] = nBins_;
node["grammage_threshold"] = dX_threshold_ / (1_g / square(1_cm));
- //! \todo add shower axis to config
-
return node;
}
@@ -187,12 +187,13 @@ namespace corsika {
// 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;
+ size_t iMaximum = 0;
+ for (size_t i = 0; i < profile_.size() - 3; ++i) {
+ double value =
+ (profile_[i + 0].at(static_cast<int>(dEdX_output::ProfileIndex::Total)) +
+ profile_[i + 1].at(static_cast<int>(dEdX_output::ProfileIndex::Total)) +
+ profile_[i + 2].at(static_cast<int>(dEdX_output::ProfileIndex::Total))) /
+ 1_GeV;
if (value > maximum) {
maximum = value;
iMaximum = i;
@@ -203,12 +204,15 @@ namespace corsika {
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);
+ profile_[iMaximum + 0].at(static_cast<int>(dEdX_output::ProfileIndex::Total)) /
+ 1_GeV,
+ profile_[iMaximum + 1].at(static_cast<int>(dEdX_output::ProfileIndex::Total)) /
+ 1_GeV,
+ profile_[iMaximum + 2].at(static_cast<int>(dEdX_output::ProfileIndex::Total)) /
+ 1_GeV);
YAML::Node summary;
- summary["total"] = getTotal() / 1_GeV;
+ summary["sum_dEdX"] = getEnergyLost() / 1_GeV;
summary["Xmax"] = Xmax;
summary["dEdXmax"] = dEdXmax;
return summary;
diff --git a/corsika/detail/modules/writers/EnergyLossWriterParquet.inl b/corsika/detail/modules/writers/EnergyLossWriterParquet.inl
index 1680e421ac51725faab57a3760c43fe6c1a8b6d9..55b301f509bc094a878e529e8ebb58360dc6568a 100644
--- a/corsika/detail/modules/writers/EnergyLossWriterParquet.inl
+++ b/corsika/detail/modules/writers/EnergyLossWriterParquet.inl
@@ -18,9 +18,13 @@
namespace corsika {
- inline EnergyLossWriterParquet::EnergyLossWriterParquet() {}
+ template <size_t NColumns>
+ inline EnergyLossWriterParquet<NColumns>::EnergyLossWriterParquet(
+ std::array<const char*, NColumns> const& colNames)
+ : columns_(colNames) {}
- inline void EnergyLossWriterParquet::startOfLibrary(
+ template <size_t NColumns>
+ inline void EnergyLossWriterParquet<NColumns>::startOfLibrary(
boost::filesystem::path const& directory) {
// setup the streamer
@@ -32,28 +36,38 @@ namespace corsika {
// build the schema
output_.addField("X", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
parquet::ConvertedType::NONE);
- output_.addField("total", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
- parquet::ConvertedType::NONE);
+ for (auto const& col : columns_) {
+ output_.addField(col, parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
+ parquet::ConvertedType::NONE);
+ }
// and build the streamer
output_.buildStreamer();
}
- inline void EnergyLossWriterParquet::write(unsigned int const showerId,
- GrammageType const grammage,
- HEPEnergyType const total) {
-
- double const dX = grammage / 1_g * square(1_cm); // g/cm2
+ template <size_t NColumns>
+ inline void EnergyLossWriterParquet<NColumns>::write(
+ unsigned int const showerId, GrammageType const grammage,
+ std::array<HEPEnergyType, NColumns> const& data) {
// and write the data into the column
- *(output_.getWriter()) << showerId << static_cast<float>(dX)
- << static_cast<float>(total / 1_GeV) << parquet::EndRow;
+ *(output_.getWriter()) << showerId
+ << static_cast<float>(grammage / 1_g * square(1_cm));
+ for (HEPEnergyType const dedx : data) {
+ *(output_.getWriter()) << static_cast<float>(dedx / 1_GeV);
+ }
+ *(output_.getWriter()) << parquet::EndRow;
}
- inline void EnergyLossWriterParquet::startOfShower(unsigned int const) {}
+ template <size_t NColumns>
+ inline void EnergyLossWriterParquet<NColumns>::startOfShower(unsigned int const) {}
- inline void EnergyLossWriterParquet::endOfShower(unsigned int const) {}
+ template <size_t NColumns>
+ inline void EnergyLossWriterParquet<NColumns>::endOfShower(unsigned int const) {}
- inline void EnergyLossWriterParquet::endOfLibrary() { output_.closeStreamer(); }
+ template <size_t NColumns>
+ inline void EnergyLossWriterParquet<NColumns>::endOfLibrary() {
+ output_.closeStreamer();
+ }
} // namespace corsika
diff --git a/corsika/detail/modules/writers/LongitudinalProfileWriterParquet.inl b/corsika/detail/modules/writers/LongitudinalProfileWriterParquet.inl
index 4873bc486cbcba7d307eed92183ae76d7820516c..2e5801fee1882cc0ee44c522ae63fc1b205a5f4b 100644
--- a/corsika/detail/modules/writers/LongitudinalProfileWriterParquet.inl
+++ b/corsika/detail/modules/writers/LongitudinalProfileWriterParquet.inl
@@ -19,14 +19,13 @@
namespace corsika {
- inline LongitudinalProfileWriterParquet::LongitudinalProfileWriterParquet(
- ShowerAxis const& showerAxis, GrammageType const dX, unsigned int const nBins)
- : output_()
- , showerAxis_(showerAxis)
- , dX_(dX) // profile binning
- , nBins_(nBins) {}
-
- inline void LongitudinalProfileWriterParquet::startOfLibrary(
+ template <size_t NColumns>
+ inline LongitudinalProfileWriterParquet<NColumns>::LongitudinalProfileWriterParquet(
+ std::array<const char*, NColumns> const& columns)
+ : columns_(columns) {}
+
+ template <size_t NColumns>
+ inline void LongitudinalProfileWriterParquet<NColumns>::startOfLibrary(
boost::filesystem::path const& directory) {
// setup the streamer
output_.initStreamer((directory / "profile.parquet").string());
@@ -37,155 +36,40 @@ namespace corsika {
// build the schema
output_.addField("X", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
parquet::ConvertedType::NONE);
- output_.addField("charged", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
- parquet::ConvertedType::NONE);
- output_.addField("hadron", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
- parquet::ConvertedType::NONE);
- output_.addField("muminus", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
- parquet::ConvertedType::NONE);
- output_.addField("muplus", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
- parquet::ConvertedType::NONE);
- output_.addField("photon", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
- parquet::ConvertedType::NONE);
- output_.addField("electron", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
- parquet::ConvertedType::NONE);
- output_.addField("positron", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
- parquet::ConvertedType::NONE);
+ for (auto const& col : columns_) {
+ output_.addField(col, parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
+ parquet::ConvertedType::NONE);
+ }
// and build the streamer
output_.buildStreamer();
}
- inline void LongitudinalProfileWriterParquet::startOfShower(unsigned int const) {
- // initialize profile
- profile_.clear();
- profile_.resize(nBins_);
- }
+ template <size_t NColumns>
+ inline void LongitudinalProfileWriterParquet<NColumns>::startOfShower(
+ unsigned int const) {}
- inline void LongitudinalProfileWriterParquet::endOfShower(unsigned int const showerId) {
- int iRow = 0;
- for (ProfileData 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::Charged)))
- << static_cast<float>(row.at(static_cast<int>(ProfileIndex::Hadron)))
- << static_cast<float>(row.at(static_cast<int>(ProfileIndex::MuMinus)))
- << static_cast<float>(row.at(static_cast<int>(ProfileIndex::MuPlus)))
- << static_cast<float>(row.at(static_cast<int>(ProfileIndex::Photon)))
- << static_cast<float>(row.at(static_cast<int>(ProfileIndex::Electron)))
- << static_cast<float>(row.at(static_cast<int>(ProfileIndex::Positron)))
- << parquet::EndRow;
- // clang-format on
- ++iRow;
- }
- }
+ template <size_t NColumns>
+ inline void LongitudinalProfileWriterParquet<NColumns>::endOfShower(
+ unsigned int const) {}
- inline void LongitudinalProfileWriterParquet::endOfLibrary() {
+ template <size_t NColumns>
+ inline void LongitudinalProfileWriterParquet<NColumns>::endOfLibrary() {
output_.closeStreamer();
}
- template <typename TTrack>
- inline void LongitudinalProfileWriterParquet::write(TTrack const& track, Code const pid,
- double const weight) {
- GrammageType const grammageStart = showerAxis_.getProjectedX(track.getPosition(0));
- GrammageType const grammageEnd = showerAxis_.getProjectedX(track.getPosition(1));
-
- // Note: particle may go also "upward", thus, grammageEnd<grammageStart
- int const binStart = std::ceil(grammageStart / dX_);
- int const binEnd = std::floor(grammageEnd / dX_);
-
- for (int bin = binStart; bin <= binEnd; ++bin) {
- if (pid == Code::Photon) {
- profile_.at(bin)[static_cast<int>(ProfileIndex::Photon)] += weight;
- } else if (pid == Code::Positron) {
- profile_.at(bin)[static_cast<int>(ProfileIndex::Positron)] += weight;
- } else if (pid == Code::Electron) {
- profile_.at(bin)[static_cast<int>(ProfileIndex::Electron)] += weight;
- } else if (pid == Code::MuPlus) {
- profile_.at(bin)[static_cast<int>(ProfileIndex::MuPlus)] += weight;
- } else if (pid == Code::MuMinus) {
- profile_.at(bin)[static_cast<int>(ProfileIndex::MuMinus)] += weight;
- } else if (is_hadron(pid)) {
- profile_.at(bin)[static_cast<int>(ProfileIndex::Hadron)] += weight;
- }
- if (is_charged(pid)) {
- profile_[bin][static_cast<int>(ProfileIndex::Charged)] += weight;
- }
- }
- }
-
- inline void LongitudinalProfileWriterParquet::write(GrammageType const Xstart,
- GrammageType const Xend,
- Code const pid,
- double const weight) {
- 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);
-
- if (pid == Code::Photon) {
- profile_.at(bin)[static_cast<int>(ProfileIndex::Photon)] += weight;
- } else if (pid == Code::Positron) {
- profile_.at(bin)[static_cast<int>(ProfileIndex::Positron)] += weight;
- } else if (pid == Code::Electron) {
- profile_.at(bin)[static_cast<int>(ProfileIndex::Electron)] += weight;
- } else if (pid == Code::MuPlus) {
- profile_.at(bin)[static_cast<int>(ProfileIndex::MuPlus)] += weight;
- } else if (pid == Code::MuMinus) {
- profile_.at(bin)[static_cast<int>(ProfileIndex::MuMinus)] += weight;
- } else if (is_hadron(pid)) {
- profile_.at(bin)[static_cast<int>(ProfileIndex::Hadron)] += weight;
- }
- if (is_charged(pid)) {
- profile_[bin][static_cast<int>(ProfileIndex::Charged)] += weight;
- }
- }
+ template <size_t NColumns>
+ inline void LongitudinalProfileWriterParquet<NColumns>::write(
+ unsigned int const showerId, GrammageType const grammage,
+ std::array<double, NColumns> const& data) {
- inline YAML::Node LongitudinalProfileWriterParquet::getSummary() const {
- // determined Xmax and dEdXmax from quadratic interpolation
-
- YAML::Node summary;
-
- for (int index = 0; index < static_cast<int>(ProfileIndex::Entries); ++index) {
- // first find highest 3-boxcar sum
- double maximum = 0;
- unsigned int iMaximum = 0;
- for (unsigned int i = 0; i < profile_.size() - 3; ++i) {
- double value = profile_[i + 0].at(index) + profile_[i + 1].at(index) +
- profile_[i + 2].at(index);
- if (value > maximum) {
- maximum = value;
- iMaximum = i;
- }
- }
-
- double const dX = dX_ / 1_g * square(1_cm);
-
- // quadratic interpolation of maximum in 3 highest points
- auto [Xmax, Nmax] = FindXmax::interpolateProfile(
- dX * (0.5 + iMaximum), dX * (1.5 + iMaximum), dX * (2.5 + iMaximum),
- profile_[iMaximum + 0].at(index), profile_[iMaximum + 1].at(index),
- profile_[iMaximum + 2].at(index));
-
- std::string const name = ProfileIndexNames[index];
- summary[name]["Xmax"] = Xmax;
- summary[name]["Nmax"] = Nmax;
+ // and write the data into the column
+ *(output_.getWriter()) << showerId
+ << static_cast<float>(grammage / 1_g * square(1_cm));
+ for (double const weight : data) {
+ *(output_.getWriter()) << static_cast<float>(weight);
}
- return summary;
+ *(output_.getWriter()) << parquet::EndRow;
}
} // namespace corsika
diff --git a/corsika/detail/modules/writers/ParticleWriterParquet.inl b/corsika/detail/modules/writers/ParticleWriterParquet.inl
index fa26e368dc956a08878885d76b18d490041c3109..74e8b45411b3c90f563d16baf3724f1787182cb0 100644
--- a/corsika/detail/modules/writers/ParticleWriterParquet.inl
+++ b/corsika/detail/modules/writers/ParticleWriterParquet.inl
@@ -44,6 +44,10 @@ namespace corsika {
output_.buildStreamer();
showerId_ = 0;
+ }
+
+ inline void ParticleWriterParquet::startOfShower(unsigned int const showerId) {
+ showerId_ = showerId;
totalEnergy_ = 0_eV;
countHadrons_ = 0;
countOthers_ = 0;
@@ -51,10 +55,6 @@ namespace corsika {
countMuons_ = 0;
}
- inline void ParticleWriterParquet::startOfShower(unsigned int const showerId) {
- showerId_ = showerId;
- }
-
inline void ParticleWriterParquet::endOfShower(unsigned int const) {}
inline void ParticleWriterParquet::endOfLibrary() { output_.closeStreamer(); }
@@ -103,6 +103,7 @@ namespace corsika {
*/
YAML::Node ParticleWriterParquet::getSummary() const {
YAML::Node summary;
+ summary["Eground"] = totalEnergy_ / 1_GeV;
summary["hadrons"] = countHadrons_;
summary["muons"] = countMuons_;
summary["em"] = countEM_;
diff --git a/corsika/modules/ObservationPlane.hpp b/corsika/modules/ObservationPlane.hpp
index aec9f663164a23bfe9216629d12c95539c44c962..21ceb790c2a8d90be35a9a6bcb6cb1078f2388cd 100644
--- a/corsika/modules/ObservationPlane.hpp
+++ b/corsika/modules/ObservationPlane.hpp
@@ -62,9 +62,6 @@ namespace corsika {
template <typename TParticle, typename TTrajectory>
LengthType getMaxStepLength(TParticle const&, TTrajectory const& vTrajectory);
- void showResults() const;
- void reset();
- HEPEnergyType getEnergyGround() const { return energy_ground_; }
YAML::Node getConfig() const;
private:
@@ -72,8 +69,6 @@ namespace corsika {
DirectionVector const xAxis_;
DirectionVector const yAxis_;
bool const deleteOnHit_;
- HEPEnergyType energy_ground_;
- unsigned int count_ground_;
};
//! @}
} // namespace corsika
diff --git a/corsika/modules/StackInspector.hpp b/corsika/modules/StackInspector.hpp
index 6723ec3d9bedb0ff3097463133cc5d1a5e8dbe7c..6e723a6ced41be1e34a970991d3e55fc40c6af88 100644
--- a/corsika/modules/StackInspector.hpp
+++ b/corsika/modules/StackInspector.hpp
@@ -48,7 +48,7 @@ namespace corsika {
bool ReportStack_;
HEPEnergyType E0_;
const HEPEnergyType dE_threshold_ = 1_eV;
- decltype(std::chrono::system_clock::now()) StartTime_;
+ std::chrono::system_clock::time_point StartTime_;
};
} // namespace corsika
diff --git a/corsika/modules/writers/EnergyLossWriter.hpp b/corsika/modules/writers/EnergyLossWriter.hpp
index 5436cb5219d9538d450993f62b697cdeeaea8e43..58cc69ebcb034600ba8755a20d57d353b50af71a 100644
--- a/corsika/modules/writers/EnergyLossWriter.hpp
+++ b/corsika/modules/writers/EnergyLossWriter.hpp
@@ -18,8 +18,13 @@
#include <vector>
#include <array>
+/**
+ * @file EnergyLossWriter.hpp
+ */
+
namespace corsika {
+ // clang-format-off
/**
* The energy loss writer can be used to pool several energy loss processes into one
* output file/stream.
@@ -27,6 +32,7 @@ namespace corsika {
* Typically many processes/modules can lead to energy losses in the shower. The
* EnergyLossWriter can be used in combination with the SubWriter class to collect all
* of them into a single output stream:
+ *
* \code {.cpp}
* # showerAxis must be a ShowerAxis object
* # the X binning can be specified.
@@ -38,16 +44,68 @@ namespace corsika {
* ...
* \endcode
*
+ * The EnergyLossWriter processes data on single-particle-level. The final output
+ * writer, e.g. EnergyLossWriterParquet, processes data on profile-level (bins in X).
* The default output option is parquet format.
*
* @tparam TOutput
*/
+ // clang-format-on
- template <typename TOutput = EnergyLossWriterParquet>
- class EnergyLossWriter : public TOutput {
+ /**
+ * Local helper namespace to store number and names of dEdX profile columns.
+ */
+ namespace dEdX_output {
+ /**
+ * Definition of longitudinal profile columns.
+ */
enum class ProfileIndex { Total, Entries };
- typedef std::array<HEPEnergyType, static_cast<int>(ProfileIndex::Entries)> Profile;
+
+ /**
+ * Number of columns (static).
+ */
+ size_t constexpr NColumns = static_cast<int>(ProfileIndex::Entries);
+
+ /**
+ * Names of columns in output.
+ */
+ static std::array<char const*, NColumns> constexpr ProfileIndexNames{{"total"}};
+
+ /**
+ * Data type to store column data.
+ */
+ typedef std::array<HEPEnergyType, NColumns> Profile;
+
+ } // namespace dEdX_output
+
+ /**
+ * The EnergyLossWriter can be used to pool the dEdX energy loss of several
+ * processes/modules into one output file/stream.
+ *
+ * Typically several processes/modules can lead to energy losses along the shower axis
+ * in the shower. The EnergyLossWriter can be used in combination with the SubWriter
+ * class to collect all of them into a single output stream:
+ *
+ * \code {.cpp}
+ * # showerAxis must be a ShowerAxis object
+ * # the X binning can be specified.
+ * EnergyLossWriter dEdX{showerAxis, 10_g / square(1_cm), 200};
+ * # add to OutputManager:
+ * output.add("energyloss", dEdX);
+ * # add SubWriters, e.g. BetheBlochPDG, CONEX:
+ * BetheBlochPDG<SubWriter<decltype(dEdX)>> long{dEdX};
+ * CONEXhybrid<SubWriter<decltype(dEdX)>> conex{..., dEdX};
+ * ...
+ * \endcode
+ *
+ * The default output option is parquet format.
+ *
+ * @tparam TOutput
+ */
+
+ template <typename TOutput = EnergyLossWriterParquet<dEdX_output::NColumns>>
+ class EnergyLossWriter : public TOutput {
public:
/**
@@ -89,7 +147,7 @@ namespace corsika {
*
* @return HEPEnergyType The total energy.
*/
- HEPEnergyType getTotal() const;
+ HEPEnergyType getEnergyLost() const;
/**
* Return a summary.
@@ -106,7 +164,7 @@ namespace corsika {
GrammageType dX_; ///< binning of profile.
size_t nBins_; ///< number of profile bins.
GrammageType dX_threshold_; ///< too short tracks are discarded.
- std::vector<Profile> profile_; // longitudinal profile
+ std::vector<dEdX_output::Profile> profile_; // longitudinal profile
}; // namespace corsika
diff --git a/corsika/modules/writers/EnergyLossWriterParquet.hpp b/corsika/modules/writers/EnergyLossWriterParquet.hpp
index 3edc4cdd614e30e26a26135871791d88c26b874a..852696982645d57533f5fb547c0b55991a463d24 100644
--- a/corsika/modules/writers/EnergyLossWriterParquet.hpp
+++ b/corsika/modules/writers/EnergyLossWriterParquet.hpp
@@ -14,18 +14,30 @@
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/media/ShowerAxis.hpp>
-#include <vector>
#include <array>
namespace corsika {
+ /**
+ * The actual writer to save dEdX data to disk.
+ *
+ * The purpose of this class is not to collect single-particle-level energy loss data.
+ * But to write entire binned profiles to disk at the end of a shower.
+ * To fill the shower data, you have to use EnergyLossWriter in combination with
+ * SubWriter. The EnergyLossWriterParquet is the default output mode of the
+ * EnergyLossWriter.
+ *
+ * @tparam NColumn -- the number of columns written to output. column names and data
+ * must be provided consistently.
+ */
+ template <size_t NColumns>
class EnergyLossWriterParquet : public BaseOutput {
public:
/**
* Construct a new writer.
*/
- EnergyLossWriterParquet();
+ EnergyLossWriterParquet(std::array<const char*, NColumns> const& colNames);
/**
* Called at the start of each library.
@@ -54,10 +66,11 @@ namespace corsika {
* Write energy lost to the file.
*/
void write(unsigned int const showerId, GrammageType const grammage,
- HEPEnergyType const total);
+ std::array<HEPEnergyType, NColumns> const& data);
private:
- ParquetStreamer output_; ///< The primary output file.
+ std::array<const char*, NColumns> columns_; ///< column names
+ ParquetStreamer output_; ///< The primary output file.
}; // namespace corsika
diff --git a/corsika/modules/writers/LongitudinalProfileWriterParquet.hpp b/corsika/modules/writers/LongitudinalProfileWriterParquet.hpp
index 0c4570c71eaa7940b4d6d30a093b33704c1c39c7..bba7000ff30f73260740c66ab35a3539f97259b6 100644
--- a/corsika/modules/writers/LongitudinalProfileWriterParquet.hpp
+++ b/corsika/modules/writers/LongitudinalProfileWriterParquet.hpp
@@ -12,56 +12,44 @@
#include <corsika/output/ParquetStreamer.hpp>
#include <corsika/framework/core/ParticleProperties.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
-#include <corsika/media/ShowerAxis.hpp>
-#include <vector>
#include <array>
-#include <string>
namespace corsika {
+ /**
+ * The actual writer to save longitudinal profile data to disk.
+ *
+ * The purpose of this class is not to collect single-particle-level profile data,
+ * but to write entire binned profiles to disk at the end of a shower.
+ * To fill the shower data, you have to use ProfileWriter in combination with
+ * SubWriter. The LongitudinalProfileWriterParquet is the default output mode of the
+ * LongitudinalWriter.
+ */
+
+ template <size_t NColumns>
class LongitudinalProfileWriterParquet : public BaseOutput {
- enum class ProfileIndex {
- Charged,
- Hadron,
- Photon,
- Electron,
- Positron,
- MuPlus,
- MuMinus,
- Entries
- };
-
- static std::array<
- char const*, static_cast<int>(ProfileIndex::Entries)> constexpr ProfileIndexNames{
- {"charged", "hadron", "photon", "electron", "positron", "muplus", "muminus"}};
-
- typedef std::array<double, static_cast<int>(ProfileIndex::Entries)> ProfileData;
-
public:
/**
* Construct a new writer.
*/
- LongitudinalProfileWriterParquet(ShowerAxis const& axis,
- GrammageType dX = 10_g /
- square(1_cm), // profile binning
- unsigned int const nBins = 200); // number of bins
+ LongitudinalProfileWriterParquet(std::array<const char*, NColumns> const& colNames);
/**
* 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.
@@ -69,33 +57,17 @@ namespace corsika {
* This must also increment the run number since we override
* the default behaviour of BaseOutput.
*/
- void endOfLibrary() final override;
-
- /**
- * Add continuous profile.
- */
- template <typename TTrack>
- void write(TTrack const& track, Code const PID, double const weight);
+ void endOfLibrary() override;
/**
- * Add binned profile.
+ * Add profile to disk.
*/
- void write(GrammageType const Xstart, GrammageType const Xend, Code const PID,
- double const weight);
-
- /**
- * Returns a summary of this output.
- */
- YAML::Node getSummary() const;
+ void write(unsigned int const showerId, GrammageType const grammage,
+ std::array<double, NColumns> const& data);
private:
- ParquetStreamer output_; ///< The parquet streamer for this process.
-
- public:
- ShowerAxis const& showerAxis_; ///< conversion between geometry and grammage
- GrammageType dX_; ///< binning of profile.
- unsigned int nBins_; ///< number of profile bins.
- std::vector<ProfileData> profile_; // longitudinal profile
+ std::array<const char*, NColumns> columns_; //!< column names
+ ParquetStreamer output_; //!< The parquet streamer for this process.
}; // namespace corsika
diff --git a/corsika/modules/writers/ParticleWriterParquet.hpp b/corsika/modules/writers/ParticleWriterParquet.hpp
index 6fa0bdd8c5626411eef561561355ed0c58fefaa3..df4ee482b287a875913efe19020cd75e8d315521 100644
--- a/corsika/modules/writers/ParticleWriterParquet.hpp
+++ b/corsika/modules/writers/ParticleWriterParquet.hpp
@@ -69,7 +69,7 @@ namespace corsika {
/**
* If plane is absorbing particles: return the total energy absorbed.
*/
- HEPEnergyType getTotalEnergy() const { return totalEnergy_; }
+ HEPEnergyType getEnergyGround() const { return totalEnergy_; }
private:
ParquetStreamer output_; ///< The primary output file.
diff --git a/examples/cascade_example.cpp b/examples/cascade_example.cpp
index d372773d9e6ce02b434e59c7dc0ec38665c2c6d5..4a33fba5047efc32ad294f501ab6392e35bac382 100644
--- a/examples/cascade_example.cpp
+++ b/examples/cascade_example.cpp
@@ -19,7 +19,6 @@
#include <corsika/output/OutputManager.hpp>
#include <corsika/modules/writers/SubWriter.hpp>
#include <corsika/modules/writers/EnergyLossWriter.hpp>
-#include <corsika/modules/writers/EnergyLossWriterParquet.hpp>
#include <corsika/media/Environment.hpp>
#include <corsika/media/HomogeneousMedium.hpp>
@@ -160,15 +159,15 @@ int main() {
EAS.run();
output.endOfShower();
- const HEPEnergyType Efinal = dEdX.getTotal();
+ const HEPEnergyType Efinal = dEdX.getEnergyLost();
CORSIKA_LOG_INFO(
"\n"
"total cut energy (GeV) : {}\n"
"relative difference (%): {}\n"
"total dEdX energy (GeV): {}\n"
"relative difference (%): {}\n",
- Efinal / 1_GeV, (Efinal / E0 - 1) * 100, dEdX.getTotal() / 1_GeV,
- dEdX.getTotal() / E0 * 100);
+ Efinal / 1_GeV, (Efinal / E0 - 1) * 100, dEdX.getEnergyLost() / 1_GeV,
+ dEdX.getEnergyLost() / E0 * 100);
output.endOfLibrary();
}
diff --git a/examples/cascade_proton_example.cpp b/examples/cascade_proton_example.cpp
index 7739e9e10cb44b4ca17109f1484eb2effd81ed82..3d878acd6c70fcac937e0352927b274198d6c269 100644
--- a/examples/cascade_proton_example.cpp
+++ b/examples/cascade_proton_example.cpp
@@ -19,7 +19,6 @@
#include <corsika/output/OutputManager.hpp>
#include <corsika/modules/writers/SubWriter.hpp>
#include <corsika/modules/writers/EnergyLossWriter.hpp>
-#include <corsika/modules/writers/EnergyLossWriterParquet.hpp>
#include <corsika/media/Environment.hpp>
#include <corsika/media/HomogeneousMedium.hpp>
@@ -129,7 +128,7 @@ int main() {
corsika::pythia8::Decay decay;
ShowerAxis const showerAxis{injectionPos, Vector{rootCS, 0_m, 0_m, -100_km}, env};
- EnergyLossWriter<EnergyLossWriterParquet> dEdX{showerAxis};
+ EnergyLossWriter dEdX{showerAxis};
output.add("energyloss", dEdX);
BetheBlochPDG<SubWriter<decltype(dEdX)>> eLoss{dEdX};
diff --git a/examples/corsika.cpp b/examples/corsika.cpp
index ffae78868fad8859d63a33cd3a250ea8d3e46329..392f7cd22eaed92a536ac4a271fa1d54491a2c7f 100644
--- a/examples/corsika.cpp
+++ b/examples/corsika.cpp
@@ -29,7 +29,7 @@
#include <corsika/output/OutputManager.hpp>
#include <corsika/modules/writers/SubWriter.hpp>
#include <corsika/modules/writers/EnergyLossWriter.hpp>
-#include <corsika/modules/writers/EnergyLossWriterParquet.hpp>
+#include <corsika/modules/writers/LongitudinalWriter.hpp>
#include <corsika/media/Environment.hpp>
#include <corsika/media/FlatExponential.hpp>
@@ -64,7 +64,6 @@
#include <CLI/Config.hpp>
#include <iomanip>
-#include <iostream>
#include <limits>
#include <string>
@@ -187,8 +186,7 @@ int main(int argc, char** argv) {
// gets all messed up
if (app.count("--pdg") == 0) {
if ((app.count("-A") == 0) || (app.count("-Z") == 0)) {
- std::cerr << "If --pdg is not provided, then both -A and -Z are required."
- << std::endl;
+ CORSIKA_LOG_ERROR("If --pdg is not provided, then both -A and -Z are required.");
return 1;
}
}
@@ -246,7 +244,7 @@ int main(int argc, char** argv) {
auto const phiRad = app["--azimuth"]->as<double>() / 180. * M_PI;
// convert Elab to Plab
- HEPMomentumType P0 = sqrt((E0 - mass) * (E0 + mass));
+ HEPMomentumType P0 = calculate_momentum(E0, mass);
// convert the momentum to the zenith and azimuth angle of the primary
auto const [px, py, pz] =
@@ -321,13 +319,13 @@ int main(int argc, char** argv) {
corsika::proposal::Interaction emCascade(env);
// NOT available for PROPOSAL due to interface trouble:
- // InteractionCounter emCascadeCounted(emCascade);
+ // InteractionCounter emCascadeCounted(emCascade);
// corsika::proposal::ContinuousProcess<SubWriter<decltype(dEdX)>> emContinuous(env);
BetheBlochPDG<SubWriter<decltype(dEdX)>> emContinuous{dEdX};
- LongitudinalProfile<corsika::LongitudinalProfileWriterParquet> longprof{
- showerAxis, 10_g / square(1_cm), 200};
- output.add("profile", longprof);
+ LongitudinalWriter profile{showerAxis, 10_g / square(1_cm), 200};
+ output.add("profile", profile);
+ LongitudinalProfile<SubWriter<decltype(profile)>> longprof{profile};
corsika::urqmd::UrQMD urqmd;
InteractionCounter urqmdCounted(urqmd);
@@ -389,7 +387,6 @@ int main(int argc, char** argv) {
string const outdir(app["--filename"]->as<std::string>());
string const labHist_file = outdir + "/inthist_lab_" + to_string(i_shower) + ".npz";
string const cMSHist_file = outdir + "/inthist_cms_" + to_string(i_shower) + ".npz";
- string const longprof_file = outdir + "/longprof_" + to_string(i_shower) + ".txt";
// setup particle stack, and add primary particle
stack.clear();
@@ -408,10 +405,14 @@ int main(int argc, char** argv) {
// run the shower
EAS.run();
- HEPEnergyType const Efinal = dEdX.getTotal();
- // +observationLevel.getTotalEnergy();
- cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl
- << "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
+ HEPEnergyType const Efinal =
+ dEdX.getEnergyLost() + observationLevel.getEnergyGround();
+
+ CORSIKA_LOG_INFO(
+ "total energy budget (GeV): {} (dEdX={} ground={}), "
+ "relative difference (%): {}",
+ Efinal / 1_GeV, dEdX.getEnergyLost() / 1_GeV,
+ observationLevel.getEnergyGround() / 1_GeV, (Efinal / E0 - 1) * 100);
// auto const hists = heModelCounted.getHistogram() + urqmdCounted.getHistogram();
auto const hists = sibyllCounted.getHistogram() + sibyllNucCounted.getHistogram() +
diff --git a/examples/em_shower.cpp b/examples/em_shower.cpp
index 8aa6d81daa37bb35889e035e5ec54dd2af6a5b60..ff94eae5cbb12a308f0e5f4352c558b18d3f596e 100644
--- a/examples/em_shower.cpp
+++ b/examples/em_shower.cpp
@@ -24,7 +24,7 @@
#include <corsika/output/OutputManager.hpp>
#include <corsika/modules/writers/SubWriter.hpp>
#include <corsika/modules/writers/EnergyLossWriter.hpp>
-#include <corsika/modules/writers/EnergyLossWriterParquet.hpp>
+#include <corsika/modules/writers/LongitudinalWriter.hpp>
#include <corsika/media/Environment.hpp>
#include <corsika/media/LayeredSphericalAtmosphereBuilder.hpp>
@@ -163,10 +163,10 @@ int main(int argc, char** argv) {
TrackWriter tracks;
output.add("tracks", tracks);
- // long. profile; columns for photon, e+, e- still need to be added
- LongitudinalProfile<corsika::LongitudinalProfileWriterParquet> longprof{
- showerAxis, 10_g / square(1_cm), 200};
- output.add("profile", longprof);
+ // long. profile
+ LongitudinalWriter profile{showerAxis, 10_g / square(1_cm), 200};
+ output.add("profile", profile);
+ LongitudinalProfile<SubWriter<decltype(profile)>> longprof{profile};
Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.}));
ObservationPlane<setup::Tracking, ParticleWriterParquet> observationLevel{
@@ -187,9 +187,12 @@ int main(int argc, char** argv) {
EAS.run();
output.endOfShower();
- const HEPEnergyType Efinal = dEdX.getTotal();
- cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl
- << "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
+ HEPEnergyType const Efinal = dEdX.getEnergyLost() + observationLevel.getEnergyGround();
+
+ CORSIKA_LOG_INFO(
+ "total energy budget (GeV): {}, "
+ "relative difference (%): {}",
+ Efinal / 1_GeV, (Efinal / E0 - 1) * 100);
output.endOfLibrary();
}
diff --git a/examples/hybrid_MC.cpp b/examples/hybrid_MC.cpp
index 904c0d082c2a12efe17a30197bdb136a3f732d43..fa635fe380582d43ebd304b93b76e4178b8c2452 100644
--- a/examples/hybrid_MC.cpp
+++ b/examples/hybrid_MC.cpp
@@ -29,7 +29,7 @@
#include <corsika/output/OutputManager.hpp>
#include <corsika/modules/writers/SubWriter.hpp>
#include <corsika/modules/writers/EnergyLossWriter.hpp>
-#include <corsika/modules/writers/EnergyLossWriterParquet.hpp>
+#include <corsika/modules/writers/LongitudinalParquet.hpp>
#include <corsika/media/Environment.hpp>
#include <corsika/media/FlatExponential.hpp>
@@ -242,11 +242,9 @@ int main(int argc, char** argv) {
CONEXhybrid conex_model(center, showerAxis, t, injectionHeight, E0,
get_PDG(Code::Proton));
- LongitudinalProfile<corsika::LongitudinalProfileWriterParquet> longprof{
- showerAxis, 10_g / square(1_cm), 200};
- output.add("profile", longprof);
-
- LongitudinalProfile longprof(showerAxis);
+ LongitudinalWriter profile{showerAxis, 10_g / square(1_cm), 200};
+ output.add("profile", profile);
+ LongitudinalProfile<SubWriter<decltype(profile)>> longprof{profile};
Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.}));
ObservationPlane<setup::Tracking> observationLevel(
diff --git a/examples/mars.cpp b/examples/mars.cpp
index bc7fa56e03fbadd0aef5b4639330cfbec2e57a0c..0b64ab765bf2056b12f5b753dccf6e4f1d0275ee 100644
--- a/examples/mars.cpp
+++ b/examples/mars.cpp
@@ -31,7 +31,7 @@
#include <corsika/output/OutputManager.hpp>
#include <corsika/modules/writers/SubWriter.hpp>
#include <corsika/modules/writers/EnergyLossWriter.hpp>
-#include <corsika/modules/writers/EnergyLossWriterParquet.hpp>
+#include <corsika/modules/writers/LongitudinalWriter.hpp>
#include <corsika/media/Environment.hpp>
#include <corsika/media/FlatExponential.hpp>
@@ -47,7 +47,6 @@
#include <corsika/modules/BetheBlochPDG.hpp>
#include <corsika/modules/LongitudinalProfile.hpp>
#include <corsika/modules/ObservationPlane.hpp>
-#include <corsika/modules/OnShellCheck.hpp>
#include <corsika/modules/StackInspector.hpp>
#include <corsika/modules/TrackWriter.hpp>
#include <corsika/modules/ParticleCut.hpp>
@@ -210,8 +209,7 @@ int main(int argc, char** argv) {
// gets all messed up
if (app.count("--pdg") == 0) {
if ((app.count("-A") == 0) || (app.count("-Z") == 0)) {
- std::cerr << "If --pdg is not provided, then both -A and -Z are required."
- << std::endl;
+ CORSIKA_LOG_ERROR("If --pdg is not provided, then both -A and -Z are required.");
return 1;
}
}
@@ -314,7 +312,7 @@ int main(int argc, char** argv) {
ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.2, env};
- EnergyLossWriter<EnergyLossWriterParquet> dEdX{showerAxis};
+ EnergyLossWriter dEdX{showerAxis};
output.add("energyloss", dEdX);
HEPEnergyType const emcut = 1_GeV;
@@ -359,8 +357,9 @@ int main(int argc, char** argv) {
// corsika::proposal::ContinuousProcess emContinuous(env);
BetheBlochPDG<SubWriter<decltype(dEdX)>> emContinuous{dEdX};
- LongitudinalProfile<corsika::LongitudinalProfileWriterParquet> profile{showerAxis};
- output.add("profile", profile);
+ LongitudinalWriter longprof{showerAxis};
+ output.add("profile", longprof);
+ LongitudinalProfile<SubWriter<decltype(longprof)>> profile{longprof};
corsika::urqmd::UrQMD urqmd;
InteractionCounter urqmdCounted{urqmd};
@@ -435,10 +434,13 @@ int main(int argc, char** argv) {
// run the shower
EAS.run();
- const HEPEnergyType Efinal = dEdX.getTotal();
- // +observationLevel.getEnergyGround();
- cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl
- << "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
+ HEPEnergyType const Efinal =
+ dEdX.getEnergyLost() + observationLevel.getEnergyGround();
+
+ CORSIKA_LOG_INFO(
+ "total energy budget (GeV): {}, "
+ "relative difference (%): {}",
+ Efinal / 1_GeV, (Efinal / E0 - 1) * 100);
auto const hists = sibyllCounted.getHistogram() + sibyllNucCounted.getHistogram() +
urqmdCounted.getHistogram();
diff --git a/examples/vertical_EAS.cpp b/examples/vertical_EAS.cpp
index fabcc0763f8ec290860b6f2330eda4d37adbcea1..ba69ffe2814e599a3f7bf6eb582a864beffc62c5 100644
--- a/examples/vertical_EAS.cpp
+++ b/examples/vertical_EAS.cpp
@@ -31,7 +31,7 @@
#include <corsika/output/OutputManager.hpp>
#include <corsika/modules/writers/SubWriter.hpp>
#include <corsika/modules/writers/EnergyLossWriter.hpp>
-#include <corsika/modules/writers/EnergyLossWriterParquet.hpp>
+#include <corsika/modules/writers/LongitudinalWriter.hpp>
#include <corsika/media/Environment.hpp>
#include <corsika/media/FlatExponential.hpp>
@@ -198,7 +198,7 @@ int main(int argc, char** argv) {
// register profile output
// construct the overall energy loss writer and register it
- EnergyLossWriter<EnergyLossWriterParquet> dEdX{showerAxis};
+ EnergyLossWriter dEdX{showerAxis};
output.add("energyloss", dEdX);
// construct the continuous energy loss model
@@ -208,8 +208,10 @@ int main(int argc, char** argv) {
ParticleCut<SubWriter<decltype(dEdX)>> cut{E0, E0, 60_GeV, 60_GeV, true, dEdX};
// setup longitudinal profile
- LongitudinalProfile<corsika::LongitudinalProfileWriterParquet> profile{showerAxis};
- output.add("profile", profile);
+ LongitudinalWriter longProf{showerAxis};
+ output.add("profile", longProf);
+
+ LongitudinalProfile<SubWriter<decltype(longProf)>> profile{longProf};
// create a track writer and register it with the output manager
TrackWriter<TrackWriterParquet> trackWriter;