Merge branch '87-add-energy-conservation-bookkeeping' into 'master'

Resolve "Add energy conservation bookkeeping"

Closes #87

See merge request !101

Documentation/Examples/CMakeLists.txt

@@ -68,6 +68,7 @@ if (Pythia8_FOUND)
     CORSIKAcascade
     ProcessEnergyLoss
     ProcessTrackWriter
+    ProcessStackInspector
     ProcessTrackingLine
     ProcessParticleCut
     ProcessHadronicElasticModel

Documentation/Examples/cascade_example.cc

@@ -29,7 +29,6 @@
 #include <corsika/process/sibyll/NuclearInteraction.h>
 
 #include <corsika/process/track_writer/TrackWriter.h>
-
 #include <corsika/process/particle_cut/ParticleCut.h>
 
 #include <corsika/units/PhysicalUnits.h>
@@ -93,27 +92,6 @@ int main() {
 
   universe.AddChild(std::move(outerMedium));
 
-  // setup processes, decays and interactions
-  tracking_line::TrackingLine tracking;
-  stack_inspector::StackInspector<setup::Stack> stackInspect(1, true);
-
-  const std::vector<particles::Code> trackedHadrons = {
-      particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
-      particles::Code::KMinus, particles::Code::K0Long,  particles::Code::K0Short};
-
-  random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
-  process::sibyll::Interaction sibyll;
-  process::sibyll::NuclearInteraction sibyllNuc(sibyll, env);
-  process::sibyll::Decay decay(trackedHadrons);
-  process::particle_cut::ParticleCut cut(20_GeV);
-
-  process::track_writer::TrackWriter trackWriter("tracks.dat");
-  process::energy_loss::EnergyLoss eLoss;
-
-  // assemble all processes into an ordered process list
-  auto sequence = stackInspect << sibyll << sibyllNuc << decay << eLoss << cut
-                               << trackWriter;
-
   // setup particle stack, and add primary particle
   setup::Stack stack;
   stack.Clear();
@@ -148,6 +126,28 @@ int main() {
         beamCode, E0, plab, pos, 0_ns, nuclA, nuclZ});
   }
 
+  // setup processes, decays and interactions
+  tracking_line::TrackingLine tracking;
+  stack_inspector::StackInspector<setup::Stack> stackInspect(1, true, E0);
+
+  const std::vector<particles::Code> trackedHadrons = {
+      particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
+      particles::Code::KMinus, particles::Code::K0Long,  particles::Code::K0Short};
+
+  random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
+  random::RNGManager::GetInstance().RegisterRandomStream("pythia");
+  process::sibyll::Interaction sibyll;
+  process::sibyll::NuclearInteraction sibyllNuc(sibyll, env);
+  process::sibyll::Decay decay(trackedHadrons);
+  process::particle_cut::ParticleCut cut(20_GeV);
+
+  process::track_writer::TrackWriter trackWriter("tracks.dat");
+  process::energy_loss::EnergyLoss eLoss;
+
+  // assemble all processes into an ordered process list
+  auto sequence = stackInspect << sibyll << sibyllNuc << decay << eLoss << cut
+                               << trackWriter;
+
   // define air shower object, run simulation
   cascade::Cascade EAS(env, tracking, sequence, stack);
   EAS.Init();

Documentation/Examples/cascade_proton_example.cc

@@ -11,6 +11,7 @@
 #include <corsika/cascade/Cascade.h>
 #include <corsika/process/ProcessSequence.h>
 #include <corsika/process/hadronic_elastic_model/HadronicElasticModel.h>
+#include <corsika/process/stack_inspector/StackInspector.h>
 #include <corsika/process/tracking_line/TrackingLine.h>
 
 #include <corsika/setup/SetupStack.h>
