Merge branch '59-provide-central-code-for-lorentz-boosts' into 'master'

Resolve "Provide central code for Lorentz-boosts"

Closes #59

See merge request !38

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

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,

Framework/Utilities/CMakeLists.txt

 
 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)
 

Framework/Utilities/COMBoost.cc

+/**
+ * (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));
+}

Framework/Utilities/COMBoost.h

+/**
+ * (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

Framework/Utilities/Dummy.cc

-#include <corsika/utl/Dummy.h>
-
-using namespace corsika::utl;
-
-// big void...

Framework/Utilities/Dummy.h

-#ifndef _include_corsika_utilties_dummy_h_
-#define _include_corsika_utilties_dummy_h_
-
-namespace corsika::utl {
-
-  // void....
-}
-
-#endif

Framework/Utilities/Test.h

-#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

Framework/Utilities/testBit.cc

-/**
-  \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);

Framework/Utilities/testCOMBoost.cc

+/**
+ * (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));
+}