diff --git a/corsika/detail/framework/utility/COMBoost.inl b/corsika/detail/framework/utility/COMBoost.inl
index 8a87616bc0eb03772c9eb806157cd1d7427d2d7e..352ee4e6c8eee3d7088a2ec0390aa64b16249380 100644
--- a/corsika/detail/framework/utility/COMBoost.inl
+++ b/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
diff --git a/corsika/detail/framework/utility/CorsikaData.inl b/corsika/detail/framework/utility/CorsikaData.inl
new file mode 100644
index 0000000000000000000000000000000000000000..e690e3e0add5dd0327681752a0f2365a6b299585
--- /dev/null
+++ b/corsika/detail/framework/utility/CorsikaData.inl
@@ -0,0 +1,22 @@
+/*
+ * (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");
+ }
+}
diff --git a/corsika/framework/utility/COMBoost.hpp b/corsika/framework/utility/COMBoost.hpp
index 96823d123fb6fb2f9c5ec0c44415dab9ec277fbd..08932465b1f54792126179d26564dfd04ccff8ef 100644
--- a/corsika/framework/utility/COMBoost.hpp
+++ b/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>
+
+
+
diff --git a/tests/framework/testCorsikaFenv.cpp b/tests/framework/testCorsikaFenv.cpp
index e4ec2eff44264cb30a33780bd1e7c387f6cac60f..a82f440333fdc5b5ccc78f38b93dcdf0b460ecba 100644
--- a/tests/framework/testCorsikaFenv.cpp
+++ b/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);
- }
+ }
+
+}
+
}
diff --git a/tests/framework/testSaveBoostHistogram.cpp b/tests/framework/testSaveBoostHistogram.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..b825816fe0af8ca112db1999e450d094dd430532
--- /dev/null
+++ b/tests/framework/testSaveBoostHistogram.cpp
@@ -0,0 +1,39 @@
+/*
+ * (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");
+}