make first almost real version of history work

Setup/CMakeLists.txt

@@ -20,7 +20,7 @@ target_link_libraries (
   CORSIKAgeometry
   NuclearStackExtension
   GeometryNodeStackExtension
-  #  CORSIKAhistory
+  CORSIKAhistory
   )
 
 target_include_directories (

Setup/SetupStack.h

@@ -8,10 +8,10 @@
 
 #pragma once
 
+#include <corsika/history/HistoryStackExtension.h>
 #include <corsika/stack/CombinedStack.h>
 #include <corsika/stack/node/GeometryNodeStackExtension.h>
 #include <corsika/stack/nuclear_extension/NuclearStackExtension.h>
-#include <corsika/history/HistoryStackExtension.h>
 
 #include <corsika/setup/SetupEnvironment.h>
 
@@ -33,7 +33,6 @@ namespace corsika::setup {
                                     ...event... >::type;
     */
 
-    
     // combine particle data stack with geometry information for tracking
     template <typename TStackIter>
     using StackWithGeometryInterface = corsika::stack::CombinedParticleInterface<

Stack/History/CMakeLists.txt

@@ -37,9 +37,15 @@ install (
 
 # ----------------
 # code unit testing
-CORSIKA_ADD_TEST(testHistory)
+CORSIKA_ADD_TEST(testHistoryStack)
 target_link_libraries (
-  testHistory
+  testHistoryStack
+  CORSIKAhistory
+  CORSIKAtesting
+  )
+CORSIKA_ADD_TEST(testHistoryView)
+target_link_libraries (
+  testHistoryView
   CORSIKAhistory
   CORSIKAtesting
   )

Stack/History/Event.hpp

@@ -14,22 +14,35 @@
 #include <iostream>
 #include <optional>
 #include <vector>
+#include <memory>
 
 namespace corsika::history {
 
   struct Event {
     size_t const projectileIndex_; //!< reference to projectile
-    std::vector<SecondaryParticle> secondaries;
-
+    std::vector<SecondaryParticle> secondaries_;
+    //std::shared_ptr<Event> parent_event_;
+    
     // meta information, could also be in a separate class
     std::optional<corsika::particles::Code>
-        targetCode; // cannot be const, value set only after construction
+        targetCode_; // cannot be const, value set only after construction
 
   public:
-    Event(size_t projectileIndex)
+    Event(const size_t projectileIndex)
         : projectileIndex_{projectileIndex} {
       std::cout << "Event created (index = " << projectileIndex_ << ")" << std::endl;
     }
+    
+    void addSecondary(units::si::HEPEnergyType energy,
+                      geometry::Vector<units::si::hepmomentum_d> momentum,
+                      particles::Code pid) {
+      secondaries_.emplace_back(energy, momentum, pid);
+    }
+    
+    void setTargetCode(const particles::Code t) { targetCode_=t; }
+    
   };
 
+  using EvtPtr = std::shared_ptr<history::Event>;
+  
 } // namespace corsika::history

Stack/History/HSecondaryView.hpp

@@ -17,32 +17,35 @@
 
 namespace corsika::history {
 
-  namespace detail {
-    template <typename TParticleIterator>
-    using SecondaryViewTypeFromIterator =
-        decltype(stack::SecondaryView{std::declval<TParticleIterator&>()});
-  }
-
-  template <typename TParticleIterator>
-  class HSecondaryView : public detail::SecondaryViewTypeFromIterator<TParticleIterator> {
+  template <typename TView>
+  class HSecondaryView : public TView {
+    
     std::shared_ptr<Event> event_;
-
-    using BaseSecondaryViewType =
-        detail::SecondaryViewTypeFromIterator<TParticleIterator>;
-
+    
+    using StackIterator = typename TView::StackIteratorValue;
+    
   public:
-    HSecondaryView(TParticleIterator& p)
-        : BaseSecondaryViewType{p}
-        , event_{std::make_shared<Event>(p.GetIndex())} {}
+    HSecondaryView(StackIterator& p)
+      : TView(p),
+	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.
+    }
 
     template <typename... Args>
     void AddSecondary(Args&&... args) {
-      auto const s = BaseSecondaryViewType::AddSecondary(
-          std::forward<Args...>(args...), event_); // what if event is not last argument?
-      event_->secondaries.emplace_back(s.GetEnergy(), s.GetMomentum(), s.GetParticleID());
+      auto s = TView::AddSecondary(
+          std::forward<Args...>(args...));
+      // store particles at production time here
+      // 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. 
+      event_->addSecondary(s.GetEnergy(), s.GetMomentum(), s.GetPID());
     }
 
-    void SetTarget(particles::Code targetCode) { event_->targetCode = targetCode; }
+    // 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/HistoryStackExtension.h

@@ -72,7 +72,7 @@ namespace corsika::history {
     void SetParticleData() {} // nullptr, already by design
 
     // create a new particle as secondary of a parent
-    void SetParticleData(HistoryDataInterface& parent) { SetParticleData(); }
+    void SetParticleData(HistoryDataInterface& /*parent*/) { SetParticleData(); }
 
     void SetEvent(const std::shared_ptr<TEvent>& v) {
       GetStackData().SetEvent(GetIndex(), v);

Stack/History/testHistory.cc → Stack/History/testHistoryStack.cc

@@ -9,8 +9,6 @@
 #include <corsika/history/HistoryStackExtension.h>
 #include <corsika/stack/CombinedStack.h>
 #include <corsika/stack/dummy/DummyStack.h>
-#include <corsika/history/Event.hpp>
-#include <corsika/history/HSecondaryView.hpp>
 
 #include <catch2/catch.hpp>
 
@@ -19,10 +17,25 @@
 using namespace corsika;
 using namespace corsika::stack;
 
+// this is our dummy environment, it only knows its trivial BaseNodeType
+class DummyEvent {
+private:
+  size_t parent_;
+  std::vector<int> secondaries_;
+
+public:
+  DummyEvent() {}
+  DummyEvent(const size_t parent) { parent_ = parent; }
+
+  size_t getParentIndex() { return parent_; }
+  void addSecondary(const int particle) { secondaries_.push_back(particle); }
+  int multiplicity() const { return secondaries_.size(); }
+};
+
 // the GeometryNode stack needs to know the type of geometry-nodes from the DummyEnv:
 template <typename TStackIter>
 using DummyHistoryDataInterface =
-    typename history::MakeHistoryDataInterface<TStackIter, history::Event>::type;
+    typename history::MakeHistoryDataInterface<TStackIter, DummyEvent>::type;
 
 // combine dummy stack with geometry information for tracking
 template <typename TStackIter>
@@ -32,10 +45,10 @@ using StackWithHistoryInterface =
 
 using TestStack =
     corsika::stack::CombinedStack<typename stack::dummy::DummyStack::StackImpl,
-                                  history::HistoryData<history::Event>,
+                                  history::HistoryData<DummyEvent>,
                                   StackWithHistoryInterface>;
 
-using EvtPtr = std::shared_ptr<history::Event>;
+using EvtPtr = std::shared_ptr<DummyEvent>;
 
 TEST_CASE("HistoryStackExtension", "[stack]") {
 
@@ -50,12 +63,4 @@ TEST_CASE("HistoryStackExtension", "[stack]") {
     EvtPtr evt = p.GetEvent();
     CHECK(evt == nullptr);
   }
-
-  SECTION("write event") {
-    history::HSecondaryView hview{p};
-
-    hview.AddSecondary(std::tuple<dummy::NoData>{noData});
-  }
-
-  // REQUIRE(pout.GetPID() == particles::Code::Electron);
 }

Stack/History/testHistoryView.cc

+/*
+ * (c) Copyright 2018 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.
+ */
+
+#include <corsika/history/HistoryStackExtension.h>
+#include <corsika/history/Event.hpp>
+#include <corsika/history/HSecondaryView.hpp>
+
+#include <corsika/stack/CombinedStack.h>
+#include <corsika/stack/dummy/DummyStack.h>
+#include <corsika/stack/nuclear_extension/NuclearStackExtension.h>
+
+#include <catch2/catch.hpp>
+
+#include <iostream>
+
+using namespace corsika;
+using namespace corsika::stack;
+using namespace corsika::geometry;
+using namespace corsika::units::si;
+
+/**
+   Need to replicate setup::SetupStack in a maximally simplified
+   way, but with real particle data
+ */
+
+// the GeometryNode stack needs to know the type of geometry-nodes from the DummyEnv:
+template <typename TStackIter>
+using HistoryDataInterface =
+    typename history::MakeHistoryDataInterface<TStackIter, history::Event>::type;
+
+// combine dummy stack with geometry information for tracking
+template <typename TStackIter>
+using StackWithHistoryInterface = corsika::stack::CombinedParticleInterface<
+    stack::nuclear_extension::ParticleDataStack::PIType, HistoryDataInterface,
+    TStackIter>;
+
+using TestStack = corsika::stack::CombinedStack<
+    typename stack::nuclear_extension::ParticleDataStack::StackImpl,
+    history::HistoryData<history::Event>, StackWithHistoryInterface>;
+
+/*
+    See Issue 161
+
+    unfortunately clang does not support this in the same way (yet) as
+    gcc, so we have to distinguish here. If clang cataches up, we
+    could remove the clang branch here and also in
+    corsika::Cascade. The gcc code is much more generic and
+    universal. If we could do the gcc version, we won't had to define
+    StackView globally, we could do it with MakeView whereever it is
+    actually needed. Keep an eye on this!
+  */
+#if defined(__clang__)
+using TestStackView = corsika::stack::SecondaryView<typename TestStack::StackImpl,
+                                                    StackWithHistoryInterface>;
+#elif defined(__GNUC__) || defined(__GNUG__)
+using TestStackView = corsika::stack::MakeView<TestStack>::type;
+#endif
+
+TEST_CASE("HistoryStackExtension", "[stack]") {
+
+  geometry::CoordinateSystem& dummyCS =
+      geometry::RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
+
+  TestStack s;
+
+  SECTION("add lone particle") {
+
+    auto p = s.AddParticle(
+        std::tuple<particles::Code, units::si::HEPEnergyType,
+                   corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
+            particles::Code::Electron, 1.5_GeV,
+            corsika::stack::MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
+            Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s});
+
+    CHECK(s.GetSize() == 1);
+    corsika::history::EvtPtr evt = p.GetEvent();
+    CHECK(evt == nullptr);
+  }
+
+  SECTION("generate event") {
+    auto p = s.AddParticle(
+        std::tuple<particles::Code, units::si::HEPEnergyType,
+                   corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
+            particles::Code::Electron, 1.5_GeV,
+            corsika::stack::MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
+            Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s});
+
+    history::HSecondaryView<TestStackView> hview(p);
+
+    hview.AddSecondary(
+        std::tuple<particles::Code, units::si::HEPEnergyType,
+                   corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
+            particles::Code::Electron, 1.5_GeV,
+            corsika::stack::MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
+            Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s});
+
+    // ... continue actual real test here an below ...
+  }
+
+  // REQUIRE(pout.GetPID() == particles::Code::Electron);
+}