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/CMakeList.txt b/Framework/Utilities/CMakeList.txt
deleted file mode 100644
index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..0000000000000000000000000000000000000000
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/Test.h b/Framework/Utilities/Test.h
deleted file mode 100644
index 41c7cecd7a9856e2c5b68f845d22a25f64a3cf0a..0000000000000000000000000000000000000000
--- a/Framework/Utilities/Test.h
+++ /dev/null
@@ -1,237 +0,0 @@
-#ifndef _utl_Test_h_
-#define _utl_Test_h_
-
-/**
- \file
- Tools to do simple testing in a readable way
-
- \author Lukas Nellen
- \author Darko Veberic
- \version $Id: Test.h 31925 2018-09-25 16:02:12Z darko $
- \date 08 Feb 2004
-
- \ingroup testing
-*/
-
-#include <utl/Triple.h>
-#include <boost/format.hpp>
-#include <boost/tuple/tuple.hpp>
-#include <boost/tuple/tuple_comparison.hpp>
-#include <boost/tuple/tuple_io.hpp>
-#include <cmath>
-
-namespace utl {
-
- //! Predicate used in STL for searching for whitespace
- struct IsSpace {
- bool operator()(const char x) const {
- return x == ' ' || x == '\r' || x == '\n' || x == '\t';
- }
- };
-
- /// Predicate for equality
- class Equal {
- public:
- template <typename T>
- bool operator()(const T& lhs, const T& rhs) const {
- return lhs == rhs;
- }
-
- static const char* Name() { return "equal"; }
- };
-
- /// Predicate for less
- class Less {
- public:
- template <typename T>
- bool operator()(const T& lhs, const T& rhs) const {
- return lhs < rhs;
- }
-
- static const char* Name() { return "less"; }
- };
-
- /// Predicate for less or equal
- class LessOrEqual {
- public:
- template <typename T>
- bool operator()(const T& lhs, const T& rhs) const {
- return lhs <= rhs;
- }
-
- static const char* Name() { return "less or equal"; }
- };
-
- /// Predicate for greater
- class Greater {
- public:
- template <typename T>
- bool operator()(const T& lhs, const T& rhs) const {
- return lhs > rhs;
- }
-
- static const char* Name() { return "greater"; }
- };
-
- /// Predicate for greater or equal
- class GreaterOrEqual {
- public:
- template <typename T>
- bool operator()(const T& lhs, const T& rhs) const {
- return lhs >= rhs;
- }
-
- static const char* Name() { return "greater or equal"; }
- };
-
- /// Predicate for approximate equality (for floating point)
- /** The default precision is 1e-6, but it can be changed at
- construction time.
- */
- class CloseTo {
- public:
- CloseTo(const double eps = 1e-6)
- : fEpsilon(eps) {}
-
- template <typename T>
- bool operator()(const T& lhs, const T& rhs) const {
- return IsCloseTo(lhs, rhs);
- }
-
- boost::format Name() const { return boost::format("close (@%g) to") % fEpsilon; }
-
- protected:
- template <typename T>
- bool IsCloseAbs(const T& lhs, const T& rhs) const {
- return std::abs(double(lhs) - double(rhs)) < fEpsilon;
- }
-
- bool IsCloseAbs(const utl::Triple& lhs, const utl::Triple& rhs) const {
- return std::sqrt(TupleDist2(lhs, rhs)) < fEpsilon;
- }
-
- template <typename T>
- bool IsCloseRel(const T& lhs, const T& rhs) const {
- return 2 * std::abs(double(lhs) - double(rhs)) /
- (std::abs(double(lhs)) + std::abs(double(rhs))) <
- fEpsilon;
- }
-
- bool IsCloseRel(const utl::Triple& lhs, const utl::Triple& rhs) const {
- return (2 * sqrt(TupleDist2(lhs, rhs)) /
- (sqrt(TupleDist2(lhs, utl::Triple(0, 0, 0))) +
- sqrt(TupleDist2(rhs, utl::Triple(0, 0, 0))))) < fEpsilon;
- }
-
- template <typename T>
- bool IsCloseTo(const T& lhs, const T& rhs) const {
- if (IsCloseAbs(lhs, rhs))
- return true;
- else
- return IsCloseRel(lhs, rhs);
- }
-
- // tuple distance
- template <typename Head, typename Tail>
- static double TupleDist2(const boost::tuples::cons<Head, Tail>& lhs,
- const boost::tuples::cons<Head, Tail>& rhs) {
- const double t = lhs.get_head() - rhs.get_head();
- return t * t + TupleDist2(lhs.get_tail(), rhs.get_tail());
- }
-
- static double TupleDist2(const boost::tuples::null_type& /*lhs*/,
- const boost::tuples::null_type& /*rhs*/) {
- return 0;
- }
-
- double fEpsilon;
- };
-
- class CloseAbs : public CloseTo {
- public:
- CloseAbs(const double eps = 1e-6)
- : CloseTo(eps) {}
-
- template <typename T>
- bool operator()(const T& lhs, const T& rhs) const {
- return IsCloseAbs(lhs, rhs);
- }
-
- boost::format Name() const {
- return boost::format("absolutely close (@%g) to") % fEpsilon;
- }
- };
-
- class CloseRel : public CloseTo {
- public:
- CloseRel(const double eps = 1e-6)
- : CloseTo(eps) {}
-
- template <typename T>
- bool operator()(const T& lhs, const T& rhs) const {
- return IsCloseRel(lhs, rhs);
- }
-
- boost::format Name() const {
- return boost::format("relatively close (@%g) to") % fEpsilon;
- }
- };
-
- template <typename Predicate>
- class Not : public Predicate {
- public:
- Not()
- : Predicate() {}
-
- Not(const double eps)
- : Predicate(eps) {}
-
- template <typename T>
- bool operator()(const T& x) const {
- return !Predicate::operator()(x);
- }
-
- template <typename T, typename U>
- bool operator()(const T& x, const U& y) const {
- return !Predicate::operator()(x, y);
- }
-
- template <typename T, typename U, typename W>
- bool operator()(const T& x, const U& y, const W& z) const {
- return !Predicate::operator()(x, y, z);
- }
-
- static boost::format Name() { return boost::format("not-%s") % Predicate().Name(); }
- };
-
- inline utl::Triple Diff(const utl::Triple& lhs, const utl::Triple& rhs) {
- return utl::Triple(lhs.get<0>() - rhs.get<0>(), lhs.get<1>() - rhs.get<1>(),
- lhs.get<2>() - rhs.get<2>());
- }
-
- /// Test condition by evaluating a predicate
- /** If the predicate evaluates to false, we print a failure message
- with the values of the left- and right-hand side and the name of
- the predicate. This information is normally not available when
- using the CPPUNIT_ASSERT macro.
- */
- template <class Predicate, typename T>
- inline bool Test(const Predicate& pred, const T& lhs, const T& rhs) {
- return pred(lhs, rhs);
- }
-
- /// Main test function
- template <class Predicate, typename T>
- inline bool Test(const T& lhs, const T& rhs) {
- return Test(Predicate(), lhs, rhs);
- }
-
- /// Test function for predicates that take an option
- template <class Predicate, typename T, typename U>
- inline bool Test(const T& lhs, const T& rhs, const U& eps) {
- return Test(Predicate(eps), lhs, rhs);
- }
-
-} // namespace utl
-
-#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..814d63e1f2f44463e2861f08d4b6ae5d3730af8e
--- /dev/null
+++ b/Framework/Utilities/testCOMBoost.cc
@@ -0,0 +1,81 @@
+/**
+ * (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-6;
+
+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();
+ };
+
+ // define projectile kinematics in lab frame
+ MassType const projectileMass = 1._GeV / cSquared;
+ Vector<momentum_d> pProjectileLab{rootCS, {0_GeV / c, 1_PeV / c, 0_GeV / c}};
+ EnergyType const eProjectileLab = energy(projectileMass, pProjectileLab);
+
+ // define target kinematics in lab frame
+ MassType const targetMass = 1_GeV / cSquared;
+ Vector<momentum_d> pTargetLab{rootCS, {0_Ns, 0_Ns, 0_Ns}};
+ EnergyType const eTargetLab = energy(targetMass, pTargetLab);
+
+ // define boost to com frame
+ COMBoost boost(eProjectileLab, pProjectileLab, targetMass);
+
+ // boost projecticle
+ auto const [eProjectileCoM, pProjectileCoM] =
+ boost.toCoM(eProjectileLab, pProjectileLab);
+
+ // boost target
+ auto const [eTargetCoM, pTargetCoM] = boost.toCoM(eTargetLab, pTargetLab);
+
+ // sum of momenta in CoM, should be 0
+ auto const sumPCoM = pProjectileCoM + pTargetCoM;
+ CHECK(sumPCoM[2] / (1_GeV / c) == Approx(0).margin(absMargin));
+
+ // mandelstam-s should be invariant under transformation
+ 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)));
+
+ // boost back...
+ auto const [eProjectileBack, pProjectileBack] =
+ boost.fromCoM(eProjectileCoM, pProjectileCoM);
+
+ // ...should yield original values before the boosts
+ CHECK(eProjectileBack / eProjectileLab == Approx(1));
+ CHECK((pProjectileBack - pProjectileLab).norm() / pProjectileLab.norm() ==
+ Approx(0).margin(absMargin));
+}