From 604198b86c67e5d6891181e739179116a7fa0b2f Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Mon, 14 Jan 2019 14:48:41 +0100
Subject: [PATCH] also allow reference members in FourVector
---
Framework/Geometry/FourVector.h | 363 +++++++++++++--------------
Framework/Geometry/testFourVector.cc | 49 ++--
Framework/Units/PhysicalUnits.h | 12 +-
3 files changed, 206 insertions(+), 218 deletions(-)
diff --git a/Framework/Geometry/FourVector.h b/Framework/Geometry/FourVector.h
index 76c5c7ec..8d2b0945 100644
--- a/Framework/Geometry/FourVector.h
+++ b/Framework/Geometry/FourVector.h
@@ -9,215 +9,198 @@
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; }
+ /**
+ FourVector supports "full" units, e.g. E in [GeV/c] and p in [GeV],
+ or also t in [s] and r in [m], etc.
- auto GetNormSqr() const {
- return GetTimeSquared() - fSpaceLike.squaredNorm();
- }
+ However, for HEP applications it is also possible to use E and p
+ both in [GeV].
- SpaceType GetNorm() const { return sqrt(abs(GetNormSqr())); }
+ The FourVector can return NormSqr and Norm, whereas Norm is
+ sqrt(abs(NormSqr)). The physical units are always calculated and
+ returned properly.
- bool IsTimelike() const {
- return GetTimeSquared() < fSpaceLike.squaredNorm();
- } // Norm2 < 0
+ FourVector can also return if it is TimeLike, SpaceLike or PhotonLike.
- bool IsSpacelike() const {
- return GetTimeSquared() > fSpaceLike.squaredNorm();
- } // Norm2 > 0
+ 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.
+ */
- bool IsPhotonlike() const {
- return GetTimeSquared() == fSpaceLike.squaredNorm();
- } // // Norm2 == 0
+ template <typename TimeType, typename SpaceVecType>
+ class FourVector {
+
+ public:
+ using SpaceType = typename std::decay<SpaceVecType>::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:
+
+ FourVector(const TimeType& eT, const SpaceVecType& eS)
+ : fTimeLike(eT)
+ , fSpaceLike(eS) {
+ }
+
+ 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;
+ SpaceVecType fSpaceLike;
+
+ /// the friends: math operators
+ template <typename T, typename U>
+ friend FourVector<typename std::decay<T>::type, typename std::decay<U>::type> operator+(const FourVector<T, U>&, const FourVector<T, U>&);
- FourVector &operator+=(const FourVector &b) {
- fTimeLike += b.fTimeLike;
- fSpaceLike += b.fSpaceLike;
- return *this;
- }
+ template <typename T, typename U>
+ friend FourVector<typename std::decay<T>::type, typename std::decay<U>::type> operator-(const FourVector<T, U>&, const FourVector<T, U>&);
- FourVector &operator-=(const FourVector &b) {
- fTimeLike -= b.fTimeLike;
- fSpaceLike -= b.fSpaceLike;
- return *this;
- }
+ template <typename T, typename U>
+ friend FourVector<typename std::decay<T>::type, typename std::decay<U>::type> operator*(const FourVector<T, U>&, const double);
- FourVector &operator*=(const double b) {
- fTimeLike *= b;
- fSpaceLike *= b;
- return *this;
- }
+ template <typename T, typename U>
+ friend FourVector<typename std::decay<T>::type, typename std::decay<U>::type> operator/(const FourVector<T, U>&, const double);
+ };
- FourVector &operator/=(const double b) {
- fTimeLike /= b;
- fSpaceLike.GetComponents() /= b; // TODO: WHY IS THIS??????
- return *this;
+ /**
+ The math operator+
+ */
+ template <typename TimeType, typename SpaceVecType>
+ inline FourVector<typename std::decay<TimeType>::type,
+ typename std::decay<SpaceVecType>::type>
+ operator+(const FourVector<TimeType, SpaceVecType>& a,
+ const FourVector<TimeType, SpaceVecType>& b) {
+ return FourVector<typename std::decay<TimeType>::type,
+ typename std::decay<SpaceVecType>::type>(
+ a.fTimeLike + b.fTimeLike, a.fSpaceLike + b.fSpaceLike);
}
-
- FourVector &operator/(const double b) {
- *this /= b;
- return *this;
+
+ /**
+ The math operator-
+ */
+ template <typename TimeType, typename SpaceVecType>
+ inline FourVector<typename std::decay<TimeType>::type,
+ typename std::decay<SpaceVecType>::type>
+ operator-(const FourVector<TimeType, SpaceVecType>& a,
+ const FourVector<TimeType, SpaceVecType>& b) {
+ return FourVector<typename std::decay<TimeType>::type,
+ typename std::decay<SpaceVecType>::type>(
+ a.fTimeLike - b.fTimeLike, a.fSpaceLike - b.fSpaceLike);
}
/**
- 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();
+ The math operator*
+ */
+ template <typename TimeType, typename SpaceVecType>
+ inline FourVector<typename std::decay<TimeType>::type,
+ typename std::decay<SpaceVecType>::type>
+ operator*(const FourVector<TimeType, SpaceVecType>& a,
+ const double b) {
+ return FourVector<typename std::decay<TimeType>::type,
+ typename std::decay<SpaceVecType>::type>(
+ a.fTimeLike * b, a.fSpaceLike * b);
}
-private:
/**
- This function is automatically compiled to use of ignore the
- extra factor of "c" for the time-like quantity
+ The math operator/
*/
- 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;
+ template <typename TimeType, typename SpaceVecType>
+ inline FourVector<typename std::decay<TimeType>::type,
+ typename std::decay<SpaceVecType>::type>
+ operator/(const FourVector<TimeType, SpaceVecType>& a,
+ const double b) {
+ return FourVector<typename std::decay<TimeType>::type,
+ typename std::decay<SpaceVecType>::type>(
+ a.fTimeLike / b, a.fSpaceLike / b);
}
-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
index 8bc17488..b0188307 100644
--- a/Framework/Geometry/testFourVector.cc
+++ b/Framework/Geometry/testFourVector.cc
@@ -26,7 +26,6 @@ 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
@@ -112,7 +111,7 @@ TEST_CASE("four vectors") {
/**
Testing the math operators
*/
-
+
SECTION("Operators and comutions") {
HEPEnergyType E1 = 100_GeV;
@@ -147,41 +146,47 @@ TEST_CASE("four vectors") {
SECTION("scale") {
double s = 10;
- FourVector p3 = p1*s;
+ 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(p3.GetNorm() / 1_GeV == Approx(sqrt(100. * 100.)));
REQUIRE(p1.GetNorm() / 1_GeV == Approx(100.));
REQUIRE(p2.GetNorm() / 1_GeV == Approx(10.));
}
}
- /*
+ /**
+ The FourVector class can be used with reference template
+ arguments. In this configuration it does not hold any data
+ itself, but rather just refers to data located elsewhere. Thus,
+ it merely provides the physical/mathematical wrapper around the
+ data.
+ */
+
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});
+ TimeType T = 10_m / corsika::units::constants::c;
+ Vector<length_d> P(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}));
+ const TimeType T_c = 10_m / corsika::units::constants::c;
+ const Vector<length_d> P_c(rootCS, {10_m, 5_m, 5_m});
+ //FourVector<TimeType&, Vector<length_d>&> p0(T_c, P_c); // this does not compile, and it shoudn't!
+ FourVector<TimeType&, Vector<length_d>&> p1(T, P);
+ FourVector<const TimeType&, const Vector<length_d>&> p2(T, P);
+ FourVector<const TimeType&, const Vector<length_d>&> p3(T_c, P_c);
+
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;
+ p1 *= 10;
+ //p2 *= 10; // this does not compile, and it shoudn't !
+ //p3 *= 10; // this does not compile, and it shoudn't !!
+
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);
+ REQUIRE(p1.GetNormSqr()/(1_m * 1_m) == Approx(10.*10. * check));
+ REQUIRE(p2.GetNorm()/1_m == Approx(10*sqrt(abs(check))));
+ REQUIRE(p3.GetNorm()/1_m == Approx(sqrt(abs(check))));
}
- */
}
diff --git a/Framework/Units/PhysicalUnits.h b/Framework/Units/PhysicalUnits.h
index 5861f9a4..450dd9d0 100644
--- a/Framework/Units/PhysicalUnits.h
+++ b/Framework/Units/PhysicalUnits.h
@@ -102,12 +102,12 @@ namespace phys {
QUANTITY_DEFINE_SCALING_LITERALS(gram, mass_d, 1e-3)
- /*
- QUANTITY_DEFINE_LITERALS(meter, length_d)
- QUANTITY_DEFINE_LITERALS(second, time_interval_d)
- QUANTITY_DEFINE_LITERALS(ampere, electric_current_d)
- QUANTITY_DEFINE_LITERALS(Kelvin, thermodynamic_temperature_d)
- */
+ /*
+ QUANTITY_DEFINE_LITERALS(meter, length_d)
+ QUANTITY_DEFINE_LITERALS(second, time_interval_d)
+ QUANTITY_DEFINE_LITERALS(ampere, electric_current_d)
+ QUANTITY_DEFINE_LITERALS(Kelvin, thermodynamic_temperature_d)
+ */
} // namespace literals
} // namespace units
--
GitLab