From a7738b9ce57feba6c227fbaa2e82a91cb220d3cf Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Sat, 19 Sep 2020 01:50:04 +0200
Subject: [PATCH] make first almost real version of history work
---
Setup/CMakeLists.txt | 2 +-
Setup/SetupStack.h | 3 +-
Stack/History/CMakeLists.txt | 10 +-
Stack/History/Event.hpp | 21 +++-
Stack/History/HSecondaryView.hpp | 41 +++----
Stack/History/HistoryStackExtension.h | 2 +-
.../{testHistory.cc => testHistoryStack.cc} | 31 ++---
Stack/History/testHistoryView.cc | 106 ++++++++++++++++++
8 files changed, 174 insertions(+), 42 deletions(-)
rename Stack/History/{testHistory.cc => testHistoryStack.cc} (69%)
create mode 100644 Stack/History/testHistoryView.cc
diff --git a/Setup/CMakeLists.txt b/Setup/CMakeLists.txt
index e198d90cd..f8490cc2f 100644
--- a/Setup/CMakeLists.txt
+++ b/Setup/CMakeLists.txt
@@ -20,7 +20,7 @@ target_link_libraries (
CORSIKAgeometry
NuclearStackExtension
GeometryNodeStackExtension
- # CORSIKAhistory
+ CORSIKAhistory
)
target_include_directories (
diff --git a/Setup/SetupStack.h b/Setup/SetupStack.h
index aeef45dce..5c22321b0 100644
--- a/Setup/SetupStack.h
+++ b/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<
diff --git a/Stack/History/CMakeLists.txt b/Stack/History/CMakeLists.txt
index 63b6b602c..d2ca77365 100644
--- a/Stack/History/CMakeLists.txt
+++ b/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
)
diff --git a/Stack/History/Event.hpp b/Stack/History/Event.hpp
index ebb5d3a2a..be8d72d4e 100644
--- a/Stack/History/Event.hpp
+++ b/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
diff --git a/Stack/History/HSecondaryView.hpp b/Stack/History/HSecondaryView.hpp
index 1abf8cf6b..4552cc2df 100644
--- a/Stack/History/HSecondaryView.hpp
+++ b/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
diff --git a/Stack/History/HistoryStackExtension.h b/Stack/History/HistoryStackExtension.h
index 58de264c1..294a53780 100644
--- a/Stack/History/HistoryStackExtension.h
+++ b/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);
diff --git a/Stack/History/testHistory.cc b/Stack/History/testHistoryStack.cc
similarity index 69%
rename from Stack/History/testHistory.cc
rename to Stack/History/testHistoryStack.cc
index ebc201ebe..fcb9934e4 100644
--- a/Stack/History/testHistory.cc
+++ b/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);
}
diff --git a/Stack/History/testHistoryView.cc b/Stack/History/testHistoryView.cc
new file mode 100644
index 000000000..673921696
--- /dev/null
+++ b/Stack/History/testHistoryView.cc
@@ -0,0 +1,106 @@
+/*
+ * (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);
+}
--
GitLab