diff --git a/corsika/detail/modules/writers/ParticleWriterParquet.inl b/corsika/detail/modules/writers/ParticleWriterParquet.inl
index d4780a71380b47c3be1b387b0c9639e1b452bb14..640c1c704b783cd4bd2bcbaab2391abd5817d931 100644
--- a/corsika/detail/modules/writers/ParticleWriterParquet.inl
+++ b/corsika/detail/modules/writers/ParticleWriterParquet.inl
@@ -15,7 +15,6 @@ namespace corsika {
inline ParticleWriterParquet::ParticleWriterParquet(bool const printZ)
: output_()
, showerId_(0)
- , totalEnergy_(0_eV)
, printZ_(printZ) {}
inline void ParticleWriterParquet::startOfLibrary(
@@ -59,18 +58,57 @@ namespace corsika {
inline void ParticleWriterParquet::startOfShower(unsigned int const showerId) {
showerId_ = showerId;
- totalEnergy_ = 0_eV;
+
countHadrons_ = 0;
countOthers_ = 0;
countEM_ = 0;
countMuons_ = 0;
+
+ kineticEnergyHadrons_ = 0_eV;
+ kineticEnergyMuons_ = 0_eV;
+ kineticEnergyEM_ = 0_eV;
+ kineticEnergyOthers_ = 0_eV;
+
+ totalEnergyHadrons_ = 0_eV;
+ totalEnergyMuons_ = 0_eV;
+ totalEnergyEM_ = 0_eV;
+ totalEnergyOthers_ = 0_eV;
}
- inline void ParticleWriterParquet::endOfShower(unsigned int const) {}
+ inline void ParticleWriterParquet::endOfShower(unsigned int const) {
+ summary_["shower_" + std::to_string(showerId_)]["hadron"]["count"] = countHadrons_;
+ summary_["shower_" + std::to_string(showerId_)]["hadron"]["kinetic_energy"] =
+ kineticEnergyHadrons_ / 1_GeV;
+ summary_["shower_" + std::to_string(showerId_)]["hadron"]["total_energy"] =
+ totalEnergyHadrons_ / 1_GeV;
+
+ summary_["shower_" + std::to_string(showerId_)]["muon"]["count"] = countMuons_;
+ summary_["shower_" + std::to_string(showerId_)]["muon"]["kinetic_energy"] =
+ kineticEnergyMuons_ / 1_GeV;
+ summary_["shower_" + std::to_string(showerId_)]["muon"]["total_energy"] =
+ totalEnergyMuons_ / 1_GeV;
+
+ summary_["shower_" + std::to_string(showerId_)]["em"]["count"] = countEM_;
+ summary_["shower_" + std::to_string(showerId_)]["em"]["kinetic_energy"] =
+ kineticEnergyEM_ / 1_GeV;
+ summary_["shower_" + std::to_string(showerId_)]["em"]["total_energy"] =
+ totalEnergyEM_ / 1_GeV;
+
+ summary_["shower_" + std::to_string(showerId_)]["other"]["count"] = countOthers_;
+ summary_["shower_" + std::to_string(showerId_)]["other"]["kinetic_energy"] =
+ kineticEnergyOthers_ / 1_GeV;
+ summary_["shower_" + std::to_string(showerId_)]["other"]["total_energy"] =
+ totalEnergyOthers_ / 1_GeV;
+ }
inline void ParticleWriterParquet::endOfLibrary() { output_.closeStreamer(); }
- inline void ParticleWriterParquet::write(Code const pid, HEPEnergyType const energy,
+ inline HEPEnergyType ParticleWriterParquet::getEnergyGround() const {
+ return totalEnergyHadrons_ + totalEnergyMuons_ + totalEnergyEM_ + totalEnergyOthers_;
+ }
+
+ inline void ParticleWriterParquet::write(Code const pid,
+ HEPEnergyType const kineticEnergy,
LengthType const x, LengthType const y,
LengthType const z, double const nx,
double const ny, double const nz,
@@ -78,7 +116,7 @@ namespace corsika {
// write the next row - we must write `shower_` first.
*(output_.getWriter()) << showerId_ << static_cast<int>(get_PDG(pid))
- << static_cast<float>(energy / 1_GeV)
+ << static_cast<float>(kineticEnergy / 1_GeV)
<< static_cast<float>(x / 1_m) << static_cast<float>(y / 1_m);
if (printZ_) { *(output_.getWriter()) << static_cast<float>(z / 1_m); }
@@ -86,30 +124,28 @@ namespace corsika {
<< static_cast<float>(nz) << static_cast<double>(t / 1_s)
<< static_cast<float>(weight) << parquet::EndRow;
- totalEnergy_ += energy;
-
if (is_hadron(pid)) {
++countHadrons_;
+ kineticEnergyHadrons_ += kineticEnergy;
+ totalEnergyHadrons_ += kineticEnergy + get_mass(pid);
} else if (is_muon(pid)) {
++countMuons_;
+ kineticEnergyMuons_ += kineticEnergy;
+ totalEnergyMuons_ += kineticEnergy + get_mass(pid);
} else if (is_em(pid)) {
++countEM_;
+ kineticEnergyEM_ += kineticEnergy;
+ totalEnergyEM_ += kineticEnergy + get_mass(pid);
} else {
++countOthers_;
+ kineticEnergyOthers_ += kineticEnergy;
+ totalEnergyOthers_ += kineticEnergy + get_mass(pid);
}
}
/**
* Return collected library-level summary for output.
*/
- inline YAML::Node ParticleWriterParquet::getSummary() const {
- YAML::Node summary;
- summary["Eground"] = totalEnergy_ / 1_GeV;
- summary["hadrons"] = countHadrons_;
- summary["muons"] = countMuons_;
- summary["em"] = countEM_;
- summary["others"] = countOthers_;
- return summary;
- }
+ inline YAML::Node ParticleWriterParquet::getSummary() const { return summary_; }
} // namespace corsika
diff --git a/corsika/modules/writers/ParticleWriterParquet.hpp b/corsika/modules/writers/ParticleWriterParquet.hpp
index 4e15a34d6420d20fd1834139e68e74f52fde6f22..4eded3f6f5b1b88a405a419f5f0e279854950e1c 100644
--- a/corsika/modules/writers/ParticleWriterParquet.hpp
+++ b/corsika/modules/writers/ParticleWriterParquet.hpp
@@ -49,7 +49,7 @@ namespace corsika {
/**
* Write a PDG/corsika::Code particle to the file.
*/
- void write(Code const pid, units::si::HEPEnergyType const energy,
+ void write(Code const pid, units::si::HEPEnergyType const kineticEnergy,
units::si::LengthType const x, units::si::LengthType const y,
units::si::LengthType const z, double const nx, double const ny,
double const nz, units::si::TimeType const time, const double weight);
@@ -62,7 +62,7 @@ namespace corsika {
/**
* If plane is absorbing particles: return the total energy absorbed.
*/
- HEPEnergyType getEnergyGround() const { return totalEnergy_; }
+ HEPEnergyType getEnergyGround() const;
private:
ParquetStreamer output_; ///< The primary output file.
@@ -73,7 +73,21 @@ namespace corsika {
double countEM_ = 0; ///< count EM particles hitting plane.
double countOthers_ = 0; ///< count other types of particles hitting plane
- HEPEnergyType totalEnergy_; ///< energy absorbed in ground.
+ HEPEnergyType kineticEnergyHadrons_ =
+ 0_eV; ///< kinetic energy of hadrons hitting plane
+ HEPEnergyType kineticEnergyMuons_ = 0_eV; ///< kinetic energy of muons hitting plane
+ HEPEnergyType kineticEnergyEM_ =
+ 0_eV; ///< kinetic energy of EM particles hitting plane.
+ HEPEnergyType kineticEnergyOthers_ =
+ 0_eV; ///< kinetic energy of other types of particles hitting plane
+
+ HEPEnergyType totalEnergyHadrons_ = 0_eV; ///< total energy of hadrons hitting plane
+ HEPEnergyType totalEnergyMuons_ = 0_eV; ///< total energy of muons hitting plane
+ HEPEnergyType totalEnergyEM_ = 0_eV; ///< total energy of EM particles hitting plane.
+ HEPEnergyType totalEnergyOthers_ =
+ 0_eV; ///< total energy of other types of particles hitting plane
+
+ YAML::Node summary_;
bool const printZ_; ///< flag to print the z coordinate
diff --git a/tests/output/testWriterObservationPlane.cpp b/tests/output/testWriterObservationPlane.cpp
index 566b7b6c0410e6f4046d447215cb486f3bb623db..ff054358388245d8159ad01c88ac3c122a1c769c 100644
--- a/tests/output/testWriterObservationPlane.cpp
+++ b/tests/output/testWriterObservationPlane.cpp
@@ -48,23 +48,26 @@ TEST_CASE("ObservationPlaneWriterParquet") {
TestWriterPlane test(true);
test.startOfLibrary(file_dir);
- test.startOfShower(0);
- // write a few particles
+ test.startOfShower(0);
test.checkWrite();
-
test.endOfShower(0);
+
test.endOfLibrary();
CHECK(boost::filesystem::exists(file_dir + "/particles.parquet"));
auto const summary = test.getSummary();
- CHECK(summary["Eground"].as<double>() == Approx(5));
- CHECK(summary["hadrons"].as<int>() == Approx(1));
- CHECK(summary["muons"].as<int>() == Approx(2));
- CHECK(summary["em"].as<int>() == Approx(1));
- CHECK(summary["others"].as<int>() == Approx(1));
+ CHECK(summary["shower_0"]["hadron"]["count"].as<double>() == Approx(1));
+ CHECK(summary["shower_0"]["muon"]["count"].as<double>() == Approx(2));
+ CHECK(summary["shower_0"]["em"]["count"].as<double>() == Approx(1));
+ CHECK(summary["shower_0"]["other"]["count"].as<double>() == Approx(1));
+
+ CHECK(summary["shower_0"]["hadron"]["kinetic_energy"].as<double>() == Approx(1));
+ CHECK(summary["shower_0"]["muon"]["kinetic_energy"].as<double>() == Approx(2));
+ CHECK(summary["shower_0"]["em"]["kinetic_energy"].as<double>() == Approx(1));
+ CHECK(summary["shower_0"]["other"]["kinetic_energy"].as<double>() == Approx(1));
// clean things up
if (boost::filesystem::exists(file_dir)) { boost::filesystem::remove_all(file_dir); }
@@ -80,23 +83,65 @@ TEST_CASE("ObservationPlaneWriterParquet") {
TestWriterPlane test(false); // do not print the z-coord
test.startOfLibrary(file_dir);
- test.startOfShower(0);
- // write a few particles
+ test.startOfShower(0);
test.checkWrite();
-
test.endOfShower(0);
+
test.endOfLibrary();
CHECK(boost::filesystem::exists(file_dir + "/particles.parquet"));
auto const summary = test.getSummary();
- CHECK(summary["Eground"].as<double>() == Approx(5));
- CHECK(summary["hadrons"].as<int>() == Approx(1));
- CHECK(summary["muons"].as<int>() == Approx(2));
- CHECK(summary["em"].as<int>() == Approx(1));
- CHECK(summary["others"].as<int>() == Approx(1));
+ CHECK(summary["shower_0"]["hadron"]["count"].as<double>() == Approx(1));
+ CHECK(summary["shower_0"]["muon"]["count"].as<double>() == Approx(2));
+ CHECK(summary["shower_0"]["em"]["count"].as<double>() == Approx(1));
+ CHECK(summary["shower_0"]["other"]["count"].as<double>() == Approx(1));
+
+ CHECK(summary["shower_0"]["hadron"]["kinetic_energy"].as<double>() == Approx(1));
+ CHECK(summary["shower_0"]["muon"]["kinetic_energy"].as<double>() == Approx(2));
+ CHECK(summary["shower_0"]["em"]["kinetic_energy"].as<double>() == Approx(1));
+ CHECK(summary["shower_0"]["other"]["kinetic_energy"].as<double>() == Approx(1));
+
+ // clean things up
+ if (boost::filesystem::exists(file_dir)) { boost::filesystem::remove_all(file_dir); }
+ }
+
+ SECTION("reset_between_showers") {
+ // Check that the counters are reset between showers
+
+ std::string const file_dir = "./output_dir_obs_plane_reset";
+
+ // preparation
+ if (boost::filesystem::exists(file_dir)) { boost::filesystem::remove_all(file_dir); }
+ boost::filesystem::create_directory(file_dir);
+
+ TestWriterPlane test(false); // do not print the z-coord
+ test.startOfLibrary(file_dir);
+
+ test.startOfShower(0);
+ test.checkWrite();
+ test.endOfShower(0);
+
+ test.startOfShower(1);
+ test.checkWrite();
+ test.endOfShower(1);
+
+ auto const summary = test.getSummary();
+
+ CHECK(summary["shower_0"]["hadron"]["count"].as<double>() ==
+ summary["shower_1"]["hadron"]["count"].as<double>());
+ CHECK(summary["shower_0"]["muon"]["count"].as<double>() ==
+ summary["shower_1"]["muon"]["count"].as<double>());
+ CHECK(summary["shower_0"]["em"]["count"].as<double>() ==
+ summary["shower_1"]["em"]["count"].as<double>());
+ CHECK(summary["shower_0"]["other"]["count"].as<double>() ==
+ summary["shower_1"]["other"]["count"].as<double>());
+
+ CHECK(boost::filesystem::exists(file_dir + "/particles.parquet"));
+
+ test.endOfLibrary();
// clean things up
if (boost::filesystem::exists(file_dir)) { boost::filesystem::remove_all(file_dir); }