started HistoryObservationPlane

Documentation/Examples/CMakeLists.txt

@@ -91,6 +91,11 @@ if (Pythia8_FOUND)
     CORSIKAunits
     CORSIKAlogging
     CORSIKArandom
+    CORSIKAhistory
+    ProcessSibyll
+    ProcessPythia8
+    ProcessUrQMD
+    ProcessSwitch
     CORSIKAcascade
     ProcessProposal
     ProcessPythia8

Stack/History/CMakeLists.txt

@@ -2,6 +2,7 @@ set (
   HISTORY_HEADERS
   Event.hpp
   HistorySecondaryView.hpp
+  HistoryObservationPlane.hpp
   HistoryStackExtension.h
   SecondaryParticle.hpp
   )
@@ -11,16 +12,23 @@ set (
   corsika/history
   )
 
-add_library (CORSIKAhistory INTERFACE)
+set (
+  HISTORY_SOURCES
+  HistoryObservationPlane.cc
+)
+
+#add_library (CORSIKAhistory INTERFACE)
+add_library (CORSIKAhistory STATIC ${HISTORY_SOURCES})
 CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKAhistory ${HISTORY_NAMESPACE} ${HISTORY_HEADERS})
 
 target_link_libraries (
   CORSIKAhistory
-  INTERFACE
+#  INTERFACE
   CORSIKAparticles
   CORSIKAgeometry
   SuperStupidStack
   NuclearStackExtension
+  C8::ext::boost
   )
 
 target_include_directories (

Stack/History/Event.hpp

@@ -12,52 +12,56 @@
 #include <corsika/history/SecondaryParticle.hpp>
 
 #include <iostream>
+#include <memory>
 #include <optional>
 #include <vector>
-#include <memory>
 
 namespace corsika::history {
 
   class Event;
   using EvtPtr = std::shared_ptr<history::Event>;
-  
+
   class Event {
-    
+
     size_t projectileIndex_ = 0; //!< reference to projectile (for secondaries_)
     std::vector<SecondaryParticle> secondaries_;
     EvtPtr parent_event_;
-    size_t sec_index_ = 0;
-    
+    size_t sec_index_ =
+        0; //!< index of the projectile of this Event in the parent_event_.secondaries_;
+
     std::optional<corsika::particles::Code>
         targetCode_; // cannot be const, value set only after construction
 
   public:
     Event() = default;
 
+    void setParentEvent(EvtPtr& evt) { parent_event_ = evt; }
+
     void setParentEventAndSecondaryIndex(EvtPtr& evt, size_t sec_index) {
-      parent_event_ = evt;
-      sec_index_ = sec_index; // this is the index of the projectile
-			      // of this Event in the parent_event_;
+      setParentEvent(evt);
+      setParentSecondaryIndex(sec_index);
     }
 
+    void setParentSecondaryIndex(size_t sec_index) { sec_index_ = sec_index; }
+
     EvtPtr parentEvent() { return parent_event_; }
 
-    void setProjectileIndex(size_t i) { projectileIndex_=i; }
+    void setProjectileIndex(size_t i) { projectileIndex_ = i; }
     size_t projectileIndex() const { return projectileIndex_; }
-    
-    template<typename TStackIterator>
-    TStackIterator projectile(TStackIterator begin) { return begin+projectileIndex_; } 
-    
+
+    template <typename TStackIterator>
+    TStackIterator projectile(TStackIterator begin) {
+      return begin + projectileIndex_;
+    }
+
     void addSecondary(units::si::HEPEnergyType energy,
                       geometry::Vector<units::si::hepmomentum_d> momentum,
                       particles::Code pid) {
       secondaries_.emplace_back(energy, momentum, pid);
     }
     std::vector<SecondaryParticle>& secondaries() { return secondaries_; }
-    
-    void setTargetCode(const particles::Code t) { targetCode_=t; }
 
-    
+    void setTargetCode(const particles::Code t) { targetCode_ = t; }
   };
-  
+
 } // namespace corsika::history

Stack/History/HistoryObservationPlane.cc

@@ -12,7 +12,6 @@
 #include <boost/histogram/ostream.hpp>
 
 #include <fstream>
-#include <iostream>
 
 using namespace corsika::units::si;
 using namespace corsika::history;

Stack/History/HistoryObservationPlane.hpp

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#pragma once
+
+#include <corsika/geometry/Plane.h>
+#include <corsika/process/ContinuousProcess.h>
+#include <corsika/setup/SetupStack.h>
+#include <corsika/setup/SetupTrajectory.h>
+#include <corsika/units/PhysicalUnits.h>
+
+#include <boost/histogram.hpp>
+
+#include <functional>
+
+namespace corsika::history {
+  namespace detail {
+    auto hist_factory() {
+      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,
+                                          boost::histogram::axis::option::growth_t>{
+              0, 10, "hadronic generation"},
+          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>{});
+      return h;
+    }
+  } // namespace detail
+
+  class HistoryObservationPlane
+      : public corsika::process::ContinuousProcess<HistoryObservationPlane> {
+  public:
+    HistoryObservationPlane(geometry::Plane const&, bool = true);
+
+    void save(std::string const&);
+
+    corsika::units::si::LengthType MaxStepLength(
+        corsika::setup::Stack::ParticleType const&,
+        corsika::setup::Trajectory const& vTrajectory);
+
+    corsika::process::EProcessReturn DoContinuous(
+        corsika::setup::Stack::ParticleType const& vParticle,
+        corsika::setup::Trajectory const& vTrajectory);
+
+  private:
+    void fillHistoryHistogram(setup::Stack::ParticleType const&);
+
+    geometry::Plane const plane_;
+    bool const deleteOnHit_;
+
+    decltype(detail::hist_factory()) histogram_ = detail::hist_factory();
+  };
+} // namespace corsika::history

