added FourVector with tests

Framework/Geometry/FourVector.h

+#ifndef _include_corsika_framework_geometry_fourvector_h_
+#define _include_corsika_framework_geometry_fourvector_h_
+
+#include <corsika/geometry/Vector.h>
+#include <corsika/units/PhysicalUnits.h>
+
+#include <iostream>
+#include <type_traits>
+
+namespace corsika::geometry {
+
+/**
+   FourVector supports "full" units, e.g. E in [GeV/c] and p in [GeV],
+   or also t in [s] and r in [m], etc.
+
+   However, for HEP applications it is also possible to use E and p
+   both in [GeV].
+
+   The FourVector can return NormSqr and Norm, whereas Norm is
+   sqrt(abs(NormSqr)). The physical units are always calculated and
+   returned properly.
+
+   FourVector can also return if it is TimeLike, SpaceLike or PhotonLike.
+
+   When a FourVector is initialized with a lvalue reference, this is
+   also used for the internal storage, which should lead to complete
+   disappearance of the FourVector class during optimization.
+ */
+
+template <typename TimeType, typename SpaceVec> class FourVector {
+
+public:
+  using SpaceType = typename std::decay<SpaceVec>::type::Quantity;
+
+  /// check the types and the physical units here:
+  static_assert(
+      std::is_same<typename std::decay<TimeType>::type, SpaceType>::value ||
+          std::is_same<typename std::decay<TimeType>::type,
+                       decltype(std::declval<SpaceType>() /
+                                corsika::units::si::meter *
+                                corsika::units::si::second)>::value,
+      "Units of time-like and space-like coordinates must either be idential "
+      "(e.g. GeV) or [E/c]=[p]");
+
+public:
+  /*
+  template <typename TT, typename SS>
+  FourVector(TT && eT, SS && eS)
+    : fTimeLike(std::forward<TT>(eT)), fSpaceLike(std::forward<SS>(eS)) {
+    std::cout << "FourVector&&\n"; }
+  */
+
+  FourVector(const TimeType &eT, const SpaceVec &eS)
+      : fTimeLike(eT), fSpaceLike(eS) {
+    // std::cout << "FourVector const &\n";
+  }
+
+  /*
+ FourVector(TimeType &eT, SpaceVec &eS)
+      : fTimeLike(eT), fSpaceLike(eS) {
+    std::cout << "FourVector &\n"; }
+  */
+
+  TimeType GetTime() { return fTimeLike; }
+
+  auto GetNormSqr() const {
+    return GetTimeSquared() - fSpaceLike.squaredNorm();
+  }
+
+  SpaceType GetNorm() const { return sqrt(abs(GetNormSqr())); }
+
+  bool IsTimelike() const {
+    return GetTimeSquared() < fSpaceLike.squaredNorm();
+  } // Norm2 < 0
+
+  bool IsSpacelike() const {
+    return GetTimeSquared() > fSpaceLike.squaredNorm();
+  } // Norm2 > 0
+
+  bool IsPhotonlike() const {
+    return GetTimeSquared() == fSpaceLike.squaredNorm();
+  } // // Norm2 == 0
+
+  FourVector &operator+=(const FourVector &b) {
+    fTimeLike += b.fTimeLike;
+    fSpaceLike += b.fSpaceLike;
+    return *this;
+  }
+
+  FourVector &operator-=(const FourVector &b) {
+    fTimeLike -= b.fTimeLike;
+    fSpaceLike -= b.fSpaceLike;
+    return *this;
+  }
+
+  FourVector &operator*=(const double b) {
+    fTimeLike *= b;
+    fSpaceLike *= b;
+    return *this;
+  }
+
+  FourVector &operator/=(const double b) {
+    fTimeLike /= b;
+    fSpaceLike.GetComponents() /= b; // TODO: WHY IS THIS??????
+    return *this;
+  }
+
+  FourVector &operator/(const double b) {
+    *this /= b;
+    return *this;
+  }
+
+  /**
+     Note that the product between two 4-vectors assumes that you use
+     the same "c" convention for both. Only the LHS vector is checked
+     for this. You cannot mix different conventions due to
+     unit-checking.
+   */
+  SpaceType operator*(const FourVector &b) {
+    if constexpr (std::is_same<typename std::decay<TimeType>::type,
+                               decltype(std::declval<SpaceType>() /
+                                        corsika::units::si::meter *
+                                        corsika::units::si::second)>::value)
+      return fTimeLike * b.fTimeLike *
+                 (corsika::units::constants::c * corsika::units::constants::c) -
+             fSpaceLike.norm();
+    else
+      return fTimeLike * fTimeLike - fSpaceLike.norm();
+  }
+
+private:
+  /**
+     This function is automatically compiled to use of ignore the
+     extra factor of "c" for the time-like quantity
+   */
+  auto GetTimeSquared() const {
+    if constexpr (std::is_same<typename std::decay<TimeType>::type,
+                               decltype(std::declval<SpaceType>() /
+                                        corsika::units::si::meter *
+                                        corsika::units::si::second)>::value)
+      return fTimeLike * fTimeLike *
+             (corsika::units::constants::c * corsika::units::constants::c);
+    else
+      return fTimeLike * fTimeLike;
+  }
+
+protected:
+  /// the data members
+  TimeType fTimeLike;
+  SpaceVec fSpaceLike;
+
+  
+  /// the friends: math operators
+  template <typename T, typename U>
+  friend FourVector<T, U> operator+(const FourVector<T, U> &,
+                                    const FourVector<T, U> &);
+
+  template <typename T, typename U>
+  friend FourVector<T, U> operator-(const FourVector<T, U> &,
+                                    const FourVector<T, U> &);
+
+  template <typename T, typename U>
+  friend FourVector<T, U> operator*(const FourVector<T, U> &, const double);
+
+    template <typename T, typename U>
+  friend FourVector<T, U> operator/(const FourVector<T, U> &, const double);
+
+};
+
+/*
+//template<typename T, typename U> FourVector(T& t, U& u) ->
+FourVector<decltype(t), decltype(u)>; template<typename T, typename U>
+FourVector(const T& t, const U& u) -> FourVector<const typename
+std::decay<T>::type, const typename std::decay<U>::type>; template<typename T,
+typename U> FourVector(T&& t, U&& u) -> FourVector<typename std::decay<T>::type,
+typename std::decay<U>::type>;
+ // template<typename T, typename U> FourVector(T&& t, U&& u) ->
+FourVector<decltype(t), decltype(u)>;
+ */
+
+/** 
+    The math operator+
+ */
+template <typename TimeType, typename SpaceVec>
+inline FourVector<TimeType, SpaceVec>
+operator+(const FourVector<TimeType, SpaceVec> &a,
+          const FourVector<TimeType, SpaceVec> &b) {
+  return FourVector<TimeType, SpaceVec>(a.fTimeLike + b.fTimeLike,
+                                        a.fSpaceLike + b.fSpaceLike);
+}
+
+ /** 
+    The math operator-
+ */
+template <typename TimeType, typename SpaceVec>
+inline FourVector<TimeType, SpaceVec>
+operator-(const FourVector<TimeType, SpaceVec> &a,
+          const FourVector<TimeType, SpaceVec> &b) {
+  return FourVector<TimeType, SpaceVec>(a.fTimeLike - b.fTimeLike,
+                                        a.fSpaceLike - b.fSpaceLike);
+}
+
+ /** 
+    The math operator*
+ */
+template <typename TimeType, typename SpaceVec>
+inline FourVector<TimeType, SpaceVec>
+operator*(const FourVector<TimeType, SpaceVec> &a, const double b) {
+  return FourVector<TimeType, SpaceVec>(a.fTimeLike * b, a.fSpaceLike * b);
+}
+
+/** 
+    The math operator/
+ */
+ template <typename TimeType, typename SpaceVec>
+inline FourVector<TimeType, SpaceVec>
+operator/(const FourVector<TimeType, SpaceVec> &a, const double b) {
+  return FourVector<TimeType, SpaceVec>(a.fTimeLike / b, a.fSpaceLike / b);
+}
+
+} // namespace corsika::geometry
+
+#endif

