added stack-inspector

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;

corsika/detail/modules/stack_inspector/StackInspector.inl

+/*
+ * (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

corsika/modules/stack_inspector/StackInspector.hpp

+/*
+ * (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>

tests/framework/testCascade.h → 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>;
 
 /*

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
   )

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);
   }
 }