From b6a749280bcf89d8f6515136ae59eb955e146378 Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Mon, 12 Oct 2020 21:11:26 +0200
Subject: [PATCH] make examples run MUCH faster and with less output

---
 Documentation/Examples/boundary_example.cc    | 32 +++++++----
 Documentation/Examples/cascade_example.cc     |  3 +
 .../Examples/cascade_proton_example.cc        |  4 ++
 Documentation/Examples/em_shower.cc           |  5 ++
 Documentation/Examples/vertical_EAS.cc        |  4 +-
 Framework/Utilities/CMakeLists.txt            |  1 +
 Framework/Utilities/COMBoost.cc               |  6 +-
 Framework/Utilities/COMBoost.h                | 18 +++---
 Processes/EnergyLoss/CMakeLists.txt           |  1 +
 Processes/EnergyLoss/EnergyLoss.cc            | 32 +++++------
 Processes/Pythia/Decay.cc                     |  9 ++-
 Processes/Pythia/Decay.h                      |  3 +-
 Processes/Pythia/Interaction.cc               | 11 ++--
 Processes/Pythia/Interaction.h                |  5 +-
 Processes/Sibyll/Decay.cc                     | 55 ++++++++++---------
 Processes/Sibyll/Decay.h                      |  5 +-
 Processes/Sibyll/Interaction.cc               | 13 +++--
 Processes/Sibyll/Interaction.h                |  3 +-
 Processes/TrackingLine/CMakeLists.txt         |  2 +-
 Processes/TrackingLine/TrackingLine.cc        |  1 +
 Processes/TrackingLine/TrackingLine.h         | 38 +++++--------
 21 files changed, 139 insertions(+), 112 deletions(-)

diff --git a/Documentation/Examples/boundary_example.cc b/Documentation/Examples/boundary_example.cc
index 5d05e2bac..a5279e89e 100644
--- a/Documentation/Examples/boundary_example.cc
+++ b/Documentation/Examples/boundary_example.cc
@@ -34,6 +34,8 @@
 
 #include <corsika/utl/CorsikaFenv.h>
 
+#include <corsika/logging/Logging.h>
+
 #include <iostream>
 #include <limits>
 #include <typeinfo>
