diff --git a/Stack/History/Event.hpp b/Stack/History/Event.hpp index e47eaf2eca14f81c9f91a486e62516d3911f7dd1..c89f8064eb8e4a527384ba87c8aee9c3601546f5 100644 --- a/Stack/History/Event.hpp +++ b/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_; } diff --git a/Stack/History/HistoryObservationPlane.cc b/Stack/History/HistoryObservationPlane.cc index ebc8ab11d989c74975c012bb32f75cdb6a442a6e..f65ce899377193f9468ba34098ec7efa4a82b7fd 100644 --- a/Stack/History/HistoryObservationPlane.cc +++ b/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); } diff --git a/Stack/History/HistoryObservationPlane.hpp b/Stack/History/HistoryObservationPlane.hpp index 8c2d5700daa63cce61d9f88705047f6aa2cf663e..a603963dba8c6ee746a0ddff1630f1cae6ff8cc9 100644 --- a/Stack/History/HistoryObservationPlane.hpp +++ b/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_;