Framework/Geometry/testFourVector.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/CoordinateSystem.h>
+#include <corsika/geometry/FourVector.h>
+#include <corsika/geometry/RootCoordinateSystem.h>
+#include <corsika/geometry/Vector.h>
+#include <corsika/units/PhysicalUnits.h>
+#include <cmath>
+
+#include <boost/type_index.hpp>
+using boost::typeindex::type_id_with_cvr;
+
+using namespace corsika::geometry;
+using namespace corsika::units::si;
+
+
+TEST_CASE("four vectors") {
+
+  // this is just needed as a baseline
+  CoordinateSystem& rootCS =
+      RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
+
+  /*
+    Test: P2 = E2 - p2 all in [GeV]
+    This is the typical HEP application
+   */
+  SECTION("Energy momentum in hep-units") {
+
+    HEPEnergyType E0 = 10_GeV;
+    Vector<hepmomentum_d> P0(rootCS, {10_GeV, 10_GeV, 10_GeV});
+
+    FourVector p0(E0, P0);
+
+    REQUIRE(p0.GetNormSqr() == -200_GeV * 1_GeV);
+    REQUIRE(p0.GetNorm() == sqrt(200_GeV * 1_GeV));
+  }
+
+  /*
+    Check space/time-like
+   */
+  SECTION("Space/time likeness") {
+
+    HEPEnergyType E0 = 20_GeV;
+    Vector<hepmomentum_d> P0(rootCS, {10_GeV, 0_GeV, 0_GeV});
+    Vector<hepmomentum_d> P1(rootCS, {10_GeV, 10_GeV, 20_GeV});
+    Vector<hepmomentum_d> P2(rootCS, {0_GeV, 20_GeV, 0_GeV});
+
+    FourVector p0(E0, P0);
+    FourVector p1(E0, P1);
+    FourVector p2(E0, P2);
+
+    CHECK(p0.IsSpacelike());
+    CHECK(!p0.IsTimelike());
+    CHECK(!p0.IsPhotonlike());
+
+    CHECK(!p1.IsSpacelike());
+    CHECK(p1.IsTimelike());
+    CHECK(!p1.IsPhotonlike());
+
+    CHECK(!p2.IsSpacelike());
+    CHECK(!p2.IsTimelike());
+    CHECK(p2.IsPhotonlike());
+  }
+
+  /*
+    Test: P2 = E2/c2 - p2 with E in [GeV/c] and P in [GeV]
+    This requires additional factors of c
+   */
+  SECTION("Energy momentum in SI-units") {
+
+    auto E1 = 100_GeV / corsika::units::constants::c;
+    Vector<hepmomentum_d> P1(rootCS, {10_GeV, 5_GeV, 15_GeV});
+
+    FourVector p1(E1, P1);
+
+    const double check = 100 * 100 - 10 * 10 - 5 * 5 - 15 * 15; // for dummies...
+
+    REQUIRE(p1.GetNormSqr() / 1_GeV / 1_GeV == Approx(check));
+    REQUIRE(p1.GetNorm() / 1_GeV == Approx(sqrt(check)));
+  }
+
+  /**
+    Test: P2 = T2/c2 - r2 with T in [s] and r in [m]
+    This requires additional factors of c
+   */
+  SECTION("Spacetime in SI-units") {
+
+    TimeType T2 = 10_m / corsika::units::constants::c;
+    Vector<length_d> P2(rootCS, {10_m, 5_m, 5_m});
+
+    const double check = 10 * 10 - 10 * 10 - 5 * 5 - 5 * 5; // for dummies...
+
+    FourVector p2(T2, P2);
+
+    REQUIRE(p2.GetNormSqr() == check * 1_m * 1_m);
+    REQUIRE(p2.GetNorm() == sqrt(abs(check)) * 1_m);
+  }
+
+  /**
+     Testing the math operators
+   */
+  
+  SECTION("Operators and comutions") {
+
+    HEPEnergyType E1 = 100_GeV;
+    Vector<hepmomentum_d> P1(rootCS, {0_GeV, 0_GeV, 0_GeV});
+
+    HEPEnergyType E2 = 0_GeV;
+    Vector<hepmomentum_d> P2(rootCS, {10_GeV, 0_GeV, 0_GeV});
+
+    FourVector p1(E1, P1);
+    const FourVector p2(E2, P2);
+
+    REQUIRE(p1.GetNorm() / 1_GeV == Approx(100.));
+    REQUIRE(p2.GetNorm() / 1_GeV == Approx(10.));
+
+    SECTION("product") {
+      FourVector p3 = p1 + p2;
+      REQUIRE(p3.GetNorm() / 1_GeV == Approx(sqrt(100. * 100. - 100.)));
+      p3 -= p2;
+      REQUIRE(p3.GetNorm() / 1_GeV == Approx(100.));
+      REQUIRE(p1.GetNorm() / 1_GeV == Approx(100.));
+      REQUIRE(p2.GetNorm() / 1_GeV == Approx(10.));
+    }
+
+    SECTION("difference") {
+      FourVector p3 = p1 - p2;
+      REQUIRE(p3.GetNorm() / 1_GeV == Approx(sqrt(100. * 100. - 100.)));
+      p3 += p2;
+      REQUIRE(p3.GetNorm() / 1_GeV == Approx(100.));
+      REQUIRE(p1.GetNorm() / 1_GeV == Approx(100.));
+      REQUIRE(p2.GetNorm() / 1_GeV == Approx(10.));
+    }
+
+    SECTION("scale") {
+      double s = 10;
+      FourVector p3 = p1*s;
+      REQUIRE(p3.GetNorm() / 1_GeV == Approx(sqrt(100. * 100. * s * s)));
+      p3 /= 10;
+      REQUIRE(p3.GetNorm() / 1_GeV == Approx(sqrt(100. * 100. )));
+      REQUIRE(p1.GetNorm() / 1_GeV == Approx(100.));
+      REQUIRE(p2.GetNorm() / 1_GeV == Approx(10.));
+    }
+  }
+
+  /*
+  SECTION("Use as wrapper") {
+
+    TimeType T1 = 10_m / corsika::units::constants::c;
+    Vector<length_d> P1(rootCS, {10_m, 5_m, 5_m});
+
+    const TimeType T2 = 10_m / corsika::units::constants::c;
+    const Vector<length_d> P2(rootCS, {10_m, 5_m, 5_m});
+
+    FourVector p1(T1, P1);
+    FourVector p2(T2, P2);
+    FourVector p3(TimeType(10_m/corsika::units::constants::c), Vector<length_d>(rootCS,
+  {10_m,10_m,10_m}));
+
+    std::cout << type_id_with_cvr<decltype(p1)>().pretty_name() << std::endl;
+    std::cout << type_id_with_cvr<decltype(p2)>().pretty_name() << std::endl;
+    std::cout << type_id_with_cvr<decltype(p3)>().pretty_name() << std::endl;
+
+    const double check = 10 * 10 - 10 * 10 - 5 * 5 - 5 * 5; // for dummies...
+    REQUIRE(p1.GetNormSqr() == check * 1_m * 1_m);
+    REQUIRE(p2.GetNormSqr() == check * 1_m * 1_m);
+    REQUIRE(p3.GetNormSqr() == check * 1_m * 1_m);
+
+    REQUIRE(p1.GetNorm() == sqrt(abs(check)) * 1_m);
+    REQUIRE(p2.GetNorm() == sqrt(abs(check)) * 1_m);
+    REQUIRE(p3.GetNorm() == sqrt(abs(check)) * 1_m);
+  }
+  */
+}