From 0cf01bf1d44481f89b02fd53fd829096b60ae123 Mon Sep 17 00:00:00 2001
From: Maximilian Reininghaus <maximilian.reininghaus@kit.edu>
Date: Tue, 6 Oct 2020 09:50:13 +0200
Subject: [PATCH] more sophisticated histogramming

---
 Stack/History/HistoryObservationPlane.cpp | 24 +++++++++++++----------
 Stack/History/HistoryObservationPlane.hpp | 11 +++++++----
 2 files changed, 21 insertions(+), 14 deletions(-)

diff --git a/Stack/History/HistoryObservationPlane.cpp b/Stack/History/HistoryObservationPlane.cpp
index eb5ca8b01..1109b0fce 100644
--- a/Stack/History/HistoryObservationPlane.cpp
+++ b/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_); }
diff --git a/Stack/History/HistoryObservationPlane.hpp b/Stack/History/HistoryObservationPlane.hpp
index 11c4b86cc..fb67ca7f9 100644
--- a/Stack/History/HistoryObservationPlane.hpp
+++ b/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&);
-- 
GitLab