@@ -85,35 +86,6 @@ int main() {
 
   const CoordinateSystem& rootCS = env.GetCoordinateSystem();
 
-  // setup processes, decays and interactions
-  tracking_line::TrackingLine tracking;
-
-  const std::vector<particles::Code> trackedHadrons = {
-      particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
-      particles::Code::KMinus, particles::Code::K0Long,  particles::Code::K0Short};
-
-  random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
-  random::RNGManager::GetInstance().RegisterRandomStream("pythia");
-  //  process::sibyll::Interaction sibyll(env);
-  process::pythia::Interaction pythia;
-  //  process::sibyll::NuclearInteraction sibyllNuc(env, sibyll);
-  //  process::sibyll::Decay decay(trackedHadrons);
-  process::pythia::Decay decay(trackedHadrons);
-  process::particle_cut::ParticleCut cut(20_GeV);
-
-  // random::RNGManager::GetInstance().RegisterRandomStream("HadronicElasticModel");
-  // process::HadronicElasticModel::HadronicElasticInteraction
-  // hadronicElastic(env);
-
-  process::track_writer::TrackWriter trackWriter("tracks.dat");
-
-  // assemble all processes into an ordered process list
-  // auto sequence = sibyll << decay << hadronicElastic << cut << trackWriter;
-  auto sequence = pythia << decay << cut << trackWriter;
-
-  // cout << "decltype(sequence)=" << type_id_with_cvr<decltype(sequence)>().pretty_name()
-  // << "\n";
-
   // setup particle stack, and add primary particle
   setup::Stack stack;
   stack.Clear();
@@ -145,6 +117,36 @@ int main() {
             beamCode, E0, plab, pos, 0_ns});
   }
 
+  // setup processes, decays and interactions
+  tracking_line::TrackingLine tracking;
+  stack_inspector::StackInspector<setup::Stack> stackInspect(1, true, E0);
+
+  const std::vector<particles::Code> trackedHadrons = {
+      particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
+      particles::Code::KMinus, particles::Code::K0Long,  particles::Code::K0Short};
+
+  random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
+  random::RNGManager::GetInstance().RegisterRandomStream("pythia");
+  //  process::sibyll::Interaction sibyll(env);
+  process::pythia::Interaction pythia;
+  //  process::sibyll::NuclearInteraction sibyllNuc(env, sibyll);
+  //  process::sibyll::Decay decay(trackedHadrons);
+  process::pythia::Decay decay(trackedHadrons);
+  process::particle_cut::ParticleCut cut(20_GeV);
+
+  // random::RNGManager::GetInstance().RegisterRandomStream("HadronicElasticModel");
+  // process::HadronicElasticModel::HadronicElasticInteraction
+  // hadronicElastic(env);
+
+  process::track_writer::TrackWriter trackWriter("tracks.dat");
+
+  // assemble all processes into an ordered process list
+  // auto sequence = sibyll << decay << hadronicElastic << cut << trackWriter;
+  auto sequence = pythia << decay << cut << trackWriter << stackInspect;
+
+  // cout << "decltype(sequence)=" << type_id_with_cvr<decltype(sequence)>().pretty_name()
+  // << "\n";
+
   // define air shower object, run simulation
   cascade::Cascade EAS(env, tracking, sequence, stack);
   EAS.Init();

Documentation/Examples/vertical_EAS.cc

@@ -12,6 +12,7 @@
 #include <corsika/process/ProcessSequence.h>
 #include <corsika/process/energy_loss/EnergyLoss.h>
 #include <corsika/process/hadronic_elastic_model/HadronicElasticModel.h>
+#include <corsika/process/stack_inspector/StackInspector.h>
 #include <corsika/process/tracking_line/TrackingLine.h>
 
 #include <corsika/setup/SetupStack.h>
@@ -27,8 +28,8 @@
 #include <corsika/process/sibyll/Interaction.h>
 #include <corsika/process/sibyll/NuclearInteraction.h>
 
-#include <corsika/process/particle_cut/ParticleCut.h>
 #include <corsika/process/track_writer/TrackWriter.h>
+#include <corsika/process/particle_cut/ParticleCut.h>
 
 #include <corsika/units/PhysicalUnits.h>
 
