Merge branch '651-particle-writer-reports-energy-per-particle-type-in-summary-file' into 'master'

Resolve "particle writer reports energy per particle-type in summary file"

Closes #651

See merge request !582

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

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
 

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); }