From 47faeb67daf16b87c680546dc7d6c99341209fd0 Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Mon, 14 Jan 2019 13:51:02 +0100
Subject: [PATCH] added FourVector with tests
---
Framework/Geometry/FourVector.h | 223 +++++++++++++++++++++++++++
Framework/Geometry/testFourVector.cc | 187 ++++++++++++++++++++++
2 files changed, 410 insertions(+)
create mode 100644 Framework/Geometry/FourVector.h
create mode 100644 Framework/Geometry/testFourVector.cc
diff --git a/Framework/Geometry/FourVector.h b/Framework/Geometry/FourVector.h
new file mode 100644
index 000000000..76c5c7ecb
--- /dev/null
+++ b/Framework/Geometry/FourVector.h
@@ -0,0 +1,223 @@
+#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
diff --git a/Framework/Geometry/testFourVector.cc b/Framework/Geometry/testFourVector.cc
new file mode 100644
index 000000000..8bc17488d
--- /dev/null
+++ b/Framework/Geometry/testFourVector.cc
@@ -0,0 +1,187 @@
+
+/**
+ * (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);
+ }
+ */
+}
--
GitLab