diff --git a/Documentation/Examples/CMakeLists.txt b/Documentation/Examples/CMakeLists.txt index 3e933ea6f97282fb6a5f0f08a55f6967b8d1a996..75b7176085f2e6d1cafdc27a88445b511b86fab8 100644 --- a/Documentation/Examples/CMakeLists.txt +++ b/Documentation/Examples/CMakeLists.txt @@ -91,6 +91,11 @@ if (Pythia8_FOUND) CORSIKAunits CORSIKAlogging CORSIKArandom + CORSIKAhistory + ProcessSibyll + ProcessPythia8 + ProcessUrQMD + ProcessSwitch CORSIKAcascade ProcessProposal ProcessPythia8 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 index d9374eee7e00ee8177cd3027b86a9fac7bc6b9cd..7971124bc99a8bf1b64ae6202a1ebac6fa81575b 100644 --- a/Stack/History/HistoryObservationPlane.cc +++ b/Stack/History/HistoryObservationPlane.cc @@ -12,7 +12,6 @@ #include <boost/histogram/ostream.hpp> #include <fstream> -#include <iostream> using namespace corsika::units::si; using namespace corsika::history; 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 diff --git a/Stack/History/SecondaryParticle.hpp b/Stack/History/SecondaryParticle.hpp index f1a2831d82fe8c7a1071116096b55e766488b71f..b6db7e8fd04684ff9a39285fe15e2b228ea74835 100644 --- a/Stack/History/SecondaryParticle.hpp +++ b/Stack/History/SecondaryParticle.hpp @@ -18,7 +18,7 @@ namespace corsika::history { /** * This class stores the non-common properties of secondaries in an event. All - * other (common) properties are available via the event itself. + * other (common) properties are available via the event itself or its projectile. */ class SecondaryParticle { units::si::HEPEnergyType const energy_;