From 2d1bbeb6d0e676c797eefb62b51ea12c36b32438 Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Sun, 20 Sep 2020 12:48:18 +0200
Subject: [PATCH] history works now, in principle and mostly, see
testHistoryView
---
Framework/StackInterface/SecondaryView.h | 4 +-
Stack/History/CMakeLists.txt | 12 +--
Stack/History/Event.hpp | 33 ++++--
Stack/History/HSecondaryView.hpp | 51 ---------
Stack/History/HistorySecondaryView.hpp | 72 +++++++++++++
Stack/History/testHistoryView.cc | 129 ++++++++++++++++++++---
6 files changed, 221 insertions(+), 80 deletions(-)
delete mode 100644 Stack/History/HSecondaryView.hpp
create mode 100644 Stack/History/HistorySecondaryView.hpp
diff --git a/Framework/StackInterface/SecondaryView.h b/Framework/StackInterface/SecondaryView.h
index dc7d33ae7..4a33b404c 100644
--- a/Framework/StackInterface/SecondaryView.h
+++ b/Framework/StackInterface/SecondaryView.h
@@ -126,13 +126,13 @@ namespace corsika::stack {
}
template <typename... Args>
- auto AddSecondary(const Args... v) {
+ StackIterator AddSecondary(const Args... v) {
StackIterator proj = GetProjectile();
return AddSecondary(proj, v...);
}
template <typename... Args>
- auto AddSecondary(StackIterator& proj, const Args... v) {
+ StackIterator AddSecondary(StackIterator& proj, const Args... v) {
// make space on stack
InnerStackType::GetStackData().IncrementSize();
innerstack_.deleted_.push_back(false);
diff --git a/Stack/History/CMakeLists.txt b/Stack/History/CMakeLists.txt
index d2ca77365..9d0365b2d 100644
--- a/Stack/History/CMakeLists.txt
+++ b/Stack/History/CMakeLists.txt
@@ -1,18 +1,18 @@
set (
- SETUP_HEADERS
+ HISTORY_HEADERS
Event.hpp
- HSecondaryView.hpp
+ HistorySecondaryView.hpp
HistoryStackExtension.h
SecondaryParticle.hpp
)
set (
- SETUP_NAMESPACE
+ HISTORY_NAMESPACE
corsika/history
)
add_library (CORSIKAhistory INTERFACE)
-CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKAhistory ${SETUP_NAMESPACE} ${SETUP_HEADERS})
+CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKAhistory ${HISTORY_NAMESPACE} ${HISTORY_HEADERS})
target_link_libraries (
CORSIKAhistory
@@ -31,8 +31,8 @@ target_include_directories (
)
install (
- FILES ${SETUP_HEADERS}
- DESTINATION include/${SETUP_NAMESPACE}
+ FILES ${HISTORY_HEADERS}
+ DESTINATION include/${HISTORY_NAMESPACE}
)
# ----------------
diff --git a/Stack/History/Event.hpp b/Stack/History/Event.hpp
index be8d72d4e..6adc08ddc 100644
--- a/Stack/History/Event.hpp
+++ b/Stack/History/Event.hpp
@@ -18,31 +18,46 @@
namespace corsika::history {
- struct Event {
- size_t const projectileIndex_; //!< reference to projectile
+ class Event;
+ using EvtPtr = std::shared_ptr<history::Event>;
+
+ class Event {
+
+ size_t projectileIndex_ = 0; //!< reference to projectile (for secondaries_)
std::vector<SecondaryParticle> secondaries_;
- //std::shared_ptr<Event> parent_event_;
+ EvtPtr parent_event_;
+ size_t sec_index_ = 0;
- // meta information, could also be in a separate class
std::optional<corsika::particles::Code>
targetCode_; // cannot be const, value set only after construction
public:
- Event(const size_t projectileIndex)
- : projectileIndex_{projectileIndex} {
- std::cout << "Event created (index = " << projectileIndex_ << ")" << std::endl;
+ Event() = default;
+
+ 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_;
}
+
+ EvtPtr parentEvent() { return parent_event_; }
+
+ void setProjectileIndex(size_t i) { projectileIndex_=i; }
+ size_t projectileIndex() const { return 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; }
+
};
-
- using EvtPtr = std::shared_ptr<history::Event>;
} // namespace corsika::history
diff --git a/Stack/History/HSecondaryView.hpp b/Stack/History/HSecondaryView.hpp
deleted file mode 100644
index 4552cc2df..000000000
--- a/Stack/History/HSecondaryView.hpp
+++ /dev/null
@@ -1,51 +0,0 @@
-/*
- * (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.
- */
-
-#pragma once
-
-#include <corsika/stack/SecondaryView.h>
-#include <corsika/history/Event.hpp>
-
-#include <memory>
-#include <type_traits>
-#include <utility>
-
-namespace corsika::history {
-
- template <typename TView>
- class HSecondaryView : public TView {
-
- std::shared_ptr<Event> event_;
-
- using StackIterator = typename TView::StackIteratorValue;
-
- public:
- 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 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());
- }
-
- // 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/HistorySecondaryView.hpp b/Stack/History/HistorySecondaryView.hpp
new file mode 100644
index 000000000..16603e948
--- /dev/null
+++ b/Stack/History/HistorySecondaryView.hpp
@@ -0,0 +1,72 @@
+/*
+ * (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.
+ */
+
+#pragma once
+
+#include <corsika/stack/SecondaryView.h>
+#include <corsika/history/Event.hpp>
+
+#include <memory>
+#include <type_traits>
+#include <utility>
+
+namespace corsika::history {
+
+ template <typename TView>
+ class HistorySecondaryView : public TView {
+
+ EvtPtr event_;
+
+ using StackIteratorValue = typename TView::StackIteratorValue;
+ using StackIterator = typename TView::StackIterator;
+
+ 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();
+ }
+ event_->setProjectileIndex(p.GetIndex());
+ // 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>
+ 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);
+
+ // store particles at production time in parent/projectile Event here
+ event_->addSecondary(sec.GetEnergy(), sec.GetMomentum(), sec.GetPID());
+
+ // 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;
+ }
+
+ // 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/testHistoryView.cc b/Stack/History/testHistoryView.cc
index 673921696..1eb9efe14 100644
--- a/Stack/History/testHistoryView.cc
+++ b/Stack/History/testHistoryView.cc
@@ -8,7 +8,7 @@
#include <corsika/history/HistoryStackExtension.h>
#include <corsika/history/Event.hpp>
-#include <corsika/history/HSecondaryView.hpp>
+#include <corsika/history/HistorySecondaryView.hpp>
#include <corsika/stack/CombinedStack.h>
#include <corsika/stack/dummy/DummyStack.h>
@@ -66,41 +66,146 @@ TEST_CASE("HistoryStackExtension", "[stack]") {
geometry::CoordinateSystem& dummyCS =
geometry::RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
+ // in this test we only use one singel stack !
TestStack s;
- SECTION("add lone particle") {
+ // add primary particle
+ auto p0 = 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});
- auto p = s.AddParticle(
+ CHECK(s.getEntries() == 1);
+ corsika::history::EvtPtr evt = p0.GetEvent();
+ CHECK(evt == nullptr);
+
+ // add secondaries, 1st generation
+ history::HistorySecondaryView<TestStackView> hview0(p0);
+
+ auto ev0 = p0.GetEvent();
+ CHECK(ev0 != nullptr);
+ // CHECK(ev0->projectile(s.begin()) == p0);
+
+ // add 5 secondaries
+ for (int i = 0; i < 5; ++i) {
+ auto sec = hview0.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});
+
+ CHECK(sec.GetEvent()->parentEvent() == ev0);
+ CHECK(sec.GetEvent()->parentEvent()->parentEvent() == nullptr);
+ }
+
+ // read 1st genertion particle particle
+ auto p1 = s.GetNextParticle();
+
+ history::HistorySecondaryView<TestStackView> hview1(p1);
+
+ auto ev1 = p1.GetEvent();
+
+ // add second generation of secondaries
+ // add 10 secondaries
+ for (int i = 0; i < 10; ++i) {
+ auto sec = hview1.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});
- CHECK(s.GetSize() == 1);
- corsika::history::EvtPtr evt = p.GetEvent();
- CHECK(evt == nullptr);
+ CHECK(sec.GetEvent()->parentEvent() == ev1);
+ CHECK(sec.GetEvent()->parentEvent()->parentEvent() == ev0);
+ CHECK(sec.GetEvent()->parentEvent()->parentEvent()->parentEvent() == nullptr);
}
- SECTION("generate event") {
- auto p = s.AddParticle(
+ // read 2nd genertion particle particle
+ auto p2 = s.GetNextParticle();
+
+ history::HistorySecondaryView<TestStackView> hview2(p2);
+
+ auto ev2 = p2.GetEvent();
+
+ // add third generation of secondaries
+ // add 15 secondaries
+ for (int i = 0; i < 15; ++i) {
+ auto sec = hview2.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});
- history::HSecondaryView<TestStackView> hview(p);
+ CHECK(sec.GetEvent()->parentEvent() == ev2);
+ CHECK(sec.GetEvent()->parentEvent()->parentEvent() == ev1);
+ CHECK(sec.GetEvent()->parentEvent()->parentEvent()->parentEvent() == ev0);
+ CHECK(sec.GetEvent()->parentEvent()->parentEvent()->parentEvent()->parentEvent() ==
+ nullptr);
+ }
+
+ // read 3rd genertion particle particle
+ auto p3 = s.GetNextParticle();
+
+ history::HistorySecondaryView<TestStackView> hview3(p3);
+
+ auto ev3 = p3.GetEvent();
- hview.AddSecondary(
+ // add fourth generation of secondaries
+ // add 20 secondaries
+ for (int i = 0; i < 20; ++i) {
+ auto sec = hview3.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 ...
+ CHECK(sec.GetEvent()->parentEvent() == ev3);
+ CHECK(sec.GetEvent()->parentEvent()->parentEvent() == ev2);
+ CHECK(sec.GetEvent()->parentEvent()->parentEvent()->parentEvent() == ev1);
+ CHECK(sec.GetEvent()->parentEvent()->parentEvent()->parentEvent()->parentEvent() ==
+ ev0);
+ CHECK(sec.GetEvent()
+ ->parentEvent()
+ ->parentEvent()
+ ->parentEvent()
+ ->parentEvent()
+ ->parentEvent() == nullptr);
}
- // REQUIRE(pout.GetPID() == particles::Code::Electron);
+ /*
+ Now, let's perform some history reading and checking based on p3
+
+ p3 should have 20 secondaries, and a projectile (with 15
+ secondaries), with another projectil (with 10 secondaries), with
+ antother projectil (5 secondaries), with NO parent
+
+ **/
+ {
+ auto test_ev3 = p3.GetEvent();
+ auto test_sec3 = test_ev3->secondaries();
+ CHECK(test_sec3.size() == 20);
+
+ auto test_proj3 = test_ev3->projectile(s.begin());
+ CHECK(test_proj3.GetEvent() == ev3);
+ auto test_ev2 = test_ev3->parentEvent();
+ auto test_sec2 = test_ev2->secondaries();
+ CHECK(test_sec2.size() == 15);
+
+ auto test_proj2 = test_ev2->projectile(s.begin());
+ CHECK(test_proj2.GetEvent() == ev2);
+ auto test_ev1 = test_ev2->parentEvent();
+ auto test_sec1 = test_ev1->secondaries();
+ CHECK(test_sec1.size() == 10);
+
+ auto test_proj1 = test_ev1->projectile(s.begin());
+ CHECK(test_proj1.GetEvent() == ev1);
+
+ CHECK(test_proj1.GetEvent()->parentEvent() == ev0);
+ CHECK(test_proj1.GetEvent()->parentEvent()->parentEvent() == nullptr);
+ }
}
--
GitLab