@@ -86,32 +87,6 @@ int main() {
 
   const CoordinateSystem& rootCS = env.GetCoordinateSystem();
 
-  // setup processes, decays and interactions
-  tracking_line::TrackingLine tracking;
-
-  const std::vector<particles::Code> trackedHadrons = {
-      particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
-      particles::Code::KMinus, particles::Code::K0Long,  particles::Code::K0Short};
-
-  random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
-  process::sibyll::Interaction sibyll;
-  process::sibyll::NuclearInteraction sibyllNuc(sibyll, env);
-  process::sibyll::Decay decay(trackedHadrons);
-  process::particle_cut::ParticleCut cut(20_GeV);
-  // random::RNGManager::GetInstance().RegisterRandomStream("HadronicElasticModel");
-  // process::HadronicElasticModel::HadronicElasticInteraction
-  // hadronicElastic(env);
-
-  process::track_writer::TrackWriter trackWriter("tracks.dat");
-  process::energy_loss::EnergyLoss eLoss;
-
-  // assemble all processes into an ordered process list
-  // cut << trackWriter;
-  auto sequence = sibyll << sibyllNuc << decay << eLoss << cut;
-
-  // cout << "decltype(sequence)=" << type_id_with_cvr<decltype(sequence)>().pretty_name()
-  // << "\n";
-
   // setup particle stack, and add primary particle
   setup::Stack stack;
   stack.Clear();
@@ -146,6 +121,26 @@ int main() {
         beamCode, E0, plab, pos, 0_ns, nuclA, nuclZ});
   }
 
+  // setup processes, decays and interactions
+  tracking_line::TrackingLine tracking;
+  stack_inspector::StackInspector<setup::Stack> stackInspect(1, true, E0);
+
+  const std::vector<particles::Code> trackedHadrons = {
+      particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
+      particles::Code::KMinus, particles::Code::K0Long,  particles::Code::K0Short};
+
+  random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
+  process::sibyll::Interaction sibyll;
+  process::sibyll::NuclearInteraction sibyllNuc(sibyll, env);
+  process::sibyll::Decay decay(trackedHadrons);
+  process::particle_cut::ParticleCut cut(20_GeV);
+
+  process::track_writer::TrackWriter trackWriter("tracks.dat");
+  process::energy_loss::EnergyLoss eLoss;
+
+  // assemble all processes into an ordered process list
+  auto sequence = sibyll << sibyllNuc << decay << eLoss << cut << stackInspect;
+
   // define air shower object, run simulation
   cascade::Cascade EAS(env, tracking, sequence, stack);
   EAS.Init();

Framework/Cascade/testCascade.cc

