diff --git a/Documentation/Examples/CMakeLists.txt b/Documentation/Examples/CMakeLists.txt
index 1592372c2c6a346def1c866508841097c4bbd1c2..cfacc09ea7df8944044862d7ca6da632f93fb01e 100644
--- a/Documentation/Examples/CMakeLists.txt
+++ b/Documentation/Examples/CMakeLists.txt
@@ -84,6 +84,7 @@ if (Pythia8_FOUND)
CORSIKAunits
CORSIKAlogging
CORSIKArandom
+ CORSIKAhistory
ProcessSibyll
ProcessPythia8
ProcessUrQMD
diff --git a/Stack/History/CMakeLists.txt b/Stack/History/CMakeLists.txt
index 9d0365b2dccbd2a280ff9c54cac7521b46a844b8..c18823b928c4321387cc6345141d3fdf665c399a 100644
--- a/Stack/History/CMakeLists.txt
+++ b/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 (
diff --git a/Stack/History/Event.hpp b/Stack/History/Event.hpp
index 6adc08ddc8fa1a77cd0868312615b2dfa707c345..2b38aae62ec2ff53623645c02b61c1b67ec5550f 100644
--- a/Stack/History/Event.hpp
+++ b/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
diff --git a/Stack/History/HistoryObservationPlane.cc b/Stack/History/HistoryObservationPlane.cc
new file mode 100644
index 0000000000000000000000000000000000000000..ebc8ab11d989c74975c012bb32f75cdb6a442a6e
--- /dev/null
+++ b/Stack/History/HistoryObservationPlane.cc
@@ -0,0 +1,61 @@
+/*
+ * (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.
+ */
+
+#include <corsika/history/HistoryObservationPlane.hpp>
+
+#include <fstream>
+
+using namespace corsika::units::si;
+using namespace corsika::history;
+using namespace corsika;
+
+HistoryObservationPlane::HistoryObservationPlane(geometry::Plane const& obsPlane,
+ bool deleteOnHit)
+ : plane_(obsPlane)
+ , deleteOnHit_(deleteOnHit) {}
+
+corsika::process::EProcessReturn HistoryObservationPlane::DoContinuous(
+ setup::Stack::ParticleType const& particle, setup::Trajectory const& trajectory) {
+ TimeType const timeOfIntersection =
+ (plane_.GetCenter() - trajectory.GetR0()).dot(plane_.GetNormal()) /
+ trajectory.GetV0().dot(plane_.GetNormal());
+
+ if (timeOfIntersection < TimeType::zero()) { return process::EProcessReturn::eOk; }
+
+ if (plane_.IsAbove(trajectory.GetR0()) == plane_.IsAbove(trajectory.GetPosition(1))) {
+ return process::EProcessReturn::eOk;
+ }
+
+ auto const pid = particle.GetPID();
+ if (particles::IsMuon(pid)) { fillHistoryHistogram(particle); }
+
+ if (deleteOnHit_) {
+ return process::EProcessReturn::eParticleAbsorbed;
+ } else {
+ return process::EProcessReturn::eOk;
+ }
+}
+
+LengthType HistoryObservationPlane::MaxStepLength(setup::Stack::ParticleType const&,
+ setup::Trajectory const& trajectory) {
+ TimeType const timeOfIntersection =
+ (plane_.GetCenter() - trajectory.GetR0()).dot(plane_.GetNormal()) /
+ trajectory.GetV0().dot(plane_.GetNormal());
+
+ if (timeOfIntersection < TimeType::zero()) {
+ return std::numeric_limits<double>::infinity() * 1_m;
+ }
+
+ auto const pointOfIntersection = trajectory.GetPosition(timeOfIntersection);
+ return (trajectory.GetR0() - pointOfIntersection).norm() * 1.0001;
+}
+
+void fillHistoryHistogram(setup::Stack::ParticleType const& muon) {
+ double const muonEnergy = muon.GetEnergy() / 1_eV;
+
+}
diff --git a/Stack/History/HistoryObservationPlane.hpp b/Stack/History/HistoryObservationPlane.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..8c2d5700daa63cce61d9f88705047f6aa2cf663e
--- /dev/null
+++ b/Stack/History/HistoryObservationPlane.hpp
@@ -0,0 +1,61 @@
+/*
+ * (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
diff --git a/Stack/History/HistorySecondaryView.hpp b/Stack/History/HistorySecondaryView.hpp
index 16603e948b3ae5a0dba9e5c9613dfaba2ee88422..365d4b4ea33943c467feca67fda2dcd7d3c5b9b5 100644
--- a/Stack/History/HistorySecondaryView.hpp
+++ b/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