diff --git a/Framework/Geometry/BaseVector.h b/Framework/Geometry/BaseVector.h
index 2cd1c3e490be90485164a060cefc708dbd64806e..c16e83bb87c2d2677e87b2b11b06bd1aa958be22 100644
--- a/Framework/Geometry/BaseVector.h
+++ b/Framework/Geometry/BaseVector.h
@@ -31,6 +31,10 @@ namespace corsika::geometry {
BaseVector(CoordinateSystem const& pCS, QuantityVector<dim> pQVector)
: qVector(pQVector)
, cs(&pCS) {}
+
+ auto const& GetCoordinateSystem() const {
+ return *cs;
+ }
};
} // namespace corsika::geometry
diff --git a/Framework/Geometry/CoordinateSystem.h b/Framework/Geometry/CoordinateSystem.h
index 6f0584f5e015aa521a0031919e0b0dda98d8d224..49deebc648bbca7a3fef5d7f8137c499c225f6fa 100644
--- a/Framework/Geometry/CoordinateSystem.h
+++ b/Framework/Geometry/CoordinateSystem.h
@@ -41,7 +41,7 @@ namespace corsika::geometry {
protected:
static auto CreateCS() { return CoordinateSystem(); }
friend corsika::geometry::RootCoordinateSystem; /// this is the only class that can
- /// creat ONE unique root CS
+ /// create ONE unique root CS
public:
static EigenTransform GetTransformation(CoordinateSystem const& c1,
diff --git a/Framework/Utilities/CMakeLists.txt b/Framework/Utilities/CMakeLists.txt
index ab6c148d9311817456198282eb1757a25264952e..5ba7827be04ef4320da6374a5b52f0fbce79be13 100644
--- a/Framework/Utilities/CMakeLists.txt
+++ b/Framework/Utilities/CMakeLists.txt
@@ -1,14 +1,13 @@
set (
UTILITIES_SOURCES
- Dummy.cc
+ COMBoost.cc
)
set (
UTILITIES_HEADERS
- Dummy.h
+ COMBoost.h
Bit.h
- Test.h
Singleton.h
)
@@ -29,8 +28,11 @@ set_target_properties (
)
# target dependencies on other libraries (also the header onlys)
-#target_link_libraries (
-# )
+target_link_libraries (
+ CORSIKAutilities
+ CORSIKAgeometry
+ CORSIKAunits
+ )
target_include_directories (
CORSIKAutilities
@@ -49,13 +51,13 @@ install (
# --------------------
# code unit testing
-# add_executable (testBit testBit.cc)
+add_executable (testCOMBoost testCOMBoost.cc)
-# target_link_libraries (
-# testBit
-# CORSIKAutilities
-# CORSIKAthirdparty # for catch2
-# )
+target_link_libraries (
+ testCOMBoost
+ CORSIKAutilities
+ CORSIKAthirdparty # for catch2
+ )
-# add_test (NAME testBit COMMAND testBit)
+add_test (NAME testCOMBoost COMMAND testCOMBoost)
diff --git a/Framework/Utilities/COMBoost.cc b/Framework/Utilities/COMBoost.cc
new file mode 100644
index 0000000000000000000000000000000000000000..4e8a9688ea07cb6ccd8a09e3f49d3f525effc360
--- /dev/null
+++ b/Framework/Utilities/COMBoost.cc
@@ -0,0 +1,78 @@
+/**
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * 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/utl/COMBoost.h>
+
+using namespace corsika::utl;
+using namespace corsika::units::si;
+
+COMBoost::COMBoost(EnergyType eProjectile, COMBoost::MomentumVector const& pProjectile,
+ MassType mTarget)
+ : fRotation(Eigen::Matrix3d::Identity())
+ , fCS(pProjectile.GetCoordinateSystem()) {
+ // calculate matrix for rotating pProjectile to z-axis first
+ // TODO: handle the case when pProjectile ~ (0, 0, -1)
+ auto const pProjNorm = pProjectile.norm();
+ auto const a = (pProjectile / pProjNorm).GetComponents().eVector;
+
+ Eigen::Vector3d const b{0, 0, 1};
+ auto const v = a.cross(b);
+
+ Eigen::Matrix3d vHat;
+ vHat << 0, -v(2), v(1), v(2), 0, -v(0), -v(1), v(0), 0;
+
+ fRotation += vHat + vHat * vHat / (1 + a.dot(b));
+
+ // calculate boost
+ double const x = pProjNorm * units::constants::c /
+ (eProjectile + mTarget * units::constants::cSquared);
+
+ /* Accurracy matters here, x = 1 - epsilon for ultra-relativistic boosts */
+ double const coshEta = 1 / std::sqrt((1 + x) * (1 - x));
+ //~ double const coshEta = 1 / std::sqrt((1-x*x));
+ double const sinhEta = -x * coshEta;
+
+ fBoost << coshEta, sinhEta, sinhEta, coshEta;
+
+ fInverseBoost << coshEta, -sinhEta, -sinhEta, coshEta;
+}
+
+std::tuple<EnergyType, corsika::geometry::QuantityVector<momentum_d>> COMBoost::toCoM(
+ EnergyType E, COMBoost::MomentumVector p) const {
+ corsika::geometry::QuantityVector<momentum_d> pComponents = p.GetComponents(fCS);
+ Eigen::Vector3d eVecRotated = fRotation * pComponents.eVector;
+ Eigen::Vector2d lab;
+
+ lab << (E * (1 / 1_GeV)), (eVecRotated(2) * (units::constants::c / 1_GeV).magnitude());
+
+ auto const boostedZ = fBoost * lab;
+ auto const E_CoM = boostedZ(0) * 1_GeV;
+
+ eVecRotated(2) = boostedZ(1) * (1_GeV / units::constants::c).magnitude();
+
+ return std::make_tuple(E_CoM,
+ corsika::geometry::QuantityVector<momentum_d>{eVecRotated});
+}
+
+std::tuple<EnergyType, COMBoost::MomentumVector> COMBoost::fromCoM(
+ EnergyType E, corsika::geometry::QuantityVector<units::si::momentum_d> pCoM) const {
+ Eigen::Vector2d com;
+ com << (E * (1 / (units::constants::c * 1_Ns))),
+ (pCoM.eVector(2) * (1 / 1_Ns).magnitude());
+
+ auto const boostedZ = fInverseBoost * com;
+ auto const E_CoM = boostedZ(0) * (1_Ns * units::constants::c);
+
+ auto pLab = pCoM;
+ pLab.eVector(2) = boostedZ(1) * (1_Ns).magnitude();
+ pLab.eVector = fRotation.transpose() * pLab.eVector;
+
+ return std::make_tuple(E_CoM, MomentumVector(fCS, pLab));
+}
diff --git a/Framework/Utilities/COMBoost.h b/Framework/Utilities/COMBoost.h
new file mode 100644
index 0000000000000000000000000000000000000000..6250d422dc256c4cad102bb93bd596e9a0df2a30
--- /dev/null
+++ b/Framework/Utilities/COMBoost.h
@@ -0,0 +1,44 @@
+/**
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * 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.
+ */
+
+#ifndef _include_corsika_utilties_comboost_h_
+#define _include_corsika_utilties_comboost_h_
+
+#include <corsika/geometry/CoordinateSystem.h>
+#include <corsika/geometry/QuantityVector.h>
+#include <corsika/geometry/Vector.h>
+#include <corsika/units/PhysicalUnits.h>
+#include <Eigen/Dense>
+#include <tuple>
+
+namespace corsika::utl {
+ class COMBoost {
+ Eigen::Matrix3d fRotation;
+ Eigen::Matrix2d fBoost, fInverseBoost;
+ corsika::geometry::CoordinateSystem const& fCS;
+ using MomentumVector = corsika::geometry::Vector<units::si::momentum_d>;
+
+ public:
+ //! construct a COMBoost given energy and momentum of projectile and mass of target
+ COMBoost(units::si::EnergyType eProjectile, MomentumVector const& pProjectile,
+ units::si::MassType mTarget);
+
+ //! transforms a 4-momentum from lab frame to the center-of-mass frame
+ std::tuple<units::si::EnergyType, geometry::QuantityVector<units::si::momentum_d>>
+ toCoM(units::si::EnergyType E, MomentumVector p) const;
+
+ //! transforms a 4-momentum from the center-of-mass frame back to lab frame
+ std::tuple<units::si::EnergyType, MomentumVector> fromCoM(
+ units::si::EnergyType E,
+ geometry::QuantityVector<units::si::momentum_d> pCoM) const;
+ };
+} // namespace corsika::utl
+
+#endif
diff --git a/Framework/Utilities/Dummy.cc b/Framework/Utilities/Dummy.cc
deleted file mode 100644
index fd578b8e09c10e3b26932dbc1dfdc1c36a1030b9..0000000000000000000000000000000000000000
--- a/Framework/Utilities/Dummy.cc
+++ /dev/null
@@ -1,5 +0,0 @@
-#include <corsika/utl/Dummy.h>
-
-using namespace corsika::utl;
-
-// big void...
diff --git a/Framework/Utilities/Dummy.h b/Framework/Utilities/Dummy.h
deleted file mode 100644
index dac2021b49df401d221c5160f825909440468e22..0000000000000000000000000000000000000000
--- a/Framework/Utilities/Dummy.h
+++ /dev/null
@@ -1,9 +0,0 @@
-#ifndef _include_corsika_utilties_dummy_h_
-#define _include_corsika_utilties_dummy_h_
-
-namespace corsika::utl {
-
- // void....
-}
-
-#endif
diff --git a/Framework/Utilities/testBit.cc b/Framework/Utilities/testBit.cc
deleted file mode 100644
index a4d36d1bcae35a23c3baaa0fa96a0c05821c8fa5..0000000000000000000000000000000000000000
--- a/Framework/Utilities/testBit.cc
+++ /dev/null
@@ -1,90 +0,0 @@
-/**
- \file
- Test Bit functions
-
- \author Hans Dembinski
- \version $Id: testBit.cc 25126 2014-02-03 22:13:10Z darko $
- \date 27 Jan 2014
-
- \ingroup testing
-*/
-
-#include <corsika/utl/Test.h>
-#include <cppunit/extensions/HelperMacros.h>
-#include <tst/Verify.h>
-#include <utl/Bit.h>
-#include <bitset>
-#include <cstdio>
-#include <iostream>
-
-using namespace tst;
-using namespace utl;
-using namespace std;
-
-/**
- \ingroup testing
-*/
-class TestBit : public CppUnit::TestFixture {
-
- CPPUNIT_TEST_SUITE(TestBit);
- CPPUNIT_TEST(TestGet);
- CPPUNIT_TEST(TestSet);
- CPPUNIT_TEST(TestMask);
- CPPUNIT_TEST_SUITE_END();
-
-public:
- void setUp() {}
-
- void tearDown() {}
-
- void TestGet() {
- const int size = sizeof(int) * 8;
- const int bc2 = 12345;
- int b2 = bc2;
- bitset<size> b1(bc2);
-
- ostringstream out1;
- ostringstream out2;
- ostringstream out3;
- for (int i = 0; i < size; ++i) {
- out1 << (b1[i] ? '^' : '.');
- out2 << (AsBitArray(bc2)[i] ? '^' : '.');
- out3 << (AsBitArray(b2)[i] ? '^' : '.');
- }
-
- CPPUNIT_ASSERT(Verify<Equal>(out1.str(), out2.str()));
- CPPUNIT_ASSERT(Verify<Equal>(out1.str(), out3.str()));
- }
-
- void TestSet() {
- const int size = sizeof(int) * 8;
- const int number = 12345;
- bitset<size> b1(number);
- int b2 = 11111;
-
- for (int i = 0; i < size; ++i) AsBitArray(b2)[i] = b1[i];
-
- CPPUNIT_ASSERT(Verify<Equal>(b2, number));
- }
-
- void TestMask() {
- const int n = (1 << 18) | (1 << 5);
- int m = 0;
-
- AsBitArray(m)[18] = true;
- AsBitArray(m)[5] = true;
- CPPUNIT_ASSERT(Verify<Equal>(n, m));
-
- for (unsigned int i = 0; i < 8 * sizeof(int); ++i) AsBitArray(m)[i] = 0;
- CPPUNIT_ASSERT(Verify<Equal>(m, 0));
-
- m = 1;
- AsBitArray(m).Mask(n, true);
- CPPUNIT_ASSERT(Verify<Equal>(m, n + 1));
-
- AsBitArray(m).Mask(n, false);
- CPPUNIT_ASSERT(Verify<Equal>(m, 1));
- }
-};
-
-CPPUNIT_TEST_SUITE_REGISTRATION(TestBit);
diff --git a/Framework/Utilities/testCOMBoost.cc b/Framework/Utilities/testCOMBoost.cc
new file mode 100644
index 0000000000000000000000000000000000000000..353fdaeee1bab24adba1e6bca62d489e8c1cff8c
--- /dev/null
+++ b/Framework/Utilities/testCOMBoost.cc
@@ -0,0 +1,72 @@
+/**
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * 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.
+ */
+
+#define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one
+ // cpp file
+#include <catch2/catch.hpp>
+
+#include <corsika/geometry/RootCoordinateSystem.h>
+#include <corsika/units/PhysicalUnits.h>
+#include <corsika/utl/COMBoost.h>
+#include <iostream>
+
+using namespace corsika::geometry;
+using namespace corsika::utl;
+using namespace corsika::units::si;
+using corsika::units::constants::c;
+using corsika::units::constants::cSquared;
+
+double constexpr absMargin = 1e-8;
+
+TEST_CASE("boosts") {
+ CoordinateSystem& rootCS =
+ RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
+
+ // relativistic energy
+ auto energy = [](MassType m, Vector<momentum_d> const& p) {
+ return sqrt(m * m * cSquared * cSquared + p.squaredNorm() * cSquared);
+ };
+
+ // mandelstam-s
+ auto s = [](EnergyType E, QuantityVector<momentum_d> const& p) {
+ return E * E / cSquared - p.squaredNorm();
+ };
+
+ MassType const projectileMass = 4_GeV / cSquared;
+ Vector<momentum_d> pProjectileLab{rootCS, {0_GeV / c, 100_GeV / c, 0_GeV / c}};
+ EnergyType const eProjectileLab = energy(projectileMass, pProjectileLab);
+
+ MassType const targetMass = 1_GeV / cSquared;
+ Vector<momentum_d> pTargetLab{rootCS, {0_Ns, 0_Ns, 0_Ns}};
+ EnergyType const eTargetLab = energy(targetMass, pTargetLab);
+
+ COMBoost boost(eProjectileLab, pProjectileLab, targetMass);
+
+ auto const [eProjectileCoM, pProjectileCoM] =
+ boost.toCoM(eProjectileLab, pProjectileLab);
+ auto const [eTargetCoM, pTargetCoM] = boost.toCoM(eTargetLab, pTargetLab);
+
+ auto const sumPCoM = pProjectileCoM + pTargetCoM;
+
+ CHECK(s(eProjectileLab + eTargetLab,
+ pProjectileLab.GetComponents() + pTargetLab.GetComponents()) /
+ (1_GeV / c) / (1_GeV / c) ==
+ Approx(s(eProjectileCoM + eTargetCoM, pProjectileCoM + pTargetCoM) / (1_GeV / c) /
+ (1_GeV / c)));
+
+ CHECK(sumPCoM[2] / (1_GeV / c) == Approx(0).margin(absMargin));
+
+ auto const [eProjectileBack, pProjectileBack] =
+ boost.fromCoM(eProjectileCoM, pProjectileCoM);
+
+ CHECK(eProjectileBack / eProjectileLab == Approx(1));
+ CHECK((pProjectileBack - pProjectileLab).norm() / (1_GeV / c) ==
+ Approx(0).margin(absMargin));
+}