@@ -131,13 +131,16 @@ public:
 };
 
 TEST_CASE("Cascade", "[Cascade]") {
+
+  HEPEnergyType E0 = 100_GeV;
+
   random::RNGManager& rmng = random::RNGManager::GetInstance();
   rmng.RegisterRandomStream("cascade");
 
   auto env = MakeDummyEnv();
   tracking_line::TrackingLine tracking;
 
-  stack_inspector::StackInspector<TestCascadeStack> stackInspect(1, true);
+  stack_inspector::StackInspector<TestCascadeStack> stackInspect(1, true, E0);
   null_model::NullModel nullModel;
 
   const GrammageType X0 = 20_g / square(1_cm);
@@ -154,7 +157,6 @@ TEST_CASE("Cascade", "[Cascade]") {
       RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
 
   stack.Clear();
-  HEPEnergyType E0 = 100_GeV;
   stack.AddParticle(
       std::tuple<particles::Code, units::si::HEPEnergyType,
                  corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{

Framework/ProcessSequence/ProcessSequence.h

@@ -136,6 +136,14 @@ namespace corsika::process {
       return ret;
     }
 
+    /**
+       The processes of type StackProcess do have an internal counter,
+       so they can be exectuted only each N steps. Often these are
+       "maintenacne processes" that do not need to run after each
+       single step of the simulations. In the CheckStep function it is
+       tested if either A or B are StackProcess and if they are due
+       for execution.
+     */
     bool CheckStep() {
       bool ret = false;
       if constexpr (std::is_base_of_v<StackProcess<T1type>, T1type> || t1ProcSeq) {
@@ -147,6 +155,9 @@ namespace corsika::process {
       return ret;
     }
 
+    /**
+       Execute the StackProcess-es in the ProcessSequence
+     */
     template <typename TStack>
     EProcessReturn DoStack(TStack& vS) {
       EProcessReturn ret = EProcessReturn::eOk;

Framework/ProcessSequence/StackProcess.h

@@ -42,6 +42,7 @@ namespace corsika::process {
     template <typename TStack>
     inline EProcessReturn DoStack(TStack&);
 
+    int GetStep() const { return fIStep; }
     bool CheckStep() { return !((++fIStep) % fNStep); }
 
   private:

Processes/CMakeLists.txt

+# general
 add_subdirectory (NullModel)
+# tracking
+add_subdirectory (TrackingLine) 
+# hadron interaction models
 add_subdirectory (Sibyll)
 if (PYTHIA8_FOUND)
   add_subdirectory (Pythia)
 endif (PYTHIA8_FOUND)
-add_subdirectory (StackInspector)
-add_subdirectory (TrackWriter)
-add_subdirectory (TrackingLine)
 add_subdirectory (HadronicElasticModel)
-add_subdirectory (EnergyLoss)
 add_subdirectory (UrQMD)
 add_subdirectory (SwitchProcess)
+
+# continuous physics
+add_subdirectory (EnergyLoss)
+add_subdirectory (TrackWriter)
+# stack processes
+add_subdirectory (StackInspector)
+# secondaries process
+# cuts, thinning, etc.
 add_subdirectory (ParticleCut)
 
+##########################################
+# add_custom_target(CORSIKAprocesses)
 add_library (CORSIKAprocesses INTERFACE)
 add_dependencies(CORSIKAprocesses ProcessNullModel)
 add_dependencies(CORSIKAprocesses ProcessSibyll)

Processes/StackInspector/StackInspector.cc

@@ -17,6 +17,8 @@
 
 #include <corsika/setup/SetupTrajectory.h>
 
+#include <chrono>
+#include <iomanip>
 #include <iostream>
 #include <limits>
 using namespace std;
@@ -27,41 +29,58 @@ using namespace corsika::units::si;
 using namespace corsika::process::stack_inspector;
 
 template <typename TStack>
-StackInspector<TStack>::StackInspector(const int nStep, const bool aReport)
-    : StackProcess<StackInspector<TStack>>(nStep)
-    , fReport(aReport)
-    , fCountStep(0) {}
+StackInspector<TStack>::StackInspector(const int vNStep, const bool vReportStack,
+                                       const HEPEnergyType vE0)
+    : StackProcess<StackInspector<TStack>>(vNStep)
+    , fReportStack(vReportStack)
+    , fE0(vE0)
+    , fStartTime(std::chrono::system_clock::now()) {}
 
 template <typename TStack>
 StackInspector<TStack>::~StackInspector() {}
 
 template <typename TStack>
-process::EProcessReturn StackInspector<TStack>::DoStack(TStack const& vS) {
-  if (!fReport) return process::EProcessReturn::eOk;
+process::EProcessReturn StackInspector<TStack>::DoStack(const TStack& vS) {
   [[maybe_unused]] int i = 0;
   HEPEnergyType Etot = 0_GeV;
 
-  for (auto& iterP : vS) {
+  for (const auto& iterP : vS) {
     HEPEnergyType E = iterP.GetEnergy();
     Etot += E;
-    geometry::CoordinateSystem& rootCS = geometry::RootCoordinateSystem::GetInstance()
-                                             .GetRootCoordinateSystem(); // for printout
-    auto pos = iterP.GetPosition().GetCoordinates(rootCS);
-    cout << "StackInspector: i=" << setw(5) << fixed << (i++) << ", id=" << setw(30)
-         << iterP.GetPID() << " E=" << setw(15) << scientific << (E / 1_GeV) << " GeV, "
-         << " pos=" << pos << " node = " << iterP.GetNode();
-    if (iterP.GetPID() == Code::Nucleus) cout << " nuc_ref=" << iterP.GetNucleusRef();
-    cout << endl;
+    if (fReportStack) {
+      geometry::CoordinateSystem& rootCS = geometry::RootCoordinateSystem::GetInstance()
+                                               .GetRootCoordinateSystem(); // for printout
+      auto pos = iterP.GetPosition().GetCoordinates(rootCS);
+      cout << "StackInspector: i=" << setw(5) << fixed << (i++) << ", id=" << setw(30)
+           << iterP.GetPID() << " E=" << setw(15) << scientific << (E / 1_GeV) << " GeV, "
+           << " pos=" << pos << " node = " << iterP.GetNode();
+      if (iterP.GetPID() == Code::Nucleus) cout << " nuc_ref=" << iterP.GetNucleusRef();
+      cout << endl;
+    }
   }
-  fCountStep++;
-  cout << "StackInspector: nStep=" << fCountStep << " stackSize=" << vS.GetSize()
-       << " Estack=" << Etot / 1_GeV << " GeV" << endl;
+
+  auto const now = std::chrono::system_clock::now();
+  const std::chrono::duration<double> elapsed_seconds = now - fStartTime;
+  std::time_t const now_time = std::chrono::system_clock::to_time_t(now);
+  double const progress = (fE0 - Etot) / fE0;
+  double const eta_seconds = elapsed_seconds.count() / progress;
+  std::time_t const eta_time = std::chrono::system_clock::to_time_t(
+      fStartTime + std::chrono::seconds((int)eta_seconds));
+
+  cout << "StackInspector: "
+       << " time=" << std::put_time(std::localtime(&now_time), "%T")
+       << ", running=" << elapsed_seconds.count() << " seconds"
+       << " (" << setw(3) << int(progress * 100) << "%)"
+       << ", nStep=" << GetStep() << ", stackSize=" << vS.GetSize()
+       << ", Estack=" << Etot / 1_GeV << " GeV"
+       << ", ETA=" << std::put_time(std::localtime(&eta_time), "%T") << endl;
   return process::EProcessReturn::eOk;
 }
 
 template <typename TStack>
 void StackInspector<TStack>::Init() {
-  fCountStep = 0;
+  fReportStack = false;
+  fStartTime = std::chrono::system_clock::now();
 }
 
 #include <corsika/cascade/testCascade.h>

Processes/StackInspector/StackInspector.h

@@ -15,6 +15,8 @@
 #include <corsika/setup/SetupTrajectory.h>
 #include <corsika/units/PhysicalUnits.h>
 
+#include <chrono>
+
 namespace corsika::process {
 
   namespace stack_inspector {
@@ -24,16 +26,25 @@ namespace corsika::process {
 
       typedef typename TStack::ParticleType Particle;
 
+      using corsika::process::StackProcess<StackInspector<TStack>>::GetStep;
+
     public:
-      StackInspector(const int nStep, const bool aReport);
+      StackInspector(const int vNStep, const bool vReportStack,
+                     const corsika::units::si::HEPEnergyType vE0);
       ~StackInspector();
 
       void Init();
-      EProcessReturn DoStack(TStack const&);
+      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) { fE0 = vE0; }
 
     private:
-      bool fReport;
-      int fCountStep = 0;
+      bool fReportStack;
+      corsika::units::si::HEPEnergyType fE0;
+      decltype(std::chrono::system_clock::now()) fStartTime;
     };
 
   } // namespace stack_inspector

Processes/StackInspector/testStackInspector.cc

@@ -49,7 +49,7 @@ TEST_CASE("StackInspector", "[processes]") {
 
   SECTION("interface") {
 
-    StackInspector<TestCascadeStack> model(1, true);
+    StackInspector<TestCascadeStack> model(1, true, E0);
 
     model.Init();
     [[maybe_unused]] const process::EProcessReturn ret = model.DoStack(stack);

Processes/SwitchProcess/CMakeLists.txt

@@ -11,14 +11,6 @@ set (
 add_library (ProcessSwitch INTERFACE)
 CORSIKA_COPY_HEADERS_TO_NAMESPACE (ProcessSwitch ${MODEL_NAMESPACE} ${MODEL_HEADERS})
 
-set_target_properties (
-  ProcessStackInspector
-  PROPERTIES
-  VERSION ${PROJECT_VERSION}
-  SOVERSION 1
-#  PUBLIC_HEADER "${MODEL_HEADERS}"
-  )
-
 # target dependencies on other libraries (also the header onlys)
 target_link_libraries (
   ProcessSwitch