diff --git a/corsika/detail/framework/geometry/FourVector.inl b/corsika/detail/framework/geometry/FourVector.inl
index 135790bfcfa83efaf1a78f218228d0c2557f5853..6c0a039e5466ef46d2d3014b7e91a6f0621c922c 100644
--- a/corsika/detail/framework/geometry/FourVector.inl
+++ b/corsika/detail/framework/geometry/FourVector.inl
@@ -95,7 +95,7 @@ namespace corsika {
 
   template <typename TTimeType, typename TSpaceVecType>
   typename FourVector<TTimeType, TSpaceVecType>::norm_type
-  FourVector<TTimeType, TSpaceVecType>::operator*(FourVector const& b) {
+      FourVector<TTimeType, TSpaceVecType>::operator*(FourVector const& b) {
     if constexpr (std::is_same<time_type, decltype(std::declval<space_type>() / meter *
                                                    second)>::value)
       return timeLike_ * b.timeLike_ * constants::cSquared - spaceLike_.norm();
diff --git a/corsika/detail/framework/geometry/QuantityVector.inl b/corsika/detail/framework/geometry/QuantityVector.inl
index 6ba02e4f8d4bd8e8cec5231d47ffb23d87931430..2ca24ee5a4a9c080276731e9652f7b03202a82ec 100644
--- a/corsika/detail/framework/geometry/QuantityVector.inl
+++ b/corsika/detail/framework/geometry/QuantityVector.inl
@@ -17,8 +17,8 @@
 namespace corsika {
 
   template <typename TDimension>
-  inline typename QuantityVector<TDimension>::quantity_type
-      QuantityVector<TDimension>::operator[](size_t const index) const {
+  inline typename QuantityVector<TDimension>::quantity_type QuantityVector<TDimension>::
+  operator[](size_t const index) const {
     return quantity_type(phys::units::detail::magnitude_tag, eigenVector_[index]);
   }
 
diff --git a/corsika/detail/media/WeightProvider.inl b/corsika/detail/media/WeightProvider.inl
index c075b236350127c9978fb7c08d87ed035292e5a8..af9849aa3ca75b9bdaeefc661d56232f1ab895f4 100644
--- a/corsika/detail/media/WeightProvider.inl
+++ b/corsika/detail/media/WeightProvider.inl
@@ -20,7 +20,7 @@ namespace corsika {
 
   template <class AConstIterator, class BConstIterator>
   typename WeightProviderIterator<AConstIterator, BConstIterator>::value_type
-  WeightProviderIterator<AConstIterator, BConstIterator>::operator*() const {
+      WeightProviderIterator<AConstIterator, BConstIterator>::operator*() const {
     return ((*aIter_) * (*bIter_)).magnitude();
   }
 
diff --git a/corsika/detail/modules/HadronicElasticModel.inl b/corsika/detail/modules/HadronicElasticModel.inl
index d0a54a711b7b1c62a48dfffbc8bcab60a8a3c5a0..bb5c16133ee16ebdde6cac89de9ad1dd84014191 100644
--- a/corsika/detail/modules/HadronicElasticModel.inl
+++ b/corsika/detail/modules/HadronicElasticModel.inl
@@ -52,7 +52,7 @@ namespace corsika {
           avgCrossSection += getCrossSection(s) * fractions[i];
         }
 
-        std::cout << "avgCrossSection: " << avgCrossSection / 1_mb << " mb" << std::endl;
+        CORSIKA_LOG_DEBUG("avgCrossSection: {} mb", avgCrossSection / 1_mb);
 
         return avgCrossSection;
       }();
@@ -116,7 +116,7 @@ namespace corsika {
     auto const s = static_pow<2>(sqrtS);
 
     auto const B = this->B(s);
-    std::cout << B << std::endl;
+    CORSIKA_LOG_DEBUG(B);
 
     ExponentialDistribution tDist(1 / B);
     auto const absT = [&]() {
@@ -133,10 +133,12 @@ namespace corsika {
       return absT;
     }();
 
-    std::cout << "HadronicElasticInteraction: s = " << s * constants::invGeVsq
-              << " GeV²; absT = " << absT * constants::invGeVsq << " GeV² (max./GeV² = "
-              << 4 * constants::invGeVsq * projectileMomentumSquaredNorm << ')'
-              << std::endl;
+    CORSIKA_LOG_DEBUG(
+        "HadronicElasticInteraction: s = {}"
+        " GeV²; absT = {} "
+        " GeV² (max./GeV² = {})",
+        s * constants::invGeVsq, absT * constants::invGeVsq,
+        4 * constants::invGeVsq * projectileMomentumSquaredNorm);
 
     auto const theta = 2 * asin(sqrt(absT / (4 * pProjectileCoMSqNorm)));
     auto const phi = phiDist(RNG_);
@@ -163,8 +165,8 @@ namespace corsika {
     auto const result =
         (2 * b_p + 2 * b_p + 4 * pow(s * constants::invGeVsq, gfEpsilon) - 4.2) *
         constants::invGeVsq;
-    std::cout << "B(" << s << ") = " << result / constants::invGeVsq << " GeV¯²"
-              << std::endl;
+    CORSIKA_LOG_DEBUG("B({}) = {}  GeV¯²", s, result / constants::invGeVsq);
+
     return result;
   }
 
@@ -180,8 +182,8 @@ namespace corsika {
         static_pow<2>(sigmaTotal) /
         (16 * constants::pi * convert_HEP_to_SI<CrossSectionType::dimension_type>(B(s)));
 
-    std::cout << "HEM sigmaTot = " << sigmaTotal / 1_mb << " mb" << std::endl;
-    std::cout << "HEM sigmaElastic = " << sigmaElastic / 1_mb << " mb" << std::endl;
+    CORSIKA_LOG_DEBUG("HEM sigmaTot = {} mb", sigmaTotal / 1_mb);
+    CORSIKA_LOG_DEBUG("HEM sigmaElastic = {} mb", sigmaElastic / 1_mb);
     return sigmaElastic;
   }
 
diff --git a/corsika/detail/modules/ObservationPlane.inl b/corsika/detail/modules/ObservationPlane.inl
index 74c29b336dc26a3837b036db5bf3af895feb96f4..13f605efca21ea9470b6d3e9d6403e00e9a03298 100644
--- a/corsika/detail/modules/ObservationPlane.inl
+++ b/corsika/detail/modules/ObservationPlane.inl
@@ -60,30 +60,23 @@ namespace corsika {
       corsika::setup::Trajectory const& trajectory) {
 
     auto const& volumeNode = particle.getNode();
-
     typedef typename std::remove_const_t<
         std::remove_reference_t<decltype(volumeNode->getModelProperties())>>
         medium_type;
-
     Intersections const intersection =
         setup::Tracking::intersect<corsika::setup::Stack::particle_type, medium_type>(
             particle, plane_, volumeNode->getModelProperties());
-
     TimeType const timeOfIntersection = intersection.getEntry();
-
     CORSIKA_LOG_TRACE("particle={}, pos={}, dir={}, plane={}, timeOfIntersection={}",
                       particle.asString(), particle.getPosition(),
                       particle.getDirection(), plane_.asString(), timeOfIntersection);
-
     if (timeOfIntersection < TimeType::zero()) {
       return std::numeric_limits<double>::infinity() * 1_m;
     }
     if (timeOfIntersection > trajectory.getDuration()) {
       return std::numeric_limits<double>::infinity() * 1_m;
     }
-
     double const fractionOfIntersection = timeOfIntersection / trajectory.getDuration();
-
     auto const pointOfIntersection = trajectory.getPosition(fractionOfIntersection);
     auto dist = (trajectory.getPosition(0) - pointOfIntersection).getNorm();
     CORSIKA_LOG_TRACE("ObservationPlane: getMaxStepLength l={} m", dist / 1_m);
diff --git a/corsika/detail/modules/OnShellCheck.inl b/corsika/detail/modules/OnShellCheck.inl
index aa6b0453fb24400246c96de206fa415dd0ade7cb..f8f0b0081657c4b7d03d0c66dc2dc71c75a6b71e 100644
--- a/corsika/detail/modules/OnShellCheck.inl
+++ b/corsika/detail/modules/OnShellCheck.inl
@@ -14,8 +14,8 @@
 
 namespace corsika {
 
-  OnShellCheck::OnShellCheck(const double vMassTolerance, const double vEnergyTolerance,
-                             const bool vError)
+  OnShellCheck::OnShellCheck(double const vMassTolerance, double const vEnergyTolerance,
+                             bool const vError)
       : mass_tolerance_(vMassTolerance)
       , energy_tolerance_(vEnergyTolerance)
       , throw_error_(vError) {
@@ -70,16 +70,18 @@ namespace corsika {
         if (abs(e_shift_relative) > energy_tolerance_) {
           logger_->warn("warning! shifted particle energy by {} %",
                         e_shift_relative * 100);
-          if (throw_error_)
+          if (throw_error_) {
             throw std::runtime_error(
                 "OnShellCheck: error! shifted energy by large amount!");
+          }
         }
 
         // reset energy
         p.setEnergy(e_shifted);
-      } else
+      } else {
         CORSIKA_LOGGER_DEBUG(logger_, "particle mass for {} OK", pid);
+      }
     }
-  }
+  } // namespace corsika
 
 } // namespace corsika
diff --git a/corsika/detail/modules/StackInspector.inl b/corsika/detail/modules/StackInspector.inl
index 517e017fceb9baf68c071399c5e67bed8a75b3b6..b657e15316efc371091071ceeaaf265ca9f9a1aa 100644
--- a/corsika/detail/modules/StackInspector.inl
+++ b/corsika/detail/modules/StackInspector.inl
@@ -46,13 +46,16 @@ namespace corsika {
       if (ReportStack_) {
         CoordinateSystemPtr const& rootCS = get_root_CoordinateSystem(); // 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();
+        CORSIKA_LOG_INFO(
+            "StackInspector: i= {:5d}"
+            ", id= {:^15}"
+            " E={:15e} GeV, "
+            " pos= {}"
+            " node = {}",
+            (i++), iterP.getPID(), (E / 1_GeV), pos, fmt::ptr(iterP.getNode()));
+
         if (iterP.getPID() == Code::Nucleus)
-          std::cout << " nuc_ref=" << iterP.getNucleusRef();
-        std::cout << std::endl;
+          CORSIKA_LOG_INFO("nuc_ref= {}", iterP.getNucleusRef());
       }
     }
 
@@ -67,13 +70,18 @@ namespace corsika {
     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;
+    CORSIKA_LOG_INFO(
+        "StackInspector: "
+        " time= {}"
+        ", running= {} seconds"
+        " ( {}%)"
+        ", nStep= {}"
+        ", stackSize= {}"
+        ", Estack= {} GeV"
+        ", ETA=",
+        std::put_time(std::localtime(&now_time), "%T"), elapsed_seconds.count(),
+        int(progress * 100), getStep(), vS.getSize(), Etot / 1_GeV,
+        std::put_time(std::localtime(&eta_time), "%T"));
   }
 
-} // namespace corsika
\ No newline at end of file
+} // namespace corsika
diff --git a/corsika/detail/modules/conex/CONEXhybrid.inl b/corsika/detail/modules/conex/CONEXhybrid.inl
index f722440d304e20a75c83a73e975fe06c1b15b36e..29c3886aab16fccc25a42ed48c8ea2a2ee7cfa15 100644
--- a/corsika/detail/modules/conex/CONEXhybrid.inl
+++ b/corsika/detail/modules/conex/CONEXhybrid.inl
@@ -11,6 +11,7 @@
 #include <corsika/modules/conex/CONEX_f.hpp>
 #include <corsika/framework/random/RNGManager.hpp>
 #include <corsika/framework/core/PhysicalConstants.hpp>
+#include <corsika/framework/core/Logging.hpp>
 
 #include <conexConfig.h>
 
@@ -35,19 +36,11 @@ namespace corsika {
             make_translation(c8cs, translation.getComponents(c8cs));
         auto const transformCS = make_rotationToZ(intermediateCS, translation);
 
-        std::cout << "translation C8/CONEX obs: " << translation.getComponents()
-                  << std::endl;
+        CORSIKA_LOG_DEBUG("translation C8/CONEX obs: ", translation.getComponents());
 
         /*
         auto const transform = CoordinateSystem::getTransformation(
             intermediateCS2, c8cs); // either this way or vice versa... TODO: test this!
-        std::cout << transform.matrix() << std::endl << std::endl;
-        std::cout << CoordinateSystem::getTransformation(intermediateCS, c8cs).matrix()
-                  << std::endl
-                  << std::endl;
-        std::cout << CoordinateSystem::getTransformation(intermediateCS2, intermediateCS)
-                         .matrix()
-                  << std::endl;
         */
         return transformCS;
       })}
@@ -63,27 +56,25 @@ namespace corsika {
       })}
       , y_sf_{showerAxis_.getDirection().cross(x_sf_)} {
 
-    std::cout << "x_sf (conexObservationCS): " << x_sf_.getComponents(conexObservationCS_)
-              << std::endl;
-    std::cout << "x_sf (C8): " << x_sf_.getComponents(center.getCoordinateSystem())
-              << std::endl;
+    CORSIKA_LOG_DEBUG("x_sf (conexObservationCS): {}",
+                      x_sf_.getComponents(conexObservationCS_));
+    CORSIKA_LOG_DEBUG("x_sf (C8): {}", x_sf_.getComponents(center.getCoordinateSystem()));
+
+    CORSIKA_LOG_DEBUG("y_sf (conexObservationCS): {}",
+                      y_sf_.getComponents(conexObservationCS_));
 
-    std::cout << "y_sf (conexObservationCS): " << y_sf_.getComponents(conexObservationCS_)
-              << std::endl;
-    std::cout << "y_sf (C8): " << y_sf_.getComponents(center.getCoordinateSystem())
-              << std::endl;
+    CORSIKA_LOG_DEBUG("y_sf (C8): {}", y_sf_.getComponents(center.getCoordinateSystem()));
 
-    std::cout << "showerAxisDirection (conexObservationCS): "
-              << showerAxis_.getDirection().getComponents(conexObservationCS_)
-              << std::endl;
-    std::cout << "showerAxisDirection (C8): "
-              << showerAxis_.getDirection().getComponents(center.getCoordinateSystem())
-              << std::endl;
+    CORSIKA_LOG_DEBUG("showerAxisDirection (conexObservationCS): {}",
+                      showerAxis_.getDirection().getComponents(conexObservationCS_));
+    CORSIKA_LOG_DEBUG(
+        "showerAxisDirection (C8): {}",
+        showerAxis_.getDirection().getComponents(center.getCoordinateSystem()));
 
-    std::cout << "showerCore (conexObservationCS): "
-              << showerCore_.getCoordinates(conexObservationCS_) << std::endl;
-    std::cout << "showerCore (C8): "
-              << showerCore_.getCoordinates(center.getCoordinateSystem()) << std::endl;
+    CORSIKA_LOG_DEBUG("showerCore (conexObservationCS): {}",
+                      showerCore_.getCoordinates(conexObservationCS_));
+    CORSIKA_LOG_DEBUG("showerCore (C8): {}",
+                      showerCore_.getCoordinates(center.getCoordinateSystem()));
 
     int randomSeeds[3] = {1234, 0, 0}; // will be overwritten later??
     int heModel = eSibyll23;
@@ -114,7 +105,10 @@ namespace corsika {
                             showerAxisConex.getX().magnitude()) *
                  180 / M_PI;
 
-    std::cout << "theta (deg) = " << theta << "; phi (deg) = " << phi << std::endl;
+    CORSIKA_LOG_DEBUG(
+        "theta (deg) = {}"
+        "; phi (deg) = {}",
+        theta, phi);
 
     int ipart = static_cast<int>(primaryPDG);
     auto rng = RNGManager::getInstance().getRandomStream("cascade");
@@ -152,8 +146,8 @@ namespace corsika {
 
     // EM particle
     auto const egs_pid = it->second;
-    std::cout << "position conexObs: " << position.getCoordinates(conexObservationCS_)
-              << std::endl;
+    CORSIKA_LOG_DEBUG("position conexObs: {}",
+                      position.getCoordinates(conexObservationCS_));
 
     auto const coords = position.getCoordinates(conexObservationCS_) / 1_m;
     double const x = coords[0].magnitude();
@@ -187,24 +181,23 @@ namespace corsika {
     double const E = energy / 1_GeV;
     double const m = mass / 1_GeV;
 
-    std::cout << "CONEXhybrid: removing " << egs_pid << " " << std::scientific << energy
-              << " GeV" << std::endl;
-
-    std::cout << "#### parameters to cegs4_() ####" << std::endl;
-    std::cout << "egs_pid = " << egs_pid << std::endl;
-    std::cout << "E = " << E << std::endl;
-    std::cout << "m = " << m << std::endl;
-    std::cout << "x = " << x << std::endl;
-    std::cout << "y = " << y << std::endl;
-    std::cout << "altitude = " << altitude << std::endl;
-    std::cout << "slantDistance = " << slantDistance << std::endl;
-    std::cout << "lateralX = " << lateralX << std::endl;
-    std::cout << "lateralY = " << lateralY << std::endl;
-    std::cout << "slantX = " << slantX << std::endl;
-    std::cout << "time = " << time << std::endl;
-    std::cout << "u = " << u << std::endl;
-    std::cout << "v = " << v << std::endl;
-    std::cout << "w = " << w << std::endl;
+    CORSIKA_LOG_DEBUG("CONEXhybrid: removing {} {:5e} GeV", egs_pid, energy);
+
+    CORSIKA_LOG_DEBUG("#### parameters to cegs4_() ####");
+    CORSIKA_LOG_DEBUG("egs_pid = {}", egs_pid);
+    CORSIKA_LOG_DEBUG("E = {}", E);
+    CORSIKA_LOG_DEBUG("m = {}", m);
+    CORSIKA_LOG_DEBUG("x = {}", x);
+    CORSIKA_LOG_DEBUG("y = {}", y);
+    CORSIKA_LOG_DEBUG("altitude = {}", altitude);
+    CORSIKA_LOG_DEBUG("slantDistance = {}", slantDistance);
+    CORSIKA_LOG_DEBUG("lateralX = {}", lateralX);
+    CORSIKA_LOG_DEBUG("lateralY = {}", lateralY);
+    CORSIKA_LOG_DEBUG("slantX = {}", slantX);
+    CORSIKA_LOG_DEBUG("time = {}", time);
+    CORSIKA_LOG_DEBUG("u = {}", u);
+    CORSIKA_LOG_DEBUG("v = {}", v);
+    CORSIKA_LOG_DEBUG("w = {}", w);
 
     ::conex::cxoptl_.dptl[10 - 1] = egs_pid;
     ::conex::cxoptl_.dptl[4 - 1] = E;
diff --git a/corsika/detail/modules/proposal/ContinuousProcess.inl b/corsika/detail/modules/proposal/ContinuousProcess.inl
index 26f604c998efa2360d4da94adab53e9fac4d28b9..b494c7b72b62a64bf2fdf3cea24f038587c58af1 100644
--- a/corsika/detail/modules/proposal/ContinuousProcess.inl
+++ b/corsika/detail/modules/proposal/ContinuousProcess.inl
@@ -149,9 +149,11 @@ namespace corsika::proposal {
   }
 
   void ContinuousProcess::showResults() const {
-    std::cout << " ******************************" << std::endl
-              << " PROCESS::ContinuousProcess: " << std::endl;
-    std::cout << " energy lost dE (GeV)      :  " << energy_lost_ / 1_GeV << std::endl;
+    CORSIKA_LOG_DEBUG(
+        " ******************************\n"
+        " PROCESS::ContinuousProcess: \n"
+        " energy lost dE (GeV)      :  {}",
+        energy_lost_ / 1_GeV);
   }
 
   void ContinuousProcess::reset() { energy_lost_ = 0_GeV; }
diff --git a/corsika/detail/modules/pythia8/Interaction.inl b/corsika/detail/modules/pythia8/Interaction.inl
index 98ec35a8e31ac728c89f096a39e2e5e7df2eb264..5395e8b1fa86308a4936f6099d34ee851ee67b67 100644
--- a/corsika/detail/modules/pythia8/Interaction.inl
+++ b/corsika/detail/modules/pythia8/Interaction.inl
@@ -21,16 +21,13 @@
 
 namespace corsika::pythia8 {
 
-  Interaction::~Interaction() {
-    std::cout << "Pythia::Interaction n=" << count_ << std::endl;
-  }
+  Interaction::~Interaction() { CORSIKA_LOG_INFO("Pythia::Interaction n= {}", count_); }
 
   Interaction::Interaction(bool const print_listing)
       : Pythia8::Pythia(CORSIKA_Pythia8_XML_DIR)
       , print_listing_(print_listing) {
 
-    std::cout << "Configuring Pythia8 from: " << CORSIKA_Pythia8_XML_DIR << std::endl;
-    std::cout << "Pythia::Interaction n=" << count_ << std::endl;
+    CORSIKA_LOG_INFO("Configuring Pythia8 from: {}", CORSIKA_Pythia8_XML_DIR);
 
     // initialize Pythia
 
@@ -72,12 +69,12 @@ namespace corsika::pythia8 {
   }
 
   void Interaction::setUnstable(Code const pCode) {
-    std::cout << "Pythia::Interaction: setting " << pCode << " unstable.." << std::endl;
+    CORSIKA_LOG_DEBUG("Pythia::Interaction: setting {} unstable..", pCode);
     Pythia8::Pythia::particleData.mayDecay(static_cast<int>(get_PDG(pCode)), true);
   }
 
   void Interaction::setStable(Code const pCode) {
-    std::cout << "Pythia::Interaction: setting " << pCode << " stable.." << std::endl;
+    CORSIKA_LOG_DEBUG("Pythia::Interaction: setting {} stable..", pCode);
     Pythia8::Pythia::particleData.mayDecay(static_cast<int>(get_PDG(pCode)), false);
   }
 
@@ -175,10 +172,12 @@ namespace corsika::pythia8 {
     HEPEnergyType const ECoM = sqrt(
         (Elab + pTotLabNorm) * (Elab - pTotLabNorm)); // binomial for numerical accuracy
 
-    std::cout << "Interaction: LambdaInt: \n"
-              << " input energy: " << particle.getEnergy() / 1_GeV << std::endl
-              << " beam can interact:" << kInteraction << std::endl
-              << " beam pid:" << particle.getPID() << std::endl;
+    CORSIKA_LOG_DEBUG(
+        "Interaction: LambdaInt: \n"
+        " input energy: {} GeV"
+        " beam can interact: {}"
+        " beam pid: {}",
+        particle.getEnergy() / 1_GeV, kInteraction, particle.getPID());
 
     // TODO: move limits into variables
     if (kInteraction && Elab >= 8.5_GeV && isValidCoMEnergy(ECoM)) {
@@ -199,17 +198,16 @@ namespace corsika::pythia8 {
             return std::get<0>(this->getCrossSection(corsikaBeamId, vTargetID, ECoM));
           });
 
-      std::cout << "Interaction: IntLength: weighted CrossSection (mb): "
-                << weightedProdCrossSection / 1_mb << std::endl
-                << "Interaction: IntLength: average mass number: "
-                << mediumComposition.getAverageMassNumber() << std::endl;
+      CORSIKA_LOG_DEBUG(
+          "Interaction: IntLength: weighted CrossSection (mb): {} "
+          "Interaction: IntLength: average mass number: {} ",
+          weightedProdCrossSection / 1_mb, mediumComposition.getAverageMassNumber());
 
       // calculate interaction length in medium
       GrammageType const int_length = mediumComposition.getAverageMassNumber() *
                                       constants::u / weightedProdCrossSection;
-      std::cout << "Interaction: "
-                << "interaction length (g/cm2): " << int_length / (0.001_kg) * 1_cm * 1_cm
-                << std::endl;
+      CORSIKA_LOG_DEBUG("Interaction: interaction length (g/cm2): {} ",
+                        int_length / (0.001_kg) * 1_cm * 1_cm);
 
       return int_length;
     }
@@ -222,10 +220,11 @@ namespace corsika::pythia8 {
 
     auto projectile = view.getProjectile();
 
-    auto const corsikaBeamId = projectile.getPID();
-    std::cout << "Pythia::Interaction: "
-              << "DoInteraction: " << corsikaBeamId << " interaction? "
-              << corsika::pythia8::Interaction::canInteract(corsikaBeamId) << std::endl;
+    const auto corsikaBeamId = projectile.getPID();
+    CORSIKA_LOG_DEBUG(
+        "Pythia::Interaction: "
+        "DoInteraction: {} interaction? ",
+        corsikaBeamId, corsika::pythia8::Interaction::canInteract(corsikaBeamId));
 
     if (is_nucleus(corsikaBeamId)) {
       // nuclei handled by different process, this should not happen
@@ -249,12 +248,15 @@ namespace corsika::pythia8 {
       auto const pTargetLab = MomentumVector(labCS, 0_GeV, 0_GeV, 0_GeV);
       FourVector const PtargLab(eTargetLab, pTargetLab);
 
-      std::cout << "Interaction: ebeam lab: " << eProjectileLab / 1_GeV << std::endl
-                << "Interaction: pbeam lab: " << pProjectileLab.getComponents() / 1_GeV
-                << std::endl;
-      std::cout << "Interaction: etarget lab: " << eTargetLab / 1_GeV << std::endl
-                << "Interaction: ptarget lab: " << pTargetLab.getComponents() / 1_GeV
-                << std::endl;
+      CORSIKA_LOG_DEBUG(
+          "Interaction: ebeam lab: {} GeV"
+          "Interaction: pbeam lab: {} GeV",
+          eProjectileLab / 1_GeV, pProjectileLab.getComponents() / 1_GeV);
+
+      CORSIKA_LOG_DEBUG(
+          "Interaction: etarget lab: {} GeV"
+          "Interaction: ptarget lab: {} GeV ",
+          eTargetLab / 1_GeV, pTargetLab.getComponents() / 1_GeV);
 
       FourVector const PprojLab(eProjectileLab, pProjectileLab);
 
@@ -270,18 +272,20 @@ namespace corsika::pythia8 {
       // boost target
       auto const PtargCoM = boost.toCoM(PtargLab);
 
-      std::cout << "Interaction: ebeam CoM: " << PprojCoM.getTimeLikeComponent() / 1_GeV
-                << std::endl
-                << "Interaction: pbeam CoM: "
-                << PprojCoM.getSpaceLikeComponents().getComponents() / 1_GeV << std::endl;
-      std::cout << "Interaction: etarget CoM: " << PtargCoM.getTimeLikeComponent() / 1_GeV
-                << std::endl
-                << "Interaction: ptarget CoM: "
-                << PtargCoM.getSpaceLikeComponents().getComponents() / 1_GeV << std::endl;
+      CORSIKA_LOG_DEBUG(
+          "Interaction: ebeam CoM: {} GeV"
+          "Interaction: pbeam CoM: {} GeV",
+          PprojCoM.getTimeLikeComponent() / 1_GeV,
+          PprojCoM.getSpaceLikeComponents().getComponents() / 1_GeV);
+
+      CORSIKA_LOG_DEBUG(
+          "Interaction: etarget CoM: {} GeV"
+          "Interaction: ptarget CoM: {} GeV",
+          PtargCoM.getTimeLikeComponent() / 1_GeV,
+          PtargCoM.getSpaceLikeComponents().getComponents() / 1_GeV);
 
-      std::cout << "Interaction: position of interaction: " << pOrig.getCoordinates()
-                << std::endl;
-      std::cout << "Interaction: time: " << tOrig << std::endl;
+      CORSIKA_LOG_DEBUG("Interaction: position of interaction: ", pOrig.getCoordinates());
+      CORSIKA_LOG_DEBUG("Interaction: time: {}", tOrig);
 
       HEPEnergyType Etot = eProjectileLab + eTargetLab;
       MomentumVector Ptot = projectile.getMomentum();
@@ -310,20 +314,23 @@ namespace corsika::pythia8 {
 
       auto const corsikaTargetId =
           mediumComposition.sampleTarget(cross_section_of_components, RNG_);
-      std::cout << "Interaction: target selected: " << corsikaTargetId << std::endl;
+      CORSIKA_LOG_DEBUG("Interaction: target selected: {}", corsikaTargetId);
 
       if (corsikaTargetId != Code::Hydrogen && corsikaTargetId != Code::Neutron &&
           corsikaTargetId != Code::Proton)
         throw std::runtime_error("DoInteraction: wrong target for PYTHIA");
 
-      std::cout << "Interaction: "
-                << " DoInteraction: E(GeV):" << eProjectileLab / 1_GeV
-                << " Ecm(GeV): " << Ecm / 1_GeV << std::endl;
+      CORSIKA_LOG_DEBUG(
+          "Interaction: "
+          " DoInteraction: E(GeV): {}"
+          " Ecm(GeV): {}",
+          eProjectileLab / 1_GeV, Ecm / 1_GeV);
 
       if (eProjectileLab < 8.5_GeV || !isValidCoMEnergy(Ecm)) {
-        std::cout << "Interaction: "
-                  << " DoInteraction: should have dropped particle.. "
-                  << "THIS IS AN ERROR" << std::endl;
+        CORSIKA_LOG_DEBUG(
+            "Interaction: "
+            " DoInteraction: should have dropped particle.. "
+            "THIS IS AN ERROR");
         throw std::runtime_error("energy too low for PYTHIA");
 
       } else {
@@ -363,9 +370,11 @@ namespace corsika::pythia8 {
           Plab_final += pnew.getMomentum();
           Elab_final += pnew.getEnergy();
         }
-        std::cout << "conservation (all GeV): "
-                  << "Elab_final=" << Elab_final / 1_GeV
-                  << ", Plab_final=" << (Plab_final / 1_GeV).getComponents() << std::endl;
+        CORSIKA_LOG_DEBUG(
+            "conservation (all GeV): "
+            "Elab_final= {}"
+            ", Plab_final= {}",
+            Elab_final / 1_GeV, (Plab_final / 1_GeV).getComponents());
       }
     }
   }
diff --git a/corsika/detail/modules/qgsjetII/Interaction.inl b/corsika/detail/modules/qgsjetII/Interaction.inl
index afdafc99f528137d4bf8baa2b06c48594c1f478e..da703d7ef604c47088c495b65674b2f025d9df6d 100644
--- a/corsika/detail/modules/qgsjetII/Interaction.inl
+++ b/corsika/detail/modules/qgsjetII/Interaction.inl
@@ -32,7 +32,7 @@ namespace corsika::qgsjetII {
     if (dataPath == "") {
       if (std::getenv("CORSIKA_DATA")) {
         data_path_ = std::string(std::getenv("CORSIKA_DATA")) + "/QGSJetII/";
-        std::cout << "Searching for QGSJetII data tables in " << data_path_ << std::endl;
+        CORSIKA_LOG_DEBUG("Searching for QGSJetII data tables in {}", data_path_);
       }
     }
 
@@ -47,7 +47,7 @@ namespace corsika::qgsjetII {
   }
 
   Interaction::~Interaction() {
-    std::cout << "QgsjetII::Interaction n=" << count_ << std::endl;
+    CORSIKA_LOG_DEBUG("QgsjetII::Interaction n= {}", count_);
   }
 
   CrossSectionType Interaction::getCrossSection(const Code beamId, const Code targetId,
@@ -76,10 +76,12 @@ namespace corsika::qgsjetII {
           throw std::runtime_error("QgsjetII target outside range. ");
       }
 
-      std::cout << "QgsjetII::getCrossSection Elab=" << Elab << " iBeam=" << iBeam
-                << " iProjectile=" << iProjectile << " iTarget=" << iTarget << std::endl;
+      CORSIKA_LOG_DEBUG(
+          "QgsjetII::getCrossSection Elab= {} GeV iBeam= {}"
+          " iProjectile= {} iTarget= {}",
+          Elab / 1_GeV, iBeam, iProjectile, iTarget);
       sigProd = qgsect_(Elab / 1_GeV, iBeam, iProjectile, iTarget);
-      std::cout << "QgsjetII::getCrossSection sigProd=" << sigProd << std::endl;
+      CORSIKA_LOG_DEBUG("QgsjetII::getCrossSection sigProd= {} mb", sigProd);
     }
 
     return sigProd * 1_mb;
@@ -103,10 +105,12 @@ namespace corsika::qgsjetII {
     // total momentum and energy
     HEPEnergyType Elab = vP.getEnergy();
 
-    std::cout << "Interaction: LambdaInt: \n"
-              << " input energy: " << vP.getEnergy() / 1_GeV << std::endl
-              << " beam can interact:" << kInteraction << std::endl
-              << " beam pid:" << vP.getPID() << std::endl;
+    CORSIKA_LOG_DEBUG(
+        "Interaction: LambdaInt: \n"
+        " input energy: {} GeV"
+        " beam can interact: {}"
+        " beam pid: {}",
+        vP.getEnergy() / 1_GeV, kInteraction, vP.getPID());
 
     if (kInteraction) {
 
@@ -131,16 +135,18 @@ namespace corsika::qgsjetII {
             return getCrossSection(corsikaBeamId, targetID, Elab, Abeam, targetA);
           });
 
-      std::cout << "Interaction: "
-                << "IntLength: weighted CrossSection (mb): "
-                << weightedProdCrossSection / 1_mb << std::endl;
+      CORSIKA_LOG_DEBUG(
+          "Interaction: "
+          "IntLength: weighted CrossSection (mb): {}",
+          weightedProdCrossSection / 1_mb);
 
       // calculate interaction length in medium
       GrammageType const int_length = mediumComposition.getAverageMassNumber() *
                                       constants::u / weightedProdCrossSection;
-      std::cout << "Interaction: "
-                << "interaction length (g/cm2): " << int_length / (0.001_kg) * 1_cm * 1_cm
-                << std::endl;
+      CORSIKA_LOG_DEBUG(
+          "Interaction: "
+          "interaction length (g/cm2): {}",
+          int_length / (0.001_kg) * 1_cm * 1_cm);
 
       return int_length;
     }
@@ -159,9 +165,10 @@ namespace corsika::qgsjetII {
     auto const projectile = view.getProjectile();
 
     const auto corsikaBeamId = projectile.getPID();
-    std::cout << "ProcessQgsjetII: "
-              << "DoInteraction: " << corsikaBeamId << " interaction? "
-              << corsika::qgsjetII::canInteract(corsikaBeamId) << std::endl;
+    CORSIKA_LOG_DEBUG(
+        "ProcessQgsjetII: "
+        "DoInteraction: {} interaction possible? {}",
+        corsikaBeamId, corsika::qgsjetII::canInteract(corsikaBeamId));
 
     if (corsika::qgsjetII::canInteract(corsikaBeamId)) {
 
@@ -187,16 +194,17 @@ namespace corsika::qgsjetII {
 
       HEPEnergyType const projectileEnergyLabPerNucleon = projectileEnergyLab / beamA;
 
-      std::cout << "Interaction: ebeam lab: " << projectileEnergyLab / 1_GeV << std::endl
-                << "Interaction: pbeam lab: "
-                << projectileMomentumLab.getComponents() / 1_GeV << std::endl;
-      std::cout << "Interaction: etarget lab: " << targetEnergyLab / 1_GeV << std::endl
-                << "Interaction: ptarget lab: "
-                << targetMomentumLab.getComponents() / 1_GeV << std::endl;
-
-      std::cout << "Interaction: position of interaction: " << pOrig.getCoordinates()
-                << std::endl;
-      std::cout << "Interaction: time: " << tOrig << std::endl;
+      CORSIKA_LOG_DEBUG(
+          "Interaction: ebeam lab: {} GeV"
+          "Interaction: pbeam lab: {} GeV",
+          projectileEnergyLab / 1_GeV, projectileMomentumLab.getComponents() / 1_GeV);
+      CORSIKA_LOG_DEBUG(
+          "Interaction: etarget lab: {} GeV"
+          "Interaction: ptarget lab: {} GeV",
+          targetEnergyLab / 1_GeV, targetMomentumLab.getComponents() / 1_GeV);
+      CORSIKA_LOG_DEBUG("Interaction: position of interaction: {}",
+                        pOrig.getCoordinates());
+      CORSIKA_LOG_DEBUG("Interaction: time: {} ", tOrig);
 
       // sample target mass number
       auto const* currentNode = projectile.getNode();
@@ -221,7 +229,7 @@ namespace corsika::qgsjetII {
 
       const auto targetCode =
           mediumComposition.sampleTarget(cross_section_of_components, rng_);
-      std::cout << "Interaction: target selected: " << targetCode << std::endl;
+      CORSIKA_LOG_DEBUG("Interaction: target selected: {}", targetCode);
 
       int targetMassNumber = 1;     // proton
       if (is_nucleus(targetCode)) { // nucleus
@@ -232,8 +240,7 @@ namespace corsika::qgsjetII {
         if (targetCode != Proton::code)
           throw std::runtime_error("QgsjetII Taget not possible.");
       }
-      std::cout << "Interaction: target qgsjetII code/A: " << targetMassNumber
-                << std::endl;
+      CORSIKA_LOG_DEBUG("Interaction: target qgsjetII code/A: ", targetMassNumber);
 
       int projectileMassNumber = 1; // "1" means "hadron"
       QgsjetIIHadronType qgsjet_hadron_type =
@@ -274,8 +281,8 @@ namespace corsika::qgsjetII {
         // else if (abs(kBeam)>6) -> throw
       }
 
-      std::cout << "Interaction: "
-                << " DoInteraction: E(GeV):" << projectileEnergyLab / 1_GeV << std::endl;
+      CORSIKA_LOG_DEBUG("Interaction: DoInteraction: E(GeV): {} GeV",
+                        projectileEnergyLab / 1_GeV);
       count_++;
       int qgsjet_hadron_type_int = static_cast<QgsjetIICodeIntType>(qgsjet_hadron_type);
       qgini_(projectileEnergyLab / 1_GeV, qgsjet_hadron_type_int, projectileMassNumber,
@@ -320,8 +327,10 @@ namespace corsika::qgsjetII {
 
             auto const energy = sqrt(momentum.getSquaredNorm() + square(nucleonMass));
             momentum.rebase(originalCS); // transform back into standard lab frame
-            std::cout << "secondary fragment> id=" << idFragm
-                      << " p=" << momentum.getComponents() << std::endl;
+            CORSIKA_LOG_DEBUG(
+                "secondary fragment> id= {}"
+                " p={}",
+                idFragm, momentum.getComponents());
             auto pnew = view.addSecondary(
                 std::make_tuple(idFragm, energy, momentum, pOrig, tOrig));
             Plab_final += pnew.getMomentum();
@@ -352,9 +361,13 @@ namespace corsika::qgsjetII {
 
           auto const energy = sqrt(momentum.getSquaredNorm() + square(nucleusMass));
           momentum.rebase(originalCS); // transform back into standard lab frame
-          std::cout << "secondary fragment> id=" << idFragm
-                    << " p=" << momentum.getComponents() << " A=" << A << " Z=" << Z
-                    << std::endl;
+          CORSIKA_LOG_DEBUG(
+              "secondary fragment> id={}"
+              " p={}"
+              " A={}"
+              " Z= {}",
+              idFragm, momentum.getComponents(), A, Z);
+
           auto pnew = view.addSecondary(
               std::make_tuple(idFragm, energy, momentum, pOrig, tOrig, A, Z));
           Plab_final += pnew.getMomentum();
@@ -370,24 +383,27 @@ namespace corsika::qgsjetII {
         auto const energy = psec.getEnergy();
 
         momentum.rebase(originalCS); // transform back into standard lab frame
-        std::cout << "secondary fragment> id="
-                  << corsika::qgsjetII::convertFromQgsjetII(psec.getPID())
-                  << " p=" << momentum.getComponents() << std::endl;
+        CORSIKA_LOG_DEBUG(
+            "secondary fragment> id= {}"
+            " p= {}",
+            corsika::qgsjetII::convertFromQgsjetII(psec.getPID()),
+            momentum.getComponents());
         auto pnew = view.addSecondary(
             std::make_tuple(corsika::qgsjetII::convertFromQgsjetII(psec.getPID()), energy,
                             momentum, pOrig, tOrig));
         Plab_final += pnew.getMomentum();
         Elab_final += pnew.getEnergy();
       }
-      std::cout << "conservation (all GeV): Ecm_final= n/a" /* << Ecm_final / 1_GeV*/
-                << std::endl
-                << "Elab_final=" << Elab_final / 1_GeV
-                << ", Plab_final=" << (Plab_final / 1_GeV).getComponents()
-                << ", N_wounded,targ="
-                << QGSJetIIFragmentsStackData::getWoundedNucleonsTarget()
-                << ", N_wounded,proj="
-                << QGSJetIIFragmentsStackData::getWoundedNucleonsProjectile()
-                << ", N_fragm,proj=" << qfs.getSize() << std::endl;
+      CORSIKA_LOG_DEBUG(
+          "conservation (all GeV): Ecm_final= n/a" /* << Ecm_final / 1_GeV*/
+          "Elab_final="
+          ", Plab_final={}"
+          ", N_wounded,targ={}"
+          ", N_wounded,proj={}"
+          ", N_fragm,proj={}",
+          Elab_final / 1_GeV, (Plab_final / 1_GeV).getComponents(),
+          QGSJetIIFragmentsStackData::getWoundedNucleonsTarget(),
+          QGSJetIIFragmentsStackData::getWoundedNucleonsProjectile(), qfs.getSize());
     }
   }
 
diff --git a/corsika/detail/modules/qgsjetII/qgsjet-II-04.inl b/corsika/detail/modules/qgsjetII/qgsjet-II-04.inl
index 130ac0173920a78bdd774914c5f285a90caa51ff..aca1429c3bfd93fac96a6519746eb88714ebbf09 100644
--- a/corsika/detail/modules/qgsjetII/qgsjet-II-04.inl
+++ b/corsika/detail/modules/qgsjetII/qgsjet-II-04.inl
@@ -11,14 +11,14 @@
 #include <corsika/modules/qgsjetII/qgsjet-II-04.hpp>
 
 #include <corsika/framework/random/RNGManager.hpp>
+#include <corsika/framework/core/Logging.hpp>
 
 #include <iostream>
 #include <random>
 
 datadir::datadir(const std::string& dir) {
   if (dir.length() > 130) {
-    std::cerr << "QGSJetII error, will cut datadir \"" << dir
-              << "\" to 130 characters: " << std::endl;
+    CORSIKA_LOG_ERROR("QGSJetII error, will cut datadir \"{}\" to 130 characters: ", {});
   }
   int i = 0;
   for (i = 0; i < std::min(130, int(dir.length())); ++i) data[i] = dir[i];
diff --git a/corsika/detail/modules/sibyll/Interaction.inl b/corsika/detail/modules/sibyll/Interaction.inl
index e2c5f001fcb12dc487558b5efad6a69e4bbdcb37..69c4d43393a8812b11794a7db93ee6f79f6117cb 100644
--- a/corsika/detail/modules/sibyll/Interaction.inl
+++ b/corsika/detail/modules/sibyll/Interaction.inl
@@ -41,8 +41,7 @@ namespace corsika::sibyll {
   }
 
   inline Interaction::~Interaction() {
-    CORSIKA_LOG_DEBUG(
-        fmt::format("Sibyll::Interaction n={}, Nnuc={}", count_, nucCount_));
+    CORSIKA_LOG_DEBUG("Sibyll::Interaction n={}, Nnuc={}", count_, nucCount_);
   }
 
   inline void Interaction::setStable(std::vector<corsika::Code> const& vParticleList) {
@@ -54,13 +53,13 @@ namespace corsika::sibyll {
   }
 
   inline void Interaction::setUnstable(const corsika::Code vCode) {
-    std::cout << "Sibyll::Interaction: setting " << vCode << " unstable.." << std::endl;
+    CORSIKA_LOG_DEBUG("Sibyll::Interaction: setting {} unstable..", vCode);
     const int s_id = abs(corsika::sibyll::convertToSibyllRaw(vCode));
     s_csydec_.idb[s_id - 1] = abs(s_csydec_.idb[s_id - 1]);
   }
 
   inline void Interaction::setStable(const corsika::Code vCode) {
-    std::cout << "Sibyll::Interaction: setting " << vCode << " stable.." << std::endl;
+    CORSIKA_LOG_DEBUG("Sibyll::Interaction: setting {} stable..", vCode);
     const int s_id = abs(corsika::sibyll::convertToSibyllRaw(vCode));
     s_csydec_.idb[s_id - 1] = (-1) * abs(s_csydec_.idb[s_id - 1]);
   }
@@ -166,17 +165,17 @@ namespace corsika::sibyll {
           });
 
       CORSIKA_LOG_DEBUG(
-          fmt::format("Interaction: "
-                      "IntLength: weighted CrossSection (mb): {} ",
-                      weightedProdCrossSection / 1_mb));
+          "Interaction: "
+          "IntLength: weighted CrossSection (mb): {} ",
+          weightedProdCrossSection / 1_mb);
 
       // calculate interaction length in medium
       GrammageType const int_length = mediumComposition.getAverageMassNumber() *
                                       constants::u / weightedProdCrossSection;
       CORSIKA_LOG_DEBUG(
-          fmt::format("Interaction: "
-                      "interaction length (g/cm2): {} ",
-                      int_length / (0.001_kg) * 1_cm * 1_cm));
+          "Interaction: "
+          "interaction length (g/cm2): {} ",
+          int_length / (0.001_kg) * 1_cm * 1_cm);
 
       return int_length;
     }
diff --git a/corsika/detail/modules/sibyll/NuclearInteraction.inl b/corsika/detail/modules/sibyll/NuclearInteraction.inl
index 52296fbd423199e5398e45d0e72ed53ae2e01622..e5d8662fbab3fec61f00846e2dfd1511c414acce 100644
--- a/corsika/detail/modules/sibyll/NuclearInteraction.inl
+++ b/corsika/detail/modules/sibyll/NuclearInteraction.inl
@@ -577,7 +577,7 @@ namespace corsika::sibyll {
     }
 
     // add inelastic interactions
-    std::cout << "calculate inelastic nucleon-nucleon interactions.." << std::endl;
+    CORSIKA_LOG_DEBUG("calculate inelastic nucleon-nucleon interactions..");
     for (int j = 0; j < nInelNucleons; ++j) {
       // TODO: sample neutron or proton
       auto pCode = Code::Proton;
diff --git a/corsika/detail/modules/urqmd/UrQMD.inl b/corsika/detail/modules/urqmd/UrQMD.inl
index 372d2d5925450b4798bd4650015e856e643d5de4..c0f530a811a12809b68ccb5caa125e1d7319ef11 100644
--- a/corsika/detail/modules/urqmd/UrQMD.inl
+++ b/corsika/detail/modules/urqmd/UrQMD.inl
@@ -223,14 +223,12 @@ namespace corsika::urqmd {
       auto const energy = sqrt(momentum.getSquaredNorm() + square(get_mass(code)));
 
       momentum.rebase(originalCS); // transform back into standard lab frame
-      std::cout << i << " " << code << " " << momentum.getComponents() << std::endl;
+      CORSIKA_LOG_DEBUG(" {} {} {} ", i, code, momentum.getComponents());
 
       projectile.addSecondary(
           std::make_tuple(code, energy, momentum, projectilePosition, projectileTime));
     }
-
-    std::cout << "UrQMD generated " << ::urqmd::sys_.npart << " secondaries!"
-              << std::endl;
+    CORSIKA_LOG_DEBUG("UrQMD generated {} secondaries!", ::urqmd::sys_.npart);
   }
 
   Code convertFromUrQMD(int vItyp, int vIso3) {
diff --git a/corsika/framework/utility/detail/COMBoost.inl b/corsika/framework/utility/detail/COMBoost.inl
deleted file mode 100644
index aa6af48349ed06df019e709eaadf2d8693b56289..0000000000000000000000000000000000000000
--- a/corsika/framework/utility/detail/COMBoost.inl
+++ /dev/null
@@ -1,125 +0,0 @@
-/*
- * (c) Copyright 2018 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/framework/geometry/CoordinateSystem.hpp>
-#include <corsika/framework/geometry/FourVector.hpp>
-#include <corsika/framework/geometry/FourVector.hpp>
-#include <corsika/framework/geometry/Vector.hpp>
-#include <corsika/framework/core/PhysicalUnits.hpp>
-#include <corsika/framework/utility/sgn.hpp>
-#include <Eigen/Dense>
-#include "../../core/PhysicalUnits.hpp"
-
-namespace corsika::utl {
-
-  auto const& COMBoost::GetRotationMatrix() const { return fRotation; }
-
-  //! transforms a 4-momentum from lab frame to the center-of-mass frame
-  template <typename FourVector>
-  FourVector COMBoost::toCoM(const FourVector& p) const {
-    using namespace corsika::units::si;
-    auto pComponents = p.GetSpaceLikeComponents().GetComponents(rotatedCS_);
-    Eigen::Vector3d eVecRotated = pComponents.eVector;
-    Eigen::Vector2d lab;
-
-    lab << (p.GetTimeLikeComponent() * (1 / 1_GeV)),
-        (eVecRotated(2) * (1 / 1_GeV).magnitude());
-
-    auto const boostedZ = boost_ * lab;
-    auto const E_CoM = boostedZ(0) * 1_GeV;
-
-    eVecRotated(2) = boostedZ(1) * (1_GeV).magnitude();
-
-    return FourVector(E_CoM,
-                      corsika::geometry::Vector<hepmomentum_d>(rotatedCS_, eVecRotated));
-  }
-
-  //! transforms a 4-momentum from the center-of-mass frame back to lab frame
-  template <typename FourVector>
-  FourVector COMBoost::fromCoM(const FourVector& p) const {
-    using namespace corsika::units::si;
-    auto pCM = p.GetSpaceLikeComponents().GetComponents(rotatedCS_);
-    auto const Ecm = p.GetTimeLikeComponent();
-
-    Eigen::Vector2d com;
-    com << (Ecm * (1 / 1_GeV)), (pCM.eVector(2) * (1 / 1_GeV).magnitude());
-
-    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;
-
-    pCM.eVector(2) = boostedZ(1) * (1_GeV).magnitude();
-
-    geometry::Vector<typename decltype(pCM)::dimension> pLab{rotatedCS_, pCM};
-    pLab.rebase(originalCS_);
-
-    FourVector f(E_lab, pLab);
-
-    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;
-  }
-
-  COMBoost::COMBoost(FourVector<HEPEnergyType, Vector<hepmomentum_d>> const& Pprojectile,
-                     const HEPMassType massTarget)
-      : fCS(Pprojectile.GetSpaceLikeComponents().GetCoordinateSystem()) {
-    auto const pProjectile = Pprojectile.GetSpaceLikeComponents();
-    auto const pProjNorm = pProjectile.norm();
-    auto const a = (pProjectile / pProjNorm).GetComponents().eVector;
-    auto const a1 = a(0), a2 = a(1);
-
-    auto const s = sgn(a(2));
-    auto const c = 1 / (1 + s * a(2));
-
-    Eigen::Matrix3d A, B;
-
-    if (s > 0) {
-      A << 1, 0, -a1,                     // comment to prevent clang-format
-          0, 1, -a2,                      // .
-          a1, a2, 1;                      // .
-      B << -a1 * a1 * c, -a1 * a2 * c, 0, // .
-          -a1 * a2 * c, -a2 * a2 * c, 0,  // .
-          0, 0, -(a1 * a1 + a2 * a2) * c; // .
-
-    } else {
-      A << 1, 0, a1,                      // comment to prevent clang-format
-          0, -1, -a2,                     // .
-          a1, a2, -1;                     // .
-      B << -a1 * a1 * c, -a1 * a2 * c, 0, // .
-          +a1 * a2 * c, +a2 * a2 * c, 0,  // .
-          0, 0, (a1 * a1 + a2 * a2) * c;  // .
-    }
-
-    fRotation = A + B;
-
-    // calculate boost
-    double const beta = pProjNorm / (Pprojectile.GetTimeLikeComponent() + massTarget);
-
-    /* Accurracy matters here, beta = 1 - epsilon for ultra-relativistic boosts */
-    double const coshEta = 1 / std::sqrt((1 + beta) * (1 - beta));
-    //~ double const coshEta = 1 / std::sqrt((1-beta*beta));
-    double const sinhEta = -beta * coshEta;
-
-    std::cout << "COMBoost (1-beta)=" << 1 - beta << " gamma=" << coshEta << std::endl;
-    std::cout << "  det = " << fRotation.determinant() - 1 << std::endl;
-
-    fBoost << coshEta, sinhEta, sinhEta, coshEta;
-
-    fInverseBoost << coshEta, -sinhEta, -sinhEta, coshEta;
-  }
-} // namespace corsika::utl
diff --git a/corsika/framework/utility/detail/CorsikaFenvFallback.inl b/corsika/framework/utility/detail/CorsikaFenvFallback.inl
deleted file mode 100644
index ba847a3a04dd63dc980b95148113f440e8d70aaf..0000000000000000000000000000000000000000
--- a/corsika/framework/utility/detail/CorsikaFenvFallback.inl
+++ /dev/null
@@ -1,18 +0,0 @@
-/*
- * (c) Copyright 2018 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.
- */
-
-#include <corsika/framework/utility/CorsikaFenv.hpp>
-#include <cfenv>
-
-extern "C" {
-#warning No enabling/disabling of floating point exceptions - platform needs better implementation
-
-int feenableexcept(int excepts) { return -1; }
-
-int fedisableexcept(int excepts) { return -1; }
-}
diff --git a/do-clang-format.py b/do-clang-format.py
index 076a5ccdf4c207ff448aef2453fa835233c0ff28..24973f6690e84137227763187cd3101dc56e4ac1 100755
--- a/do-clang-format.py
+++ b/do-clang-format.py
@@ -15,6 +15,8 @@ parser.add_argument('--apply', action="store_true",
     help="Apply clang-format to files which need changes.")
 parser.add_argument("--all", action="store_true",
     help="Check all files below current path instead of new/modified.")
+parser.add_argument("--docker", action="store_true",
+    help="Use corsika/devel:clang-8 container to run clang-format. Make sure you are in group \"docker\". ")
 
 args = parser.parse_args()
 
@@ -70,11 +72,22 @@ else:
 cmd = "clang-format"
 if "CLANG_FORMAT" in os.environ:
   cmd = os.environ["CLANG_FORMAT"]
-cmd +=  " -style=file"
+if args.docker: 
+  USER=os.environ["USER"]
+  UID=os.getuid()
+  GID=os.getgid()
+  PWD=os.getcwd()
+  # note, currently in container it is clang-8 
+  cmd = "docker container run --rm -v {}:/corsika -w /corsika -u {}:{} corsika/devel:clang-8 clang-format-8".format(PWD,UID,GID)
+cmd += " -style=file"
+
+version = subp.check_output(cmd.split() + ["--version"]).decode("utf-8")
+print (version)
+
 if args.apply:
     for filename in filelist:        
         subp.check_call(cmd.split() + ["-i", filename])
-
+        
 else:
     # only print files which need formatting
     files_need_formatting = 0
diff --git a/tests/framework/testFunctionTimer.cpp b/tests/framework/testFunctionTimer.cpp
index f0bc38e67ec68a7aa85c898f007e1fd405b4dc1c..46d97f710f5793a5f27bbc4ed255c49cbdad627b 100644
--- a/tests/framework/testFunctionTimer.cpp
+++ b/tests/framework/testFunctionTimer.cpp
@@ -12,7 +12,6 @@
 #include <catch2/catch.hpp>
 
 #include <chrono>
-#include <iostream>
 #include <thread>
 
 using namespace corsika;
@@ -39,15 +38,15 @@ TEST_CASE("FunctionTimer", "[Timer]") {
 
     auto test = corsika::FunctionTimer(testFunc);
 
-    std::cout << test() << std::endl;
-    std::cout << test.getTime().count() << std::endl;
+    CORSIKA_LOG_DEBUG(test());
+    CORSIKA_LOG_DEBUG(test.getTime().count());
   }
 
   SECTION("Measure runtime of a class functor") {
     TestClass testC;
     auto test = corsika::FunctionTimer(testC);
 
-    std::cout << test() << std::endl;
-    std::cout << test.getTime().count() << std::endl;
+    CORSIKA_LOG_DEBUG(test());
+    CORSIKA_LOG_DEBUG(test.getTime().count());
   }
 }
diff --git a/tests/framework/testInteractionCounter.cpp b/tests/framework/testInteractionCounter.cpp
index bf73fc578316eb3c7e948c22671f6a1aad220ddc..556d4733b60b82bf295ff47463a24d3775c48553 100644
--- a/tests/framework/testInteractionCounter.cpp
+++ b/tests/framework/testInteractionCounter.cpp
@@ -83,7 +83,7 @@ TEST_CASE("InteractionCounter", "[process]") {
       auto const file = GENERATE(as<std::string>{}, "testInteractionCounter_file1",
                                  "testInteractionCounter_file2");
 
-      std::cout << file + ".npz vs " << refDataDir + "/" + file + "_REF.npz" << std::endl;
+      CORSIKA_LOG_INFO("{0}.npz vs {1}/{0}_REF.npz", file, refDataDir);
 
       // compare to binary reference data
       // note that this currenly compares the whole files byte by byte. If the new
diff --git a/tests/framework/testProcessSequence.cpp b/tests/framework/testProcessSequence.cpp
index d2612bb8860497c34aa1550b1eb227ba2d108c68..3d60d2de15a57a54a0e0565ad4b568f2a63a944c 100644
--- a/tests/framework/testProcessSequence.cpp
+++ b/tests/framework/testProcessSequence.cpp
@@ -39,7 +39,10 @@ public:
       : v_(v)
       , step_(step) {
 
-    cout << "globalCount: " << globalCount << ", v_: " << v_ << std::endl;
+    CORSIKA_LOG_DEBUG(
+        "globalCount: {} "
+        ", v_: {} ",
+        globalCount, v_);
     globalCount++;
   }
 
@@ -48,7 +51,7 @@ public:
   template <typename D, typename T>
   inline ProcessReturn doContinuous(D& d, T&, bool const flag) const {
     flag_ = flag;
-    cout << "ContinuousProcess1::DoContinuous" << endl;
+    CORSIKA_LOG_TRACE("ContinuousProcess1::DoContinuous");
     checkCont |= 1;
     for (int i = 0; i < nData; ++i) d.data_[i] += 0.933;
     return ProcessReturn::Ok;
@@ -73,7 +76,10 @@ public:
   ContinuousProcess2(int const v, LengthType const step)
       : v_(v)
       , step_(step) {
-    cout << "globalCount: " << globalCount << ", v_: " << v_ << std::endl;
+    CORSIKA_LOG_DEBUG(
+        "globalCount: {}"
+        ", v_: {}",
+        globalCount, v_);
     globalCount++;
   }
 
@@ -82,7 +88,7 @@ public:
   template <typename D, typename T>
   inline ProcessReturn doContinuous(D& d, T&, bool const flag) const {
     flag_ = flag;
-    cout << "ContinuousProcess2::DoContinuous" << endl;
+    CORSIKA_LOG_DEBUG("ContinuousProcess2::DoContinuous");
     checkCont |= 2;
     for (int i = 0; i < nData; ++i) d.data_[i] += 0.111;
     return ProcessReturn::Ok;
@@ -107,7 +113,10 @@ public:
   ContinuousProcess3(int const v, LengthType const step)
       : v_(v)
       , step_(step) {
-    cout << "globalCount: " << globalCount << ", v_: " << v_ << std::endl;
+    CORSIKA_LOG_DEBUG(
+        "globalCount: {}"
+        ", v_: {} ",
+        globalCount, v_);
     globalCount++;
   }
 
@@ -116,7 +125,7 @@ public:
   template <typename D, typename T>
   inline ProcessReturn doContinuous(D& d, T&, bool const flag) const {
     flag_ = flag;
-    cout << "ContinuousProcess3::DoContinuous" << endl;
+    CORSIKA_LOG_DEBUG("ContinuousProcess3::DoContinuous");
     checkCont |= 4;
     for (int i = 0; i < nData; ++i) d.data_[i] += 0.333;
     return ProcessReturn::Ok;
@@ -140,7 +149,11 @@ class Process1 : public InteractionProcess<Process1> {
 public:
   Process1(int const v)
       : v_(v) {
-    cout << "globalCount: " << globalCount << ", v_: " << v_ << std::endl;
+    CORSIKA_LOG_DEBUG(
+        "globalCount: {}"
+        ", v_: {}",
+        globalCount, v_);
+    ;
     globalCount++;
   }
 
@@ -163,7 +176,10 @@ class Process2 : public InteractionProcess<Process2> {
 public:
   Process2(int const v)
       : v_(v) {
-    cout << "globalCount: " << globalCount << ", v_: " << v_ << std::endl;
+    CORSIKA_LOG_DEBUG(
+        "globalCount: {}"
+        ", v_: {}",
+        globalCount, v_);
     globalCount++;
   }
 
@@ -171,11 +187,11 @@ public:
   inline void doInteraction(TView& v) const {
     checkInteract |= 2;
     for (int i = 0; i < nData; ++i) v.parent().data_[i] /= 1.1;
-    cout << "Process2::doInteraction" << endl;
+    CORSIKA_LOG_DEBUG("Process2::doInteraction");
   }
   template <typename Particle>
   GrammageType getInteractionLength(Particle&) const {
-    cout << "Process2::GetInteractionLength" << endl;
+    CORSIKA_LOG_DEBUG("Process2::GetInteractionLength");
     return 20_g / (1_cm * 1_cm);
   }
 
@@ -187,7 +203,10 @@ class Process3 : public InteractionProcess<Process3> {
 public:
   Process3(int const v)
       : v_(v) {
-    cout << "globalCount: " << globalCount << ", v_: " << v_ << std::endl;
+    CORSIKA_LOG_DEBUG(
+        "globalCount: {}"
+        ", v_: {}",
+        globalCount, v_);
     globalCount++;
   }
 
@@ -195,11 +214,11 @@ public:
   inline void doInteraction(TView& v) const {
     checkInteract |= 4;
     for (int i = 0; i < nData; ++i) v.parent().data_[i] *= 1.01;
-    cout << "Process3::doInteraction" << endl;
+    CORSIKA_LOG_DEBUG("Process3::doInteraction");
   }
   template <typename Particle>
   GrammageType getInteractionLength(Particle&) const {
-    cout << "Process3::GetInteractionLength" << endl;
+    CORSIKA_LOG_DEBUG("Process3::GetInteractionLength");
     return 30_g / (1_cm * 1_cm);
   }
 
@@ -211,13 +230,16 @@ class Process4 : public BaseProcess<Process4> {
 public:
   Process4(int const v)
       : v_(v) {
-    cout << "globalCount: " << globalCount << ", v_: " << v_ << std::endl;
+    CORSIKA_LOG_DEBUG(
+        "globalCount: {}"
+        ", v_: {}",
+        globalCount, v_);
     globalCount++;
   }
 
   template <typename D, typename T>
   inline ProcessReturn doContinuous(D& d, T&, bool const) const {
-    std::cout << "Base::doContinuous" << std::endl;
+    CORSIKA_LOG_DEBUG("Base::doContinuous");
     checkCont |= 8;
     for (int i = 0; i < nData; ++i) { d.data_[i] /= 1.2; }
     return ProcessReturn::Ok;
@@ -234,7 +256,7 @@ private:
 class Decay1 : public DecayProcess<Decay1> {
 public:
   Decay1(int const) {
-    cout << "Decay1()" << endl;
+    CORSIKA_LOG_DEBUG("Decay1()");
     globalCount++;
   }
 
@@ -251,7 +273,7 @@ public:
 class Decay2 : public DecayProcess<Decay2> {
 public:
   Decay2(int const) {
-    cout << "Decay2()" << endl;
+    CORSIKA_LOG_DEBUG("Decay2()");
     globalCount++;
   }
 
@@ -354,7 +376,10 @@ TEST_CASE("ProcessSequence General", "ProcessSequence") {
     auto sequence2 = make_sequence(cp1, m2, m3);
     GrammageType const tot = sequence2.getInteractionLength(particle);
     InverseGrammageType const tot_inv = sequence2.getInverseInteractionLength(particle);
-    cout << "lambda_tot=" << tot << "; lambda_tot_inv=" << tot_inv << endl;
+    CORSIKA_LOG_DEBUG(
+        "lambda_tot={}"
+        "; lambda_tot_inv={}",
+        tot, tot_inv);
 
     CHECK(tot / 1_g * square(1_cm) == 12);
     CHECK(tot_inv * 1_g / square(1_cm) == 1. / 12);
@@ -373,7 +398,10 @@ TEST_CASE("ProcessSequence General", "ProcessSequence") {
     auto sequence2 = make_sequence(cp1, m2, m3, d3);
     TimeType const tot = sequence2.getLifetime(particle);
     InverseTimeType const tot_inv = sequence2.getInverseLifetime(particle);
-    cout << "lambda_tot=" << tot << "; lambda_tot_inv=" << tot_inv << endl;
+    CORSIKA_LOG_DEBUG(
+        "lambda_tot={}"
+        "; lambda_tot_inv={}",
+        tot, tot_inv);
 
     CHECK(tot / 1_s == 1);
     CHECK(tot_inv * 1_s == 1.);
@@ -419,9 +447,9 @@ TEST_CASE("ProcessSequence General", "ProcessSequence") {
     CORSIKA_LOG_INFO("step2, l={}, i={}", LengthType(step2),
                      ContinuousProcessIndex(step2).getIndex());
 
-    cout << "-->init sequence2" << endl;
+    CORSIKA_LOG_DEBUG("-->init sequence2");
     globalCount = 0;
-    cout << "-->docontinuous" << endl;
+    CORSIKA_LOG_DEBUG("-->docont");
 
     // validation data
     double test_data[nData] = {0};
@@ -431,16 +459,16 @@ TEST_CASE("ProcessSequence General", "ProcessSequence") {
     track = DummyTrajectory();
 
     int const nLoop = 5;
-    cout << "Running loop with n=" << nLoop << endl;
+    CORSIKA_LOG_DEBUG("Running loop with n={}", nLoop);
     for (int iLoop = 0; iLoop < nLoop; ++iLoop) {
       for (int i = 0; i < nData; ++i) { test_data[i] += 0.933 + 0.111; }
       sequence2.doContinuous(particle, track, ContinuousProcessIndex(1));
     }
     for (int i = 0; i < nData; i++) {
-      cout << "data_[" << i << "]=" << particle.data_[i] << endl;
+      CORSIKA_LOG_DEBUG("data_[{}]={}", i, particle.data_[i]);
       CHECK(particle.data_[i] == Approx(test_data[i]).margin(1e-9));
     }
-    cout << "done" << endl;
+    CORSIKA_LOG_DEBUG("done");
   }
 
   SECTION("StackProcess") {
diff --git a/tests/modules/testCONEX.cpp b/tests/modules/testCONEX.cpp
index 621ae626a95b5fe6b1e73953b5d0d94c5a71c4ba..34718d69bb94e93c7e032e42df8f7e8e9c541f6e 100644
--- a/tests/modules/testCONEX.cpp
+++ b/tests/modules/testCONEX.cpp
@@ -97,13 +97,15 @@ TEST_CASE("CONEXSourceCut") {
 
   auto const emPosition = showerCore + showerAxis.getDirection() * (-20_km);
 
-  std::cout << "position injection: "
-            << injectionPos.getCoordinates(conex.getObserverCS()) << " "
-            << injectionPos.getCoordinates(rootCS) << std::endl;
-  std::cout << "position core: " << showerCore.getCoordinates(conex.getObserverCS())
-            << " " << showerCore.getCoordinates(rootCS) << std::endl;
-  std::cout << "position EM: " << emPosition.getCoordinates(conex.getObserverCS()) << " "
-            << emPosition.getCoordinates(rootCS) << std::endl;
+  CORSIKA_LOG_DEBUG("position injection: {} {}",
+                    injectionPos.getCoordinates(conex.getObserverCS()),
+                    injectionPos.getCoordinates(rootCS));
+  CORSIKA_LOG_DEBUG("position core: {} {}",
+                    showerCore.getCoordinates(conex.getObserverCS()),
+                    showerCore.getCoordinates(rootCS));
+  CORSIKA_LOG_DEBUG("position EM: {} {}",
+                    emPosition.getCoordinates(conex.getObserverCS()),
+                    emPosition.getCoordinates(rootCS));
 
   conex.addParticle(Code::Proton, Eem, 0_eV, emPosition, momentum.normalized(), 0_s);
   // supperimpose a photon
diff --git a/tests/modules/testPythia8.cpp b/tests/modules/testPythia8.cpp
index 8ddd613ad21e26f50166a14fcb11a5f7580385cb..8e178affc3f7c26c5a34df8caf2777386dba8d24 100644
--- a/tests/modules/testPythia8.cpp
+++ b/tests/modules/testPythia8.cpp
@@ -27,8 +27,6 @@ TEST_CASE("Pythia", "[processes]") {
 
   SECTION("linking pythia") {
     using namespace Pythia8;
-    using std::cout;
-    using std::endl;
 
     // Generator. Process selection. LHC initialization. Histogram.
     Pythia pythia;
@@ -50,15 +48,15 @@ TEST_CASE("Pythia", "[processes]") {
     event.append(321, 1, 0, 0, 0., 0., 100., sqrt(pz * pz + m * m), m);
 
     if (!pythia.next())
-      cout << "decay failed!" << endl;
+      CORSIKA_LOG_CRITICAL("decay failed!");
     else
-      cout << "particles after decay: " << event.size() << endl;
+      CORSIKA_LOG_DEBUG("particles after decay: {}", event.size());
     event.list();
 
     // loop over final state
     for (int i = 0; i < pythia.event.size(); ++i)
       if (pythia.event[i].isFinal()) {
-        cout << "particle: id=" << pythia.event[i].id() << endl;
+        CORSIKA_LOG_DEBUG("particle: id= {}", pythia.event[i].id());
       }
   }
 
@@ -94,7 +92,7 @@ auto sumMomentum(TStackView const& view, CoordinateSystemPtr const& vCS) {
   return sum;
 }
 
-TEST_CASE("pythia process") {
+TEST_CASE("PythiaInterface", "[processes]") {
 
   logging::set_level(logging::level::info);
   corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
diff --git a/tests/modules/testQGSJetII.cpp b/tests/modules/testQGSJetII.cpp
index 3ddbb254c133050e07b5ec257683767df8508be7..86bb8fce99a1cbed5a57d0eb2018e56278b5c065 100644
--- a/tests/modules/testQGSJetII.cpp
+++ b/tests/modules/testQGSJetII.cpp
@@ -19,7 +19,6 @@
 #include <string>
 #include <cstdlib>
 #include <experimental/filesystem>
-#include <iostream>
 
 using namespace corsika;
 
@@ -49,10 +48,11 @@ TEST_CASE("CORSIKA_DATA", "[processes]") {
     CHECK(data != 0);
     CHECK(std::experimental::filesystem::is_directory(
         std::experimental::filesystem::path(std::string(data) + "/QGSJetII")));
-    std::cout << "data: " << data << " isDir: "
-              << std::experimental::filesystem::is_directory(std::string(data) +
-                                                             "/QGSJetII")
-              << std::endl;
+    CORSIKA_LOG_INFO(
+        "data: {}"
+        " isDir: {}"
+        "/QGSJetII",
+        data, std::experimental::filesystem::is_directory(std::string(data)));
   }
 }