@@ -60,7 +62,9 @@ struct MyBoundaryCrossingProcess
   EProcessReturn DoBoundaryCrossing(Particle& p,
                                     typename Particle::BaseNodeType const& from,
                                     typename Particle::BaseNodeType const& to) {
-    std::cout << "boundary crossing! from: " << &from << "; to: " << &to << std::endl;
+
+    C8LOG_INFO("MyBoundaryCrossingProcess: crossing! from: {} to: {} ", fmt::ptr(&from),
+               fmt::ptr(&to));
 
     auto const& name = particles::GetName(p.GetPID());
     auto const start = p.GetPosition().GetCoordinates();
@@ -82,7 +86,9 @@ private:
 //
 int main() {
 
-  std::cout << "boundary_example" << std::endl;
+  logging::SetLevel(logging::level::info);
+
+  C8LOG_INFO("boundary_example");
 
   feenableexcept(FE_INVALID);
   // initialize random number sequence(s)
@@ -121,7 +127,7 @@ int main() {
   process::sibyll::Interaction sibyll;
   process::sibyll::Decay decay;
 
-  process::particle_cut::ParticleCut cut(20_GeV, true, true);
+  process::particle_cut::ParticleCut cut(50_GeV, true, true);
 
   process::track_writer::TrackWriter trackWriter("tracks.dat");
   MyBoundaryCrossingProcess<true> boundaryCrossing("crossings.dat");
@@ -132,9 +138,9 @@ int main() {
   // setup particle stack, and add primary particles
   setup::Stack stack;
   stack.Clear();
-  const Code beamCode = Code::Proton;
-  const HEPMassType mass = particles::GetMass(Code::Proton);
-  const HEPEnergyType E0 = 50_TeV;
+  const Code beamCode = Code::MuPlus;
+  const HEPMassType mass = particles::GetMass(beamCode);
+  const HEPEnergyType E0 = 100_GeV;
 
   std::uniform_real_distribution distTheta(0., 180.);
   std::uniform_real_distribution distPhi(0., 360.);
@@ -155,9 +161,11 @@ int main() {
     auto const [px, py, pz] =
         momentumComponents(theta / 180. * M_PI, phi / 180. * M_PI, P0);
     auto plab = corsika::stack::MomentumVector(rootCS, {px, py, pz});
-    cout << "input particle: " << beamCode << endl;
-    cout << "input angles: theta=" << theta << " phi=" << phi << endl;
-    cout << "input momentum: " << plab.GetComponents() / 1_GeV << endl;
+    C8LOG_INFO(
+        "input particle: {} "
+        "input angles: theta={} phi={}"
+        "input momentum: {} GeV",
+        beamCode, theta, phi, plab.GetComponents() / 1_GeV);
     Point pos(rootCS, 0_m, 0_m, 0_m);
     stack.AddParticle(
         std::tuple<particles::Code, units::si::HEPEnergyType,
@@ -170,10 +178,10 @@ int main() {
 
   EAS.Run();
 
-  cout << "Result: E0=" << E0 / 1_GeV << endl;
+  C8LOG_INFO("Result: E0={}GeV", E0 / 1_GeV);
   cut.ShowResults();
   const HEPEnergyType Efinal =
       cut.GetCutEnergy() + cut.GetInvEnergy() + cut.GetEmEnergy();
-  cout << "total energy (GeV): " << Efinal / 1_GeV << endl
-       << "relative difference (%): " << (Efinal / E0 - 1.) * 100 << endl;
+  C8LOG_INFO("Total energy (GeV): {} relative difference (%): {}", Efinal / 1_GeV,
+             (Efinal / E0 - 1.) * 100);
 }
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index 993714c08..45e26dc85 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -35,6 +35,7 @@
 #include <corsika/random/RNGManager.h>
 
 #include <corsika/utl/CorsikaFenv.h>
+#include <corsika/logging/Logging.h>
 
 #include <iostream>
 #include <limits>
@@ -56,6 +57,8 @@ using namespace corsika::units::si;
 //
 int main() {
 
+  logging::SetLevel(logging::level::info);
+
   std::cout << "cascade_example" << std::endl;
 
   const LengthType height_atmosphere = 112.8_km;
diff --git a/Documentation/Examples/cascade_proton_example.cc b/Documentation/Examples/cascade_proton_example.cc
index b7d3970c3..235ee2d9e 100644
--- a/Documentation/Examples/cascade_proton_example.cc
+++ b/Documentation/Examples/cascade_proton_example.cc
@@ -38,6 +38,8 @@
 
 #include <corsika/utl/CorsikaFenv.h>
 
+#include <corsika/logging/Logging.h>
+
 #include <iostream>
 #include <limits>
 #include <typeinfo>
@@ -59,6 +61,8 @@ using namespace corsika::units::si;
 //
 int main() {
 
+  logging::SetLevel(logging::level::info);
+
   std::cout << "cascade_proton_example" << std::endl;
 
   feenableexcept(FE_INVALID);
diff --git a/Documentation/Examples/em_shower.cc b/Documentation/Examples/em_shower.cc
index 5a528975a..95dfedfc4 100644
--- a/Documentation/Examples/em_shower.cc
+++ b/Documentation/Examples/em_shower.cc
@@ -29,6 +29,8 @@
 #include <corsika/utl/CorsikaFenv.h>
 #include <corsika/process/interaction_counter/InteractionCounter.hpp>
 
+#include <corsika/logging/Logging.h>
+
 #include <iomanip>
 #include <iostream>
 #include <limits>
@@ -54,6 +56,9 @@ void registerRandomStreams() {
 }
 
 int main(int argc, char** argv) {
+
+  logging::SetLevel(logging::level::info);
+
   if (argc != 2) {
     std::cerr << "usage: em_shower <energy/GeV>" << std::endl;
     return 1;
diff --git a/Documentation/Examples/vertical_EAS.cc b/Documentation/Examples/vertical_EAS.cc
index 2e37e11a5..9b0cf6182 100644
--- a/Documentation/Examples/vertical_EAS.cc
+++ b/Documentation/Examples/vertical_EAS.cc
@@ -78,7 +78,7 @@ void registerRandomStreams(const int seed) {
 
 int main(int argc, char** argv) {
 
-  logging::SetLevel(logging::level::trace);
+  logging::SetLevel(logging::level::info);
 
   C8LOG_INFO("vertical_EAS");
 
@@ -199,7 +199,7 @@ int main(int argc, char** argv) {
 
   decaySibyll.PrintDecayConfig();
 
-  process::particle_cut::ParticleCut cut{50_GeV, false, true};
+  process::particle_cut::ParticleCut cut{60_GeV, false, true};
   process::proposal::Interaction proposal(env, cut.GetECut());
   process::proposal::ContinuousProcess em_continuous(env, cut.GetECut());
   process::interaction_counter::InteractionCounter proposalCounted(proposal);
diff --git a/Framework/Utilities/CMakeLists.txt b/Framework/Utilities/CMakeLists.txt
index 5618ab4a8..95170d4a3 100644
--- a/Framework/Utilities/CMakeLists.txt
+++ b/Framework/Utilities/CMakeLists.txt
@@ -37,6 +37,7 @@ set (
   UTILITIES_DEPENDS
   CORSIKAgeometry
   CORSIKAunits
+  CORSIKAlogging
   C8::ext::boost # so far only for MetaProgramming
   C8::ext::eigen3 # for COMboost
   )
diff --git a/Framework/Utilities/COMBoost.cc b/Framework/Utilities/COMBoost.cc
index c6688199a..5bb4be9da 100644
--- a/Framework/Utilities/COMBoost.cc
+++ b/Framework/Utilities/COMBoost.cc
@@ -12,6 +12,7 @@
 #include <corsika/units/PhysicalUnits.h>
 #include <corsika/utl/COMBoost.h>
 #include <corsika/utl/sgn.h>
+#include <corsika/logging/Logging.h>
 
 #include <cmath>
 
@@ -38,9 +39,8 @@ COMBoost::COMBoost(FourVector<HEPEnergyType, Vector<hepmomentum_d>> const& Pproj
 
   setBoost(coshEta, sinhEta);
 
-  std::cout << "COMBoost (1-beta)=" << 1 - sinhEta / coshEta << " gamma=" << coshEta
-            << std::endl;
-  std::cout << "  det = " << boost_.determinant() - 1 << std::endl;
+  C8LOG_TRACE("COMBoost (1-beta)={}, gamma={}, det={}", 1 - sinhEta / coshEta, coshEta,
+              boost_.determinant() - 1);
 }
 
 COMBoost::COMBoost(geometry::Vector<units::si::hepmomentum_d> const& momentum,
diff --git a/Framework/Utilities/COMBoost.h b/Framework/Utilities/COMBoost.h
index 0c78d49ac..73b7c9300 100644
--- a/Framework/Utilities/COMBoost.h
+++ b/Framework/Utilities/COMBoost.h
@@ -11,6 +11,7 @@
 #include <corsika/geometry/CoordinateSystem.h>
 #include <corsika/geometry/FourVector.h>
 #include <corsika/units/PhysicalUnits.h>
+#include <corsika/logging/Logging.h>
 
 #include <Eigen/Dense>
 
@@ -69,10 +70,10 @@ namespace corsika::utl {
       Eigen::Vector2d com;
       com << (Ecm * (1 / 1_GeV)), (pCM.eVector(2) * (1 / 1_GeV).magnitude());
 
-      std::cout << "COMBoost::fromCoM Ecm=" << Ecm / 1_GeV << " GeV, "
-                << " pcm = " << pCM / 1_GeV << " (norm = " << pCM.norm() / 1_GeV
-                << " GeV), invariant mass = " << p.GetNorm() / 1_GeV << " GeV"
-                << std::endl;
+      C8LOG_TRACE(
+          "COMBoost::fromCoM Ecm={} GeV"
+          " pcm={} GeV (norm = {} GeV), invariant mass={} GeV",
+          Ecm / 1_GeV, pCM / 1_GeV, pCM.norm() / 1_GeV, p.GetNorm() / 1_GeV);
 
       auto const boostedZ = inverseBoost_ * com;
       auto const E_lab = boostedZ(0) * 1_GeV;
@@ -84,10 +85,11 @@ namespace corsika::utl {
 
       FourVector f(E_lab, pLab);
 
-      std::cout << "COMBoost::fromCoM --> Elab=" << E_lab / 1_GeV << "GeV, "
-                << " plab = " << pLab.GetComponents() << " (norm =" << pLab.norm() / 1_GeV
-                << " GeV), invariant mass = " << f.GetNorm() / 1_GeV << " GeV"
-                << std::endl;
+      C8LOG_TRACE("COMBoost::fromCoM --> Elab={} GeV",
+                  " plab={} GeV (norm={} GeV) "
+                  " GeV), invariant mass = {}",
+                  E_lab / 1_GeV, f.GetNorm() / 1_GeV, pLab.GetComponents(),
+                  pLab.norm() / 1_GeV);
 
       return f;
     }
diff --git a/Processes/EnergyLoss/CMakeLists.txt b/Processes/EnergyLoss/CMakeLists.txt
index 62fca1ca0..3df3978d5 100644
--- a/Processes/EnergyLoss/CMakeLists.txt
+++ b/Processes/EnergyLoss/CMakeLists.txt
@@ -32,6 +32,7 @@ target_link_libraries (
   CORSIKAgeometry
   CORSIKAenvironment
   CORSIKAsetup
+  CORSIKAlogging
   )
 
 target_include_directories (
diff --git a/Processes/EnergyLoss/EnergyLoss.cc b/Processes/EnergyLoss/EnergyLoss.cc
index 107fa813a..42efc08a5 100644
--- a/Processes/EnergyLoss/EnergyLoss.cc
+++ b/Processes/EnergyLoss/EnergyLoss.cc
@@ -13,6 +13,8 @@
 #include <corsika/setup/SetupStack.h>
 #include <corsika/setup/SetupTrajectory.h>
 
+#include <corsika/logging/Logging.h>
+
 #include <corsika/geometry/Line.h>
 
 #include <cmath>
@@ -85,19 +87,18 @@ HEPEnergyType EnergyLoss::BetheBloch(SetupParticle const& p, GrammageType const
   double const beta2 = (gamma2 - 1) / gamma2; // 1-1/gamma2    (1-1/gamma)*(1+1/gamma);
                                               // (gamma_2-1)/gamma_2 = (1-1/gamma2);
   double constexpr c2 = 1;                    // HEP convention here c=c2=1
-  cout << "BetheBloch beta2=" << beta2 << " gamma2=" << gamma2 << endl;
+  C8LOG_DEBUG("BetheBloch beta2={}, gamma2={}", beta2, gamma2);
   [[maybe_unused]] double const eta2 = beta2 / (1 - beta2);
   HEPMassType const Wmax =
       2 * me * c2 * beta2 * gamma2 / (1 + 2 * gamma * me / m + me2 / m2);
   // approx, but <<1%    HEPMassType const Wmax = 2*me*c2*beta2*gamma2;      for HEAVY
   // PARTICLES Wmax ~ 2me v2 for non-relativistic particles
-  cout << "BetheBloch Wmax=" << Wmax << endl;
+  C8LOG_DEBUG("BetheBloch Wmax={}", Wmax);
 
   // Sternheimer parameterization, density corrections towards high energies
   // NOTE/TODO: when Cbar is 0 it needs to be approximated from parameterization ->
   // MISSING
-  cout << "BetheBloch p.GetMomentum().GetNorm()/m=" << p.GetMomentum().GetNorm() / m
-       << endl;
+  C8LOG_DEBUG("BetheBloch p.GetMomentum().GetNorm()/m{}=", p.GetMomentum().GetNorm() / m);
   double const x = log10(p.GetMomentum().GetNorm() / m);
   double delta = 0;
   if (x >= x1) {
@@ -107,7 +108,7 @@ HEPEnergyType EnergyLoss::BetheBloch(SetupParticle const& p, GrammageType const
   } else if (x < x0) { // and IF conductor (otherwise, this is 0)
     delta = delta0 * pow(100, 2 * (x - x0));
   }
-  cout << "BetheBloch delta=" << delta << endl;
+  C8LOG_DEBUG("BetheBloch delta={}", delta);
 
   // with further low energies correction, accurary ~1% down to beta~0.05 (1MeV for p)
 
@@ -141,9 +142,6 @@ HEPEnergyType EnergyLoss::BetheBloch(SetupParticle const& p, GrammageType const
   double const y2 = Z * Z * alpha * alpha / beta2;
   double const bloch = -y2 * (1.202 - y2 * (1.042 - 0.855 * y2 + 0.343 * y2 * y2));
 
-  // cout << "BetheBloch Erel=" << Erel << " barkas=" << barkas << " bloch=" << bloch <<
-  // endl;
-
   double const aux = 2 * me * c2 * beta2 * gamma2 * Wmax / (Ieff * Ieff);
   return -K * Z2 * ZoverA / beta2 *
          (0.5 * log(aux) - beta2 - Cadj / Z - delta / 2 + barkas + bloch) * dX;
@@ -168,15 +166,14 @@ process::EProcessReturn EnergyLoss::DoContinuous(SetupParticle& p, SetupTrack co
 
   GrammageType const dX =
       p.GetNode()->GetModelProperties().IntegratedGrammage(t, t.GetLength());
-  cout << "EnergyLoss " << p.GetPID() << ", z=" << p.GetChargeNumber()
-       << ", dX=" << dX / 1_g * square(1_cm) << "g/cm2" << endl;
+  C8LOG_DEBUG("EnergyLoss pid={}, z={}, dX={} g/cm2", p.GetPID(), p.GetChargeNumber(),
+              dX / 1_g * square(1_cm));
   HEPEnergyType dE = TotalEnergyLoss(p, dX);
   auto E = p.GetEnergy();
   const auto Ekin = E - p.GetMass();
   auto Enew = E + dE;
-  cout << "EnergyLoss  dE=" << dE / 1_MeV << "MeV, "
-       << " E=" << E / 1_GeV << "GeV,  Ekin=" << Ekin / 1_GeV << ", Enew=" << Enew / 1_GeV
-       << "GeV" << endl;
+  C8LOG_DEBUG("EnergyLoss  dE={} MeV, E={} GeV, Ekin={} GeV, Enew={} GeV", dE / 1_MeV,
+              E / 1_GeV, Ekin / 1_GeV, Enew / 1_GeV);
   p.SetEnergy(Enew);
   MomentumUpdate(p, Enew);
   FillProfile(t, dE);
@@ -226,8 +223,8 @@ void EnergyLoss::FillProfile(SetupTrack const& vTrack, const HEPEnergyType dE) {
   const int binStart = grammageStart / dX_;
   const int binEnd = grammageEnd / dX_;
 
-  std::cout << "energy deposit of " << -dE << " between " << grammageStart << " and "
-            << grammageEnd << std::endl;
+  C8LOG_DEBUG("energy deposit of -dE={} between {} and {}", -dE, grammageStart,
+              grammageEnd);
 
   auto energyCount = HEPEnergyType::zero();
 
@@ -237,8 +234,7 @@ void EnergyLoss::FillProfile(SetupTrack const& vTrack, const HEPEnergyType dE) {
       profile_[bin] += increment;
       energyCount += increment;
 
-      std::cout << "filling bin " << bin << " with weight " << weight << ": " << increment
-                << std::endl;
+      C8LOG_DEBUG("filling bin {} with weight {} : {} ", bin, weight, increment);
     }
   };
 
@@ -250,7 +246,7 @@ void EnergyLoss::FillProfile(SetupTrack const& vTrack, const HEPEnergyType dE) {
 
   for (int bin = binStart + 1; bin < binEnd; ++bin) { fill(bin, dX_); }
 
-  std::cout << "total energy added to histogram: " << energyCount << std::endl;
+  C8LOG_DEBUG("total energy added to histogram: {} ", energyCount);
 }
 
 void EnergyLoss::PrintProfile() const {
diff --git a/Processes/Pythia/Decay.cc b/Processes/Pythia/Decay.cc
index 4caf02487..92920dfba 100644
--- a/Processes/Pythia/Decay.cc
+++ b/Processes/Pythia/Decay.cc
@@ -28,7 +28,8 @@ using Track = Trajectory;
 
 namespace corsika::process::pythia {
 
-  Decay::Decay() {
+  Decay::Decay(const bool print_listing)
+      : print_listing_(print_listing) {
 
     // set random number generator in pythia
     Pythia8::RndmEngine* rndm = new corsika::process::pythia::Random();
@@ -231,8 +232,10 @@ namespace corsika::process::pythia {
     else
       cout << "Pythia::Decay: particles after decay: " << event.size() << endl;
 
-    // list final state
-    event.list();
+    if (print_listing_) {
+      // list final state
+      event.list();
+    }
 
     // loop over final state
     for (int i = 0; i < event.size(); ++i)
diff --git a/Processes/Pythia/Decay.h b/Processes/Pythia/Decay.h
index 733b9eae1..51a43b7e3 100644
--- a/Processes/Pythia/Decay.h
+++ b/Processes/Pythia/Decay.h
@@ -23,9 +23,10 @@ namespace corsika::process {
     class Decay : public corsika::process::DecayProcess<Decay> {
       int fCount = 0;
       bool handleAllDecays_ = true;
+      bool print_listing_ = false;
 
     public:
-      Decay();
+      Decay(const bool print_listing = false);
       Decay(std::set<particles::Code>);
       ~Decay();
 
diff --git a/Processes/Pythia/Interaction.cc b/Processes/Pythia/Interaction.cc
index da9ef31e9..7a0dbfa64 100644
--- a/Processes/Pythia/Interaction.cc
+++ b/Processes/Pythia/Interaction.cc
@@ -29,9 +29,9 @@ namespace corsika::process::pythia {
 
   typedef corsika::geometry::Vector<corsika::units::si::hepmomentum_d> MomentumVector;
 
-  Interaction::~Interaction() {}
+  Interaction::Interaction(const bool print_listing)
+      : print_listing_(print_listing) {
 
-  Interaction::Interaction() {
     cout << "Pythia::Interaction n=" << fCount << endl;
 
     using random::RNGManager;
@@ -365,8 +365,11 @@ namespace corsika::process::pythia {
 
         // link to pythia stack
         Pythia8::Event& event = fPythia.event;
-        // print final state
-        event.list();
+
+        if (print_listing_) {
+          // list final state
+          event.list();
+        }
 
         MomentumVector Plab_final(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
         HEPEnergyType Elab_final = 0_GeV;
diff --git a/Processes/Pythia/Interaction.h b/Processes/Pythia/Interaction.h
index 67a6ce423..ef2bdfbcc 100644
--- a/Processes/Pythia/Interaction.h
+++ b/Processes/Pythia/Interaction.h
@@ -22,10 +22,11 @@ namespace corsika::process::pythia {
 
     int fCount = 0;
     bool fInitialized = false;
+    bool print_listing_ = false;
 
   public:
-    Interaction();
-    ~Interaction();
+    Interaction(const bool print_listing = false);
+    ~Interaction() = default;
 
     void SetParticleListStable(std::vector<particles::Code> const&);
     void SetUnstable(const corsika::particles::Code);
diff --git a/Processes/Sibyll/Decay.cc b/Processes/Sibyll/Decay.cc
index 330dcfd27..a85541adc 100644
--- a/Processes/Sibyll/Decay.cc
+++ b/Processes/Sibyll/Decay.cc
@@ -27,7 +27,8 @@ using Particle = corsika::setup::Stack::ParticleType;
 
 namespace corsika::process::sibyll {
 
-  Decay::Decay() {
+  Decay::Decay(const bool sibyll_printout_on)
+      : sibyll_listing_(sibyll_printout_on) {
     // switch off decays to avoid internal decay chains
     SetAllStable();
     // handle all decays by default
@@ -40,7 +41,7 @@ namespace corsika::process::sibyll {
     SetAllStable();
   }
 
-  Decay::~Decay() { C8LOG_DEBUG(fmt::format("Sibyll::Decay n={}", fCount)); }
+  Decay::~Decay() { C8LOG_DEBUG("Sibyll::Decay n={}", count_); }
 
   bool Decay::CanHandleDecay(const particles::Code vParticleCode) {
     using namespace corsika::particles;
@@ -60,7 +61,7 @@ namespace corsika::process::sibyll {
 
   void Decay::SetHandleDecay(const particles::Code vParticleCode) {
     handleAllDecays_ = false;
-    C8LOG_DEBUG(fmt::format("Sibyll::Decay: set to handle decay of {}", vParticleCode));
+    C8LOG_DEBUG("Sibyll::Decay: set to handle decay of {}", vParticleCode);
     if (Decay::CanHandleDecay(vParticleCode))
       handledDecays_.insert(vParticleCode);
     else
@@ -102,13 +103,13 @@ namespace corsika::process::sibyll {
   }
 
   void Decay::SetUnstable(const particles::Code vCode) {
-    C8LOG_DEBUG(fmt::format("Sibyll::Decay: setting {} unstable. ", vCode));
+    C8LOG_DEBUG("Sibyll::Decay: setting {} unstable. ", vCode);
     const int s_id = abs(process::sibyll::ConvertToSibyllRaw(vCode));
     s_csydec_.idb[s_id - 1] = abs(s_csydec_.idb[s_id - 1]);
   }
 
   void Decay::SetStable(const particles::Code vCode) {
-    C8LOG_DEBUG(fmt::format("Sibyll::Decay: setting {} stable. ", vCode));
+    C8LOG_DEBUG("Sibyll::Decay: setting {} stable. ", vCode);
     const int s_id = abs(process::sibyll::ConvertToSibyllRaw(vCode));
     s_csydec_.idb[s_id - 1] = (-1) * abs(s_csydec_.idb[s_id - 1]);
   }
@@ -124,19 +125,17 @@ namespace corsika::process::sibyll {
   void Decay::PrintDecayConfig(const particles::Code vCode) {
     const int sibCode = process::sibyll::ConvertToSibyllRaw(vCode);
     const int absSibCode = abs(sibCode);
-    C8LOG_DEBUG(
-        fmt::format("Decay: Sibyll decay configuration: {} is {}", vCode,
-                    (s_csydec_.idb[absSibCode - 1] <= 0) ? "stable" : "unstable"));
+    C8LOG_DEBUG("Decay: Sibyll decay configuration: {} is {}", vCode,
+                (s_csydec_.idb[absSibCode - 1] <= 0) ? "stable" : "unstable");
   }
 
   void Decay::PrintDecayConfig() {
-    C8LOG_DEBUG(fmt::format("Sibyll::Decay: decay configuration:"));
+    C8LOG_DEBUG("Sibyll::Decay: decay configuration:");
     if (handleAllDecays_) {
-      C8LOG_DEBUG(fmt::format(
-          "     all particles known to Sibyll are handled by Sibyll::Decay!"));
+      C8LOG_DEBUG("     all particles known to Sibyll are handled by Sibyll::Decay!");
     } else {
       for (auto& pCode : handledDecays_) {
-        C8LOG_DEBUG(fmt::format("      Decay of {}  is handled by Sibyll!", pCode));
+        C8LOG_DEBUG("      Decay of {}  is handled by Sibyll!", pCode);
       }
     }
   }
@@ -157,17 +156,17 @@ namespace corsika::process::sibyll {
 
       const auto mkin =
           (E * E - vP.GetMomentum().squaredNorm()); // delta_mass(vP.GetMomentum(), E, m);
-      C8LOG_DEBUG(fmt::format("Sibyll::Decay: code: {} ", vP.GetPID()));
-      C8LOG_DEBUG(fmt::format("Sibyll::Decay: MinStep: t0: {} ", t0));
-      C8LOG_DEBUG(fmt::format("Sibyll::Decay: MinStep: energy: {} GeV ", E / 1_GeV));
-      C8LOG_DEBUG(fmt::format("Sibyll::Decay: momentum: {} GeV ",
-                              vP.GetMomentum().GetComponents() / 1_GeV));
-      C8LOG_DEBUG(fmt::format("Sibyll::Decay: momentum: shell mass-kin. inv. mass {} {}",
-                              mkin / 1_GeV / 1_GeV, m / 1_GeV * m / 1_GeV));
+      C8LOG_DEBUG("Sibyll::Decay: code: {} ", vP.GetPID());
+      C8LOG_DEBUG("Sibyll::Decay: MinStep: t0: {} ", t0);
+      C8LOG_DEBUG("Sibyll::Decay: MinStep: energy: {} GeV ", E / 1_GeV);
+      C8LOG_DEBUG("Sibyll::Decay: momentum: {} GeV ",
+                  vP.GetMomentum().GetComponents() / 1_GeV);
+      C8LOG_DEBUG("Sibyll::Decay: momentum: shell mass-kin. inv. mass {} {}",
+                  mkin / 1_GeV / 1_GeV, m / 1_GeV * m / 1_GeV);
       auto sib_id = process::sibyll::ConvertToSibyllRaw(vP.GetPID());
-      C8LOG_DEBUG(fmt::format("Sibyll::Decay: sib mass: {}", get_sibyll_mass2(sib_id)));
-      C8LOG_DEBUG(fmt::format("Sibyll::Decay: MinStep: gamma:  {}", gamma));
-      C8LOG_DEBUG(fmt::format("Sibyll::Decay: MinStep: tau {} s: ", lifetime / 1_s));
+      C8LOG_DEBUG("Sibyll::Decay: sib mass: {}", get_sibyll_mass2(sib_id));
+      C8LOG_DEBUG("Sibyll::Decay: MinStep: gamma:  {}", gamma);
+      C8LOG_DEBUG("Sibyll::Decay: MinStep: tau {} s: ", lifetime / 1_s);
 
       return lifetime;
     } else
@@ -186,7 +185,7 @@ namespace corsika::process::sibyll {
     if (!IsDecayHandled(pCode))
       throw std::runtime_error("STOP! Sibyll not configured to execute this decay!");
 
-    fCount++;
+    count_++;
     SibStack ss;
     ss.Clear();
 
@@ -206,12 +205,14 @@ namespace corsika::process::sibyll {
     PrintDecayConfig(pCode);
 
     // call sibyll decay
-    C8LOG_DEBUG(fmt::format("Decay: calling Sibyll decay routine.."));
+    C8LOG_DEBUG("Decay: calling Sibyll decay routine..");
     decsib_();
 
-    // print output
-    int print_unit = 6;
-    sib_list_(print_unit);
+    if (sibyll_listing_) {
+      // print output
+      int print_unit = 6;
+      sib_list_(print_unit);
+    }
 
     // reset to stable
     SetStable(pCode);
diff --git a/Processes/Sibyll/Decay.h b/Processes/Sibyll/Decay.h
index 42867ed2b..e5862344e 100644
--- a/Processes/Sibyll/Decay.h
+++ b/Processes/Sibyll/Decay.h
@@ -20,11 +20,12 @@ namespace corsika::process {
   namespace sibyll {
 
     class Decay : public corsika::process::DecayProcess<Decay> {
-      int fCount = 0;
+      int count_ = 0;
       bool handleAllDecays_ = true;
+      bool sibyll_listing_ = false;
 
     public:
-      Decay();
+      Decay(const bool sibyll_listing = false);
       Decay(std::set<particles::Code>);
       ~Decay();
 
diff --git a/Processes/Sibyll/Interaction.cc b/Processes/Sibyll/Interaction.cc
index 8ba3ab187..f020922f0 100644
--- a/Processes/Sibyll/Interaction.cc
+++ b/Processes/Sibyll/Interaction.cc
@@ -33,7 +33,8 @@ namespace corsika::process::sibyll {
 
   bool Interaction::initialized_ = false;
 
-  Interaction::Interaction() {
+  Interaction::Interaction(const bool sibyll_printout_on)
+      : sibyll_listing_(sibyll_printout_on) {
     using random::RNGManager;
 
     // initialize Sibyll
@@ -302,10 +303,12 @@ namespace corsika::process::sibyll {
       // running sibyll, filling stack
       sibyll_(kBeam, targetSibCode, sqs);
 
-      // print final state
-      int print_unit = 6;
-      sib_list_(print_unit);
-      nucCount_ += get_nwounded() - 1;
+      if (sibyll_listing_) {
+        // print final state
+        int print_unit = 6;
+        sib_list_(print_unit);
+        nucCount_ += get_nwounded() - 1;
+      }
 
       // add particles from sibyll to stack
       // link to sibyll stack
diff --git a/Processes/Sibyll/Interaction.h b/Processes/Sibyll/Interaction.h
index 3bcf68807..52a43309d 100644
--- a/Processes/Sibyll/Interaction.h
+++ b/Processes/Sibyll/Interaction.h
@@ -21,9 +21,10 @@ namespace corsika::process::sibyll {
     int count_ = 0;
     int nucCount_ = 0;
     static bool initialized_; ///! flag to assure init is done only once
+    bool sibyll_listing_;
 
   public:
-    Interaction();
+    Interaction(const bool sibyll_printout_on = false);
     ~Interaction();
 
     void SetAllStable();
diff --git a/Processes/TrackingLine/CMakeLists.txt b/Processes/TrackingLine/CMakeLists.txt
index 8411487ff..dd626f877 100644
--- a/Processes/TrackingLine/CMakeLists.txt
+++ b/Processes/TrackingLine/CMakeLists.txt
@@ -21,7 +21,6 @@ set_target_properties (
   PROPERTIES
   VERSION ${PROJECT_VERSION}
   SOVERSION 1
-#  PUBLIC_HEADER "${MODEL_HEADERS}"
   )
 
 # target dependencies on other libraries (also the header onlys)
@@ -33,6 +32,7 @@ target_link_libraries (
   CORSIKAunits
   CORSIKAenvironment
   CORSIKAgeometry
+  CORSIKAlogging
   )
 
 target_include_directories (
diff --git a/Processes/TrackingLine/TrackingLine.cc b/Processes/TrackingLine/TrackingLine.cc
index ee76db62c..23efd2776 100644
--- a/Processes/TrackingLine/TrackingLine.cc
+++ b/Processes/TrackingLine/TrackingLine.cc
@@ -12,6 +12,7 @@
 #include <corsika/geometry/Sphere.h>
 #include <corsika/geometry/Vector.h>
 #include <corsika/process/tracking_line/TrackingLine.h>
+#include <corsika/logging/Logging.h>
 
 #include <limits>
 #include <stdexcept>
diff --git a/Processes/TrackingLine/TrackingLine.h b/Processes/TrackingLine/TrackingLine.h
index 48bdc3a45..80f2b7032 100644
--- a/Processes/TrackingLine/TrackingLine.h
+++ b/Processes/TrackingLine/TrackingLine.h
@@ -14,6 +14,8 @@
 #include <corsika/geometry/Trajectory.h>
 #include <corsika/geometry/Vector.h>
 #include <corsika/units/PhysicalUnits.h>
+#include <corsika/logging/Logging.h>
+
 #include <optional>
 #include <type_traits>
 #include <utility>
@@ -38,7 +40,7 @@ namespace corsika::process {
     class TrackingLine {
 
     public:
-      TrackingLine(){};
+      TrackingLine() = default;
 
       template <typename Particle> // was Stack previously, and argument was
                                    // Stack::StackIterator
@@ -49,27 +51,23 @@ namespace corsika::process {
             p.GetMomentum() / p.GetEnergy() * corsika::units::constants::c;
 
         auto const currentPosition = p.GetPosition();
-        std::cout << "TrackingLine pid: " << p.GetPID()
-                  << " , E = " << p.GetEnergy() / 1_GeV << " GeV" << std::endl;
-        std::cout << "TrackingLine pos: "
-                  << currentPosition.GetCoordinates()
-                  // << " [" << p.GetNode()->GetModelProperties().GetName() << "]"
-                  << std::endl;
-        std::cout << "TrackingLine   E: " << p.GetEnergy() / 1_GeV << " GeV" << std::endl;
-        std::cout << "TrackingLine   p: " << p.GetMomentum().GetComponents() / 1_GeV
-                  << " GeV " << std::endl;
-        std::cout << "TrackingLine   v: " << velocity.GetComponents() << std::endl;
+        C8LOG_DEBUG(
+            "TrackingLine pid: {}"
+            " , E = {} GeV",
+            p.GetPID(), p.GetEnergy() / 1_GeV);
+        C8LOG_DEBUG("TrackingLine pos: {}", currentPosition.GetCoordinates());
+        C8LOG_DEBUG("TrackingLine   E: {} GeV", p.GetEnergy() / 1_GeV);
+        C8LOG_DEBUG("TrackingLine   p: {} GeV", p.GetMomentum().GetComponents() / 1_GeV);
+        C8LOG_DEBUG("TrackingLine   v: {} ", velocity.GetComponents());
 
         // to do: include effect of magnetic field
         geometry::Line line(currentPosition, velocity);
 
         auto const* currentLogicalVolumeNode = p.GetNode();
-        //~ auto const* currentNumericalVolumeNode =
-        //~ fEnvironment.GetUniverse()->GetContainingNode(currentPosition);
         auto const numericallyInside =
             currentLogicalVolumeNode->GetVolume().Contains(currentPosition);
 
-        std::cout << "numericallyInside = " << (numericallyInside ? "true" : "false");
+        C8LOG_DEBUG("numericallyInside = {} ", numericallyInside);
 
         auto const& children = currentLogicalVolumeNode->GetChildNodes();
         auto const& excluded = currentLogicalVolumeNode->GetExcludedNodes();
@@ -85,14 +83,11 @@ namespace corsika::process {
 
           if (auto opt = TimeOfIntersection(line, sphere); opt.has_value()) {
             auto const [t1, t2] = *opt;
-            std::cout << "intersection times: " << t1 / 1_s << "; "
-                      << t2 / 1_s
-                      // << " " << vtn.GetModelProperties().GetName()
-                      << std::endl;
+            C8LOG_DEBUG("intersection times: {} s; {} s", t1 / 1_s, t2 / 1_s);
             if (t1.magnitude() > 0)
               intersections.emplace_back(t1, &vtn);
             else if (t2.magnitude() > 0)
-              std::cout << "inside other volume" << std::endl;
+              C8LOG_DEBUG("inside other volume");
           }
         };
 
@@ -123,10 +118,7 @@ namespace corsika::process {
           min = minIter->first;
         }
 
-        std::cout << " t-intersect: "
-                  << min
-                  // << " " << minIter->second->GetModelProperties().GetName()
-                  << std::endl;
+        C8LOG_DEBUG(" t-intersect: {} ", min);
 
         return std::make_tuple(geometry::Trajectory<geometry::Line>(line, min),
                                velocity.norm() * min, minIter->second);
-- 
GitLab