simple HistoryObservationPlane

Stack/History/Event.hpp

@@ -42,6 +42,9 @@ namespace corsika::history {
 
     template <typename TStackIterator>
     TStackIterator projectile(TStackIterator begin) {
+      // MR: This is dangerous. You can pass any iterator though it must
+      // be stack.begin() to yield the correct projectile
+
       return begin + projectileIndex_;
     }
 

Stack/History/HistoryObservationPlane.cc

@@ -14,10 +14,12 @@ using namespace corsika::units::si;
 using namespace corsika::history;
 using namespace corsika;
 
-HistoryObservationPlane::HistoryObservationPlane(geometry::Plane const& obsPlane,
+HistoryObservationPlane::HistoryObservationPlane(setup::Stack const& stack,
+                                                 geometry::Plane const& obsPlane,
                                                  bool deleteOnHit)
-    : plane_(obsPlane)
-    , deleteOnHit_(deleteOnHit) {}
+    : stack_{stack}
+    , plane_{obsPlane}
+    , deleteOnHit_{deleteOnHit} {}
 
 corsika::process::EProcessReturn HistoryObservationPlane::DoContinuous(
     setup::Stack::ParticleType const& particle, setup::Trajectory const& trajectory) {
@@ -55,7 +57,18 @@ LengthType HistoryObservationPlane::MaxStepLength(setup::Stack::ParticleType con
   return (trajectory.GetR0() - pointOfIntersection).norm() * 1.0001;
 }
 
-void fillHistoryHistogram(setup::Stack::ParticleType const& muon) {
+void HistoryObservationPlane::fillHistoryHistogram(
+    setup::Stack::ParticleType const& muon) {
   double const muonEnergy = muon.GetEnergy() / 1_eV;
 
+  auto parent = stack_.begin() + muon.GetEvent()->projectileIndex();
+
+  auto* event = muon.GetEvent().get();
+
+  int intCounter = 0;
+  while (event) {
+    event = event->parentEvent().get();
+    intCounter++;
+  }
+  histogram_(intCounter);
 }

Stack/History/HistoryObservationPlane.hpp

@@ -21,7 +21,7 @@
 namespace corsika::history {
   namespace detail {
     auto hist_factory() {
-      auto h = boost::histogram::make_histogram(
+      /*auto h = boost::histogram::make_histogram(
           boost::histogram::axis::regular<double, boost::histogram::axis::transform::log>{
               130, 1e8, 1e21, "muon energy/eV"},
           boost::histogram::axis::integer<int, boost::histogram::use_default,
@@ -30,7 +30,12 @@ namespace corsika::history {
           boost::histogram::axis::regular<double, boost::histogram::axis::transform::log>{
               130, 1e8, 1e21, "hadronic energy/eV"},
           boost::histogram::axis::category<int, boost::histogram::use_default,
-                                           boost::histogram::axis::option::growth_t>{});
+                                           boost::histogram::axis::option::growth_t>{});*/
+
+      auto h = boost::histogram::make_histogram(
+          boost::histogram::axis::integer<int, boost::histogram::use_default,
+                                          boost::histogram::axis::option::growth_t>{
+              0, 10, "hadronic generation"});
       return h;
     }
   } // namespace detail
@@ -38,9 +43,9 @@ namespace corsika::history {
   class HistoryObservationPlane
       : public corsika::process::ContinuousProcess<HistoryObservationPlane> {
   public:
-    HistoryObservationPlane(geometry::Plane const&, bool = true);
+    HistoryObservationPlane(setup::Stack const&, geometry::Plane const&, bool = true);
 
-    void save(std::string const&);
+    //~ void save(std::string const&);
 
     corsika::units::si::LengthType MaxStepLength(
         corsika::setup::Stack::ParticleType const&,
@@ -53,6 +58,7 @@ namespace corsika::history {
   private:
     void fillHistoryHistogram(setup::Stack::ParticleType const&);
 
+    setup::Stack const& stack_;
     geometry::Plane const plane_;
     bool const deleteOnHit_;