event structure

Stack/History/CMakeLists.txt

@@ -26,6 +26,8 @@ target_link_libraries (
 #  INTERFACE
   CORSIKAparticles
   CORSIKAgeometry
+  CORSIKAprocesssequence # for HistoryObservationPlane
+  CORSIKAsetup # for HistoryObservationPlane
   SuperStupidStack
   NuclearStackExtension
   C8::ext::boost
@@ -50,6 +52,7 @@ target_link_libraries (
   testHistoryStack
   CORSIKAhistory
   CORSIKAtesting
+  DummyStack
   )
 CORSIKA_ADD_TEST(testHistoryView)
 target_link_libraries (

Stack/History/Event.hpp

@@ -19,15 +19,13 @@
 namespace corsika::history {
 
   class Event;
-  using EvtPtr = std::shared_ptr<history::Event>;
+  using EventPtr = std::shared_ptr<history::Event>;
 
   class Event {
 
-    size_t projectileIndex_ = 0; //!< reference to projectile (for secondaries_)
+    size_t projectileIndex_ = 0; //!< index of projectile on stack
     std::vector<SecondaryParticle> secondaries_;
-    EvtPtr parent_event_;
-    size_t sec_index_ =
-        0; //!< index of the projectile of this Event in the parent_event_.secondaries_;
+    EventPtr parent_event_;
 
     std::optional<corsika::particles::Code>
         targetCode_; // cannot be const, value set only after construction
@@ -35,16 +33,9 @@ namespace corsika::history {
   public:
     Event() = default;
 
-    void setParentEvent(EvtPtr& evt) { parent_event_ = evt; }
+    void setParentEvent(EventPtr const& evt) { parent_event_ = evt; }
 
-    void setParentEventAndSecondaryIndex(EvtPtr& evt, size_t sec_index) {
-      setParentEvent(evt);
-      setParentSecondaryIndex(sec_index);
-    }
-
-    void setParentSecondaryIndex(size_t sec_index) { sec_index_ = sec_index; }
-
-    EvtPtr parentEvent() { return parent_event_; }
+    EventPtr parentEvent() { return parent_event_; }
 
     void setProjectileIndex(size_t i) { projectileIndex_ = i; }
     size_t projectileIndex() const { return projectileIndex_; }
@@ -54,11 +45,13 @@ namespace corsika::history {
       return begin + projectileIndex_;
     }
 
-    void addSecondary(units::si::HEPEnergyType energy,
-                      geometry::Vector<units::si::hepmomentum_d> momentum,
-                      particles::Code pid) {
+    size_t addSecondary(units::si::HEPEnergyType energy,
+                        geometry::Vector<units::si::hepmomentum_d> const& momentum,
+                        particles::Code pid) {
       secondaries_.emplace_back(energy, momentum, pid);
+      return secondaries_.size() - 1;
     }
+
     std::vector<SecondaryParticle>& secondaries() { return secondaries_; }
 
     void setTargetCode(const particles::Code t) { targetCode_ = t; }

Stack/History/HistorySecondaryView.hpp

@@ -20,7 +20,7 @@ namespace corsika::history {
   template <typename TView>
   class HistorySecondaryView : public TView {
 
-    EvtPtr event_;
+    EventPtr event_;
 
     using StackIteratorValue = typename TView::StackIteratorValue;
     using StackIterator = typename TView::StackIterator;
@@ -30,31 +30,24 @@ namespace corsika::history {
         : 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
-      // actual back reference to the parent event. This needs to be added.
+      event_->setParentEvent(p.GetEvent());
     }
 
     template <typename... Args>
     StackIterator AddSecondary(Args&&... args) {
-      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>();
-      sec_event->setParentEventAndSecondaryIndex(event_, event_->secondaries().size());
-      sec.SetEvent(sec_event);*/
+      auto stack_sec = TView::AddSecondary(std::forward<Args...>(args...));
 
-      // store particles at production time in parent/projectile Event here
-      event_->addSecondary(sec.GetEnergy(), sec.GetMomentum(), sec.GetPID());
+      // store particles at production time in Event here
+      auto const sec_index = event_->addSecondary(
+          stack_sec.GetEnergy(), stack_sec.GetMomentum(), stack_sec.GetPID());
+      stack_sec.SetParentEventIndex(sec_index);
+      stack_sec.SetEvent(event_);
 
       // RU: consider if we can call
       // TView::AddSecondary twice instead: 1. for particle at production time
       // , 2. dynamic particle ... not sure, but would be extremely flexible.
 
-      return sec;
+      return stack_sec;
     }
   };
 

Stack/History/HistoryStackExtension.h

@@ -11,7 +11,7 @@
 #include <corsika/stack/Stack.h>
 
 #include <memory>
-#include <tuple>
+#include <utility>
 #include <vector>
 
 namespace corsika::history {
@@ -25,28 +25,41 @@ namespace corsika::history {
    */
   template <typename TEvent>
   class HistoryData {
+    using EventPtr =
+        std::shared_ptr<TEvent>; //!< Pointer to the event where this particle was created
+    using ParentEventIndex = int; //!< index to TEvent::secondaries_
+    using DataType = std::pair<EventPtr, ParentEventIndex>;
 
   public:
     // these functions are needed for the Stack interface
-    void Clear() { event_.clear(); }
-    unsigned int GetSize() const { return event_.size(); }
-    unsigned int GetCapacity() const { return event_.size(); }
-    void Copy(const int i1, const int i2) { event_[i2] = event_[i1]; }
-    void Swap(const int i1, const int i2) { std::swap(event_[i1], event_[i2]); }
+    void Clear() { historyData_.clear(); }
+    unsigned int GetSize() const { return historyData_.size(); }
+    unsigned int GetCapacity() const { return historyData_.size(); }
+    void Copy(const int i1, const int i2) { historyData_[i2] = historyData_[i1]; }
+    void Swap(const int i1, const int i2) {
+      std::swap(historyData_[i1], historyData_[i2]);
+    }
 
     // custom data access function
-    void SetEvent(const int i, std::shared_ptr<TEvent> v) { event_[i] = std::move(v); }
-    std::shared_ptr<TEvent> GetEvent(const int i) const { return event_[i]; }
+    void SetEvent(const int i, EventPtr v) { historyData_[i].first = std::move(v); }
+    EventPtr GetEvent(const int i) const { return historyData_[i].first; }
+
+    void SetParentEventIndex(const int i, ParentEventIndex v) {
+      historyData_[i].second = std::move(v);
+    }
+    ParentEventIndex GetParentEventIndex(const int i) const {
+      return historyData_[i].second;
+    }
 
     // these functions are also needed by the Stack interface
-    void IncrementSize() { event_.push_back(nullptr); }
+    void IncrementSize() { historyData_.push_back(DataType{}); }
     void DecrementSize() {
-      if (event_.size() > 0) { event_.pop_back(); }
+      if (historyData_.size() > 0) { historyData_.pop_back(); }
     }
 
     // custom private data section
   private:
-    std::vector<std::shared_ptr<TEvent>> event_;
+    std::vector<DataType> historyData_;
   };
 
   /**
@@ -68,7 +81,7 @@ namespace corsika::history {
 
   public:
     // create a new particle from scratch
-    void SetParticleData() {} // nullptr, already by design
+    void SetParticleData() { GetStackData().SetParentEventIndex(GetIndex(), -1); }
 
     // create a new particle as secondary of a parent
     void SetParticleData(HistoryDataInterface& /*parent*/) { SetParticleData(); }
@@ -77,9 +90,17 @@ namespace corsika::history {
       GetStackData().SetEvent(GetIndex(), v);
     }
 
+    void SetParentEventIndex(int index) {
+      GetStackData().SetParentEventIndex(GetIndex(), index);
+    }
+
     std::shared_ptr<TEvent> GetEvent() const {
       return GetStackData().GetEvent(GetIndex());
     }
+
+    int GetParentEventIndex() const {
+      return GetStackData().GetParentEventIndex(GetIndex());
+    }
   };
 
   template <typename T, typename TEvent>
@@ -89,9 +110,8 @@ namespace corsika::history {
 
 } // namespace corsika::history
 
-  
-  // for user-friendlyness we create the HistoryDataInterface type
-  // with the histoy::Event data content right here: 
+// for user-friendlyness we create the HistoryDataInterface type
+// with the histoy::Event data content right here:
 
 #include <corsika/history/Event.hpp>
 
@@ -99,8 +119,8 @@ namespace corsika::history {
 
   template <typename TStackIter>
   using HistoryEventDataInterface =
-  typename history::MakeHistoryDataInterface<TStackIter, history::Event>::type;  
+      typename history::MakeHistoryDataInterface<TStackIter, history::Event>::type;
 
   using HistoryEventData = history::HistoryData<history::Event>;
-  
+
 } // namespace corsika::history