implemented Lorentz transform, accuracy not so good for huge E_proj/m_target

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

[Ralf Ulrich]

We should use corsika::units::hep for all energy-momentum operations. This is also done on the stack now.

[Maximilian Reininghaus]

I don't like your units::hep solution. It does not really prevent mixing with the SI units. I think we should rather add another dimension for microscopic energy/momentum/mass in PhysUnits, even if it is redundant in principle. We rarely want to cancel eV with macroscopic units.

[Ralf Ulrich]

True, mixing with SI units is absolutely possible right now, this is slightly problematic. At least in typical application this is ok.

Indeed, we could add a complete new dimension to phys/units (I think then we finally take over full control and responsibility over this package). That should not difficult to add. And for user code this will almost make no difference: we just use units::hep like now, there only will be extra conversion factors between units::hep <-> units::si, we might even define them already now.

[Felix Riehn]

So whats the verdict? I just ran into this problem trying to use the new COMBoost for the sibyll stack (which is on hep units). Do I convert in place or are you planning to change the boost to hep soon?

[Maximilian Reininghaus]

I think for the time being you can just remove all the annoying c factors in COMBoost and switch to units::hep there. I will look into the units issue again next year.

[Felix Riehn]

OK. will give it a try.