Stack/History/HistorySecondaryView.hpp

@@ -27,16 +27,11 @@ namespace corsika::history {
 
   public:
     HistorySecondaryView(StackIteratorValue& p)
-        : TView(p)
-        , event_{p.GetEvent()} {
-      if (event_ == nullptr) {
-        // aha, this particle has no [registered] anchestor
-        // thus, create Event here:
-        p.SetEvent(std::make_shared<Event>());
-        std::cout << "Event created for index=" << p.GetIndex() << std::endl;
-        event_ = p.GetEvent();
-      }
+        : TView{p}
+        , event_{std::make_shared<Event>()} {
       event_->setProjectileIndex(p.GetIndex());
+      event_->
+
       //	event_{std::make_shared<Event>(p.GetIndex())} {
       // p.SetEvent(event_); // here an entry on the main particle stack obtains its Event
       // RU: what seems to missing to me right now, at 2am..., is the
@@ -48,9 +43,9 @@ namespace corsika::history {
       auto sec = TView::AddSecondary(std::forward<Args...>(args...));
       // generate new Event for all secondaries to link them to
       // anchestor (aka projectile, here).
-      auto sec_event = std::make_shared<Event>();
+      /*auto sec_event = std::make_shared<Event>();
       sec_event->setParentEventAndSecondaryIndex(event_, event_->secondaries().size());
-      sec.SetEvent(sec_event);
+      sec.SetEvent(sec_event);*/
 
       // store particles at production time in parent/projectile Event here
       event_->addSecondary(sec.GetEnergy(), sec.GetMomentum(), sec.GetPID());
@@ -61,12 +56,6 @@ namespace corsika::history {
 
       return sec;
     }
-
-    // it is probably better to have one method "GetEvent()" and then one can always call
-    // GetEvent().set/getWhatever(...)
-    void SetTarget(const particles::Code targetCode) {
-      event_->setTargetCode(targetCode);
-    }
   };
 
 } // namespace corsika::history

Stack/History/SecondaryParticle.hpp

@@ -18,7 +18,7 @@ namespace corsika::history {
 
   /**
    * This class stores the non-common properties of secondaries in an event. All
-   * other (common) properties are available via the event itself.
+   * other (common) properties are available via the event itself or its projectile.
    */
   class SecondaryParticle {
     units::si::HEPEnergyType const energy_;