From a0fb1ca3ad29e50685041cb187093dc0bc783af6 Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Tue, 20 Oct 2020 15:01:26 +0200
Subject: [PATCH] added stack-inspector
---
.../modules/sibyll/NuclearInteraction.inl | 2 +-
.../stack_inspector/StackInspector.inl | 92 +++++++++++++++++++
.../stack_inspector/StackInspector.hpp | 53 +++++++++++
.../{testCascade.h => testCascade.hpp} | 9 +-
tests/modules/CMakeLists.txt | 4 +-
tests/modules/testStackInspector.cpp | 35 ++++---
6 files changed, 169 insertions(+), 26 deletions(-)
create mode 100644 corsika/detail/modules/stack_inspector/StackInspector.inl
create mode 100644 corsika/modules/stack_inspector/StackInspector.hpp
rename tests/framework/{testCascade.h => testCascade.hpp} (78%)
diff --git a/corsika/detail/modules/sibyll/NuclearInteraction.inl b/corsika/detail/modules/sibyll/NuclearInteraction.inl
index 5f08a662f..5f1c86225 100644
--- a/corsika/detail/modules/sibyll/NuclearInteraction.inl
+++ b/corsika/detail/modules/sibyll/NuclearInteraction.inl
@@ -493,7 +493,7 @@ namespace corsika::sibyll {
fragm_(kATarget, kAProj, nIntProj, impactPar, nFragments, AFragments);
// this should not occur but well :)
- if (nFragments > GetMaxNFragments())
+ if (nFragments > (int)GetMaxNFragments())
throw std::runtime_error("Number of nuclear fragments in NUCLIB exceeded!");
std::cout << "number of fragments: " << nFragments << std::endl;
diff --git a/corsika/detail/modules/stack_inspector/StackInspector.inl b/corsika/detail/modules/stack_inspector/StackInspector.inl
new file mode 100644
index 000000000..3a49d683d
--- /dev/null
+++ b/corsika/detail/modules/stack_inspector/StackInspector.inl
@@ -0,0 +1,92 @@
+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * 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/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/RootCoordinateSystem.hpp>
+#include <corsika/modules/stack_inspector/StackInspector.hpp>
+
+#include <corsika/setup/SetupTrajectory.hpp>
+
+#include <chrono>
+#include <iomanip>
+#include <iostream>
+#include <limits>
+
+namespace corsika::stack_inspector {
+
+ template <typename TStack>
+ StackInspector<TStack>::StackInspector(const int vNStep, const bool vReportStack,
+ const corsika::units::si::HEPEnergyType vE0)
+ : StackProcess<StackInspector<TStack>>(vNStep)
+ , ReportStack_(vReportStack)
+ , E0_(vE0)
+ , StartTime_(std::chrono::system_clock::now()) {}
+
+ template <typename TStack>
+ StackInspector<TStack>::~StackInspector() {}
+
+ template <typename TStack>
+ corsika::EProcessReturn StackInspector<TStack>::DoStack(const TStack& vS) {
+
+ using namespace units::si;
+
+ [[maybe_unused]] int i = 0;
+ HEPEnergyType Etot = 0_GeV;
+
+ for (const auto& iterP : vS) {
+ units::si::HEPEnergyType E = iterP.GetEnergy();
+ Etot += E;
+ if (ReportStack_) {
+ corsika::CoordinateSystem& rootCS =
+ corsika::RootCoordinateSystem::GetInstance()
+ .GetRootCoordinateSystem(); // for printout
+ auto pos = iterP.GetPosition().GetCoordinates(rootCS);
+ std::cout << "StackInspector: i=" << std::setw(5) << std::fixed << (i++)
+ << ", id=" << std::setw(30) << iterP.GetPID() << " E=" << std::setw(15)
+ << std::scientific << (E / 1_GeV) << " GeV, "
+ << " pos=" << pos << " node = " << iterP.GetNode();
+ if (iterP.GetPID() == Code::Nucleus)
+ std::cout << " nuc_ref=" << iterP.GetNucleusRef();
+ std::cout << std::endl;
+ }
+ }
+
+ auto const now = std::chrono::system_clock::now();
+ const std::chrono::duration<double> elapsed_seconds = now - StartTime_;
+ std::time_t const now_time = std::chrono::system_clock::to_time_t(now);
+ auto const dE = E0_ - Etot;
+ if (dE < dE_threshold_) return corsika::EProcessReturn::eOk;
+ double const progress = dE / E0_;
+
+ double const eta_seconds = elapsed_seconds.count() / progress;
+ std::time_t const eta_time = std::chrono::system_clock::to_time_t(
+ StartTime_ + std::chrono::seconds((int)eta_seconds));
+
+ std::cout << "StackInspector: "
+ << " time=" << std::put_time(std::localtime(&now_time), "%T")
+ << ", running=" << elapsed_seconds.count() << " seconds"
+ << " (" << std::setw(3) << int(progress * 100) << "%)"
+ << ", nStep=" << GetStep() << ", stackSize=" << vS.GetSize()
+ << ", Estack=" << Etot / 1_GeV << " GeV"
+ << ", ETA=" << std::put_time(std::localtime(&eta_time), "%T") << std::endl;
+
+ return corsika::EProcessReturn::eOk;
+ }
+
+ template <typename TStack>
+ void StackInspector<TStack>::Init() {
+ ReportStack_ = false;
+ StartTime_ = std::chrono::system_clock::now();
+ }
+
+} // namespace corsika::stack_inspector
\ No newline at end of file
diff --git a/corsika/modules/stack_inspector/StackInspector.hpp b/corsika/modules/stack_inspector/StackInspector.hpp
new file mode 100644
index 000000000..91fd5271f
--- /dev/null
+++ b/corsika/modules/stack_inspector/StackInspector.hpp
@@ -0,0 +1,53 @@
+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * 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/framework/sequence/StackProcess.hpp>
+#include <corsika/setup/SetupTrajectory.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+
+#include <chrono>
+
+namespace corsika::stack_inspector {
+
+ template <typename TStack>
+ class StackInspector : public corsika::StackProcess<StackInspector<TStack>> {
+
+ typedef typename TStack::ParticleType Particle;
+
+ using corsika::StackProcess<StackInspector<TStack>>::GetStep;
+
+ public:
+ StackInspector(const int vNStep, const bool vReportStack,
+ const corsika::units::si::HEPEnergyType vE0);
+ ~StackInspector();
+
+ void Init();
+ EProcessReturn DoStack(const TStack&);
+
+ /**
+ * To set a new E0, for example when a new shower event is started
+ */
+ void SetE0(const corsika::units::si::HEPEnergyType vE0) { E0_ = vE0; }
+
+ private:
+ bool ReportStack_;
+ corsika::units::si::HEPEnergyType E0_;
+ const corsika::units::si::HEPEnergyType dE_threshold_ = std::invoke([]() {
+ using namespace units::si;
+ return 1_eV;
+ });
+ decltype(std::chrono::system_clock::now()) StartTime_;
+ };
+
+} // namespace corsika::stack_inspector
+
+#include <corsika/detail/modules/stack_inspector/StackInspector.inl>
diff --git a/tests/framework/testCascade.h b/tests/framework/testCascade.hpp
similarity index 78%
rename from tests/framework/testCascade.h
rename to tests/framework/testCascade.hpp
index 715012b9b..895a95820 100644
--- a/tests/framework/testCascade.h
+++ b/tests/framework/testCascade.hpp
@@ -9,11 +9,10 @@
#pragma once
#include <corsika/media/Environment.hpp>
-
#include <corsika/setup/SetupStack.hpp>
-using TestEnvironmentInterface = corsika::environment::IEmpty;
-using TestEnvironmentType = corsika::environment::Environment<TestEnvironmentInterface>;
+using TestEnvironmentType =
+ corsika::Environment<corsika::IMediumModel>;
template <typename T>
using SetupGeometryDataInterface =
@@ -22,11 +21,11 @@ using SetupGeometryDataInterface =
// combine particle data stack with geometry information for tracking
template <typename StackIter>
using StackWithGeometryInterface =
- corsika::CombinedParticleInterface<corsika::detail::ParticleDataStack::PIType,
+ corsika::CombinedParticleInterface<corsika::setup::detail::ParticleDataStack::PIType,
SetupGeometryDataInterface, StackIter>;
using TestCascadeStack =
- corsika::CombinedStack<typename corsika::detail::ParticleDataStack::StackImpl,
+ corsika::CombinedStack<typename corsika::setup::detail::ParticleDataStack::StackImpl,
GeometryData<TestEnvironmentType>, StackWithGeometryInterface>;
/*
diff --git a/tests/modules/CMakeLists.txt b/tests/modules/CMakeLists.txt
index 8572b7346..41a0b6ca6 100644
--- a/tests/modules/CMakeLists.txt
+++ b/tests/modules/CMakeLists.txt
@@ -7,8 +7,8 @@ set (test_modules_sources
testParticleCut.cpp
testPythia8.cpp
#testQGSJetII.cpp
- #testStackInspector.cpp
- #testSwitchProcess.cpp
+ testStackInspector.cpp
+ testSwitchProcess.cpp
#testTrackingLine.cpp
#testUrQMD.cpp
)
diff --git a/tests/modules/testStackInspector.cpp b/tests/modules/testStackInspector.cpp
index 1371f8bfa..640992f9f 100644
--- a/tests/modules/testStackInspector.cpp
+++ b/tests/modules/testStackInspector.cpp
@@ -10,39 +10,38 @@
#include <catch2/catch.hpp>
-#include <corsika/process/stack_inspector/StackInspector.h>
+#include <corsika/modules/stack_inspector/StackInspector.hpp>
-#include <corsika/geometry/Point.h>
-#include <corsika/geometry/RootCoordinateSystem.h>
-#include <corsika/geometry/Vector.h>
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/geometry/RootCoordinateSystem.hpp>
+#include <corsika/framework/geometry/Vector.hpp>
-#include <corsika/units/PhysicalUnits.h>
+#include <corsika/framework/core/PhysicalUnits.hpp>
-#include <corsika/cascade/testCascade.h>
+#include <../framework/testCascade.hpp>
using namespace corsika::units::si;
-using namespace corsika::process::stack_inspector;
+using namespace corsika::stack_inspector;
using namespace corsika;
-using namespace corsika::geometry;
TEST_CASE("StackInspector", "[processes]") {
auto const& rootCS =
- geometry::RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
- geometry::Point const origin(rootCS, {0_m, 0_m, 0_m});
- geometry::Vector<units::si::SpeedType::dimension_type> v(rootCS, 0_m / second,
+ RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
+ Point const origin(rootCS, {0_m, 0_m, 0_m});
+ Vector<units::si::SpeedType::dimension_type> v(rootCS, 0_m / second,
0_m / second, 1_m / second);
- geometry::Line line(origin, v);
- geometry::Trajectory<geometry::Line> track(line, 10_s);
+ Line line(origin, v);
+ Trajectory<Line> track(line, 10_s);
TestCascadeStack stack;
stack.Clear();
HEPEnergyType E0 = 100_GeV;
stack.AddParticle(
- std::tuple<particles::Code, units::si::HEPEnergyType,
- corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
- particles::Code::Electron, E0,
- corsika::stack::MomentumVector(rootCS, {0_GeV, 0_GeV, -1_GeV}),
+ std::tuple<corsika::Code, units::si::HEPEnergyType,
+ corsika::MomentumVector, Point, units::si::TimeType>{
+ Code::Electron, E0,
+ corsika::MomentumVector(rootCS, {0_GeV, 0_GeV, -1_GeV}),
Point(rootCS, {0_m, 0_m, 10_km}), 0_ns});
SECTION("interface") {
@@ -50,6 +49,6 @@ TEST_CASE("StackInspector", "[processes]") {
StackInspector<TestCascadeStack> model(1, true, E0);
model.Init();
- [[maybe_unused]] const process::EProcessReturn ret = model.DoStack(stack);
+ [[maybe_unused]] const corsika::EProcessReturn ret = model.DoStack(stack);
}
}
--
GitLab