more sophisticated histogramming

Stack/History/HistoryObservationPlane.cpp

@@ -11,7 +11,8 @@
 
 #include <boost/histogram/ostream.hpp>
 
-#include <fstream>
+#include <iomanip>
+#include <iostream>
 
 using namespace corsika::units::si;
 using namespace corsika::history;
@@ -65,17 +66,20 @@ LengthType HistoryObservationPlane::MaxStepLength(setup::Stack::ParticleType con
 
 void HistoryObservationPlane::fillHistoryHistogram(
     setup::Stack::ParticleType const& muon) {
-  // double const muonEnergy = muon.GetEnergy() / 1_eV;
-  // auto parent = stack_.at(muon.GetEvent()->projectileIndex());
-  Event const* event = muon.GetEvent().get();
+  double const muon_energy = muon.GetEnergy() / 1_GeV;
 
-  int intCounter = 0;
+  int genctr{0};
+  Event const* event = muon.GetEvent().get();
   while (event) {
-    if (event->eventType() == EventType::Interaction) intCounter++;
-    event = event->parentEvent().get();
+    auto const projectile = stack_.cfirst() + event->projectileIndex();
+    if (event->eventType() == EventType::Interaction) {
+      genctr++;
+      double const projEnergy = projectile.GetEnergy() / 1_GeV;
+      int const pdg = static_cast<int>(particles::GetPDG(projectile.GetPID()));
+
+      histogram_(muon_energy, projEnergy, pdg);
+    }
+    event = event->parentEvent().get(); // projectile.GetEvent().get();
   }
-
-  histogram_(intCounter);
 }
 
-void HistoryObservationPlane::print() { C8LOG_INFO(histogram_); }

Stack/History/HistoryObservationPlane.hpp

@@ -23,9 +23,12 @@ namespace corsika::history {
     inline auto hist_factory() {
       namespace bh = boost::histogram;
       namespace bha = bh::axis;
-      auto h =
-          bh::make_histogram(bha::integer<int, bh::use_default, bha::option::growth_t>{
-              0, 10, "hadronic generation"});
+      auto h = bh::make_histogram(
+          bha::regular<float, bha::transform::log>{11 * 5, 1e0, 1e11, "muon energy"},
+          bha::regular<float, bha::transform::log>{11 * 5, 1e0, 1e11,
+                                                   "projectile energy"},
+          bha::category<int, bh::use_default, bha::option::growth_t>{
+              {211, -211, 2212, -2212}, "projectile PDG"});
       return h;
     }
   } // namespace detail
@@ -43,7 +46,7 @@ namespace corsika::history {
         corsika::setup::Stack::ParticleType const& vParticle,
         corsika::setup::Trajectory const& vTrajectory);
 
-    void print();
+    auto const& histogram() const { return histogram_; }
 
   private:
     void fillHistoryHistogram(setup::Stack::ParticleType const&);