Added missing classes

corsika/detail/framework/utility/COMBoost.inl

@@ -10,110 +10,111 @@
 
 #pragma once
 
+#include <cmath>
+
 #include <Eigen/Dense>
 #include <corsika/framework/core/PhysicalUnits.hpp>
 #include <corsika/framework/geometry/CoordinateSystem.hpp>
 #include <corsika/framework/geometry/FourVector.hpp>
 #include <corsika/framework/geometry/Vector.hpp>
+#include <corsika/framework/logging/Logging.h>
 #include <corsika/framework/utility/sgn.hpp>
 
+// using namespace corsika::units::si;
+
 namespace corsika {
 
-  auto const& COMBoost::GetRotationMatrix() const { return fRotation; }
+  COMBoost::COMBoost(FourVector<HEPEnergyType, Vector<hepmomentum_d>> const& Pprojectile,
+                     const HEPMassType massTarget)
+      : originalCS_{Pprojectile.GetSpaceLikeComponents().GetCoordinateSystem()}
+      , rotatedCS_{originalCS_.RotateToZ(Pprojectile.GetSpaceLikeComponents())} {
+    auto const pProjectile = Pprojectile.GetSpaceLikeComponents();
+    auto const pProjNormSquared = pProjectile.squaredNorm();
+    auto const pProjNorm = sqrt(pProjNormSquared);
+
+    auto const eProjectile = Pprojectile.GetTimeLikeComponent();
+    auto const massProjectileSquared = eProjectile * eProjectile - pProjNormSquared;
+    auto const s =
+        massTarget * massTarget + massProjectileSquared + 2 * eProjectile * massTarget;
+
+    auto const sqrtS = sqrt(s);
+    auto const sinhEta = -pProjNorm / sqrtS;
+    auto const coshEta = sqrt(1 + pProjNormSquared / s);
+
+    setBoost(coshEta, sinhEta);
+
+    C8LOG_TRACE("COMBoost (1-beta)={}, gamma={}, det={}", 1 - sinhEta / coshEta, coshEta,
+                boost_.determinant() - 1);
+  }
+
+  COMBoost::COMBoost(geometry::Vector<units::si::hepmomentum_d> const& momentum,
+                     units::si::HEPEnergyType mass)
+      : originalCS_{momentum.GetCoordinateSystem()}
+      , rotatedCS_{originalCS_.RotateToZ(momentum)} {
+    auto const squaredNorm = momentum.squaredNorm();
+    auto const norm = sqrt(squaredNorm);
+    auto const sinhEta = -norm / mass;
+    auto const coshEta = sqrt(1 + squaredNorm / (mass * mass));
+    setBoost(coshEta, sinhEta);
+  }
 
-  //! transforms a 4-momentum from lab frame to the center-of-mass frame
   template <typename FourVector>
-  inline FourVector COMBoost::toCoM(const FourVector& p) const {
-    auto pComponents = p.GetSpaceLikeComponents().GetComponents(fCS);
-    Eigen::Vector3d eVecRotated = fRotation * pComponents.eVector;
+  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 = fBoost * lab;
+    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::Vector<hepmomentum_d>(fCS, eVecRotated));
+    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>
-  inline FourVector COMBoost::fromCoM(const FourVector& p) const {
+  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 << (p.GetTimeLikeComponent() * (1 / 1_GeV)),
-        (p.GetSpaceLikeComponents().GetComponents().eVector(2) * (1 / 1_GeV).magnitude());
+    com << (Ecm * (1 / 1_GeV)), (pCM.eVector(2) * (1 / 1_GeV).magnitude());
 
-    auto const plab = p.GetSpaceLikeComponents().GetComponents();
-    std::cout << "COMBoost::fromCoM Ecm=" << p.GetTimeLikeComponent() / 1_GeV << " GeV, "
-              << " pcm = " << plab / 1_GeV << " (norm = " << plab.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 = fInverseBoost * com;
+    auto const boostedZ = inverseBoost_ * com;
     auto const E_lab = boostedZ(0) * 1_GeV;
 
-    auto pLab = p.GetSpaceLikeComponents().GetComponents();
-    pLab.eVector(2) = boostedZ(1) * (1_GeV).magnitude();
-    pLab.eVector = fRotation.transpose() * pLab.eVector;
+    pCM.eVector(2) = boostedZ(1) * (1_GeV).magnitude();
+
+    geometry::Vector<typename decltype(pCM)::dimension> pLab{rotatedCS_, pCM};
+    pLab.rebase(originalCS_);
 
-    FourVector f(E_lab, corsika::Vector(fCS, pLab));
+    FourVector f(E_lab, pLab);
 
-    std::cout << "COMBoost::fromCoM --> Elab=" << E_lab / 1_GeV << "GeV, "
-              << " pcm = " << pLab / 1_GeV << " (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;
   }
 
-  inline 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;
+  void COMBoost::setBoost(double coshEta, double sinhEta) {
+    boost_ << coshEta, sinhEta, sinhEta, coshEta;
+    inverseBoost_ << coshEta, -sinhEta, -sinhEta, coshEta;
   }
 
+  CoordinateSystem const& COMBoost::GetRotatedCS() const { return rotatedCS_; }
+
 } // namespace corsika

corsika/detail/framework/utility/CorsikaData.inl

+/*
+ * (c) Copyright 2020 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 <cstdlib>
+#include <stdexcept>
+#include <string>
+
+std::string corsika::CorsikaData(std::string const& key) {
+  if (auto const* p = std::getenv("CORSIKA_DATA"); p != nullptr) {
+    auto const path = std::string(p) + "/" + key;
+    return path;
+  } else {
+    throw std::runtime_error("CORSIKA_DATA not set");
+  }
+}

corsika/framework/utility/COMBoost.hpp

@@ -8,12 +8,13 @@
 
 #pragma once
 
-#include <iostream>
 
-#include <Eigen/Dense>
-#include <corsika/framework/core/PhysicalUnits.hpp>
 #include <corsika/framework/geometry/CoordinateSystem.hpp>
 #include <corsika/framework/geometry/FourVector.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/logging/Logging.hpp>
+
+#include <Eigen/Dense>
 
 namespace corsika {
 
@@ -23,26 +24,38 @@ namespace corsika {
    */
 
   class COMBoost {
-    Eigen::Matrix3d fRotation;
-    Eigen::Matrix2d fBoost, fInverseBoost;
-    corsika::CoordinateSystem const& fCS;
+    Eigen::Matrix2d boost_, inverseBoost_;
+    corsika::CoordinateSystem const &originalCS_, rotatedCS_;
+
+    void setBoost(double coshEta, double sinhEta);
 
   public:
-    //! construct a COMBoost given four-vector of prjectile and mass of target
-    COMBoost(const corsika::FourVector<HEPEnergyType, corsika::Vector<hepmomentum_d>>&
-                 Pprojectile,
-             const HEPEnergyType massTarget);
+    //! construct a COMBoost given four-vector of projectile and mass of target
+    COMBoost(
+        const corsika::geometry::FourVector<
+            corsika::units::si::HEPEnergyType,
+            corsika::geometry::Vector<corsika::units::si::hepmomentum_d>>& Pprojectile,
+        const corsika::units::si::HEPEnergyType massTarget);
 
-    inline auto const& GetRotationMatrix() const;
+    //! construct a COMBoost to boost into the rest frame given a 3-momentum and mass
+    COMBoost(Vector<units::si::hepmomentum_d> const& momentum,
+             units::si::HEPEnergyType mass);
 
     //! transforms a 4-momentum from lab frame to the center-of-mass frame
     template <typename FourVector>
-    inline FourVector toCoM(const FourVector& p) const;
+    FourVector toCoM(const FourVector& p) const;
+
 
     //! transforms a 4-momentum from the center-of-mass frame back to lab frame
     template <typename FourVector>
-    inline FourVector fromCoM(const FourVector& p) const;
+    FourVector fromCoM(const FourVector& p) const;
+
+
+    CoordinateSystem const& GetRotatedCS() const;
   };
 } // namespace corsika
 
 #include <corsika/detail/framework/utility/COMBoost.inl>
+
+
+

tests/framework/testCorsikaFenv.cpp

@@ -21,7 +21,11 @@ static void handle_fpe(int /*signo*/) { gRESULT = 0; }
 
 TEST_CASE("CorsikaFenv", "[fenv]") {
 
-  feenableexcept(FE_ALL_EXCEPT);
+  SECTION("Enable all exceptions")
+  {
+    feenableexcept(FE_ALL_EXCEPT);
+  }
+  
   signal(SIGFPE, handle_fpe);
 
   SECTION("exception") {
@@ -29,5 +33,8 @@ TEST_CASE("CorsikaFenv", "[fenv]") {
     [[maybe_unused]] auto trigger = std::log(0.);
     std::cout << "trigger: " << trigger << std::endl;
     CHECK(gRESULT == 0);
-  }
+  }  
+  
+}
+
 }

tests/framework/testSaveBoostHistogram.cpp

+/*
+ * (c) Copyright 2020 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 <catch2/catch.hpp>
+#include <corsika/utl/SaveBoostHistogram.hpp>
+
+#include <random>
+
+namespace bh = boost::histogram;
+
+TEST_CASE("SaveHistogram") {
+  std::mt19937 rng;
+  std::normal_distribution n1{2., 6.};
+  std::exponential_distribution e{1.};
+  std::uniform_int_distribution u{1, 10};
+  std::uniform_real_distribution<double> r{-3, 3};
+
+  auto h = bh::make_histogram(
+      bh::axis::regular{5, 0, 10, "normal"}, bh::axis::regular{3, 0, 4, "exponential"},
+      bh::axis::category<int>{{2, 3, 5, 7}, "integer category"},
+      bh::axis::regular<double, bh::use_default, bh::use_default,
+                        bh::axis::option::growth_t>{10, -1, 1, "integer category"});
+
+  for (int i{0}; i < 100'000; ++i) {
+    auto const a = n1(rng);
+    auto const b = e(rng);
+    auto const c = u(rng);
+    auto const d = r(rng);
+
+    h(a, b, c, d);
+  }
+
+  corsika::save_hist(h, "hist.npz");
+}