diff --git a/corsika/detail/framework/geometry/FourVector.inl b/corsika/detail/framework/geometry/FourVector.inl index 05fee3c86d1767363af3c72ba1ff26ed431f6ee4..f3de79d9e81518f034ead3b309c4f7186997a9e3 100644 --- a/corsika/detail/framework/geometry/FourVector.inl +++ b/corsika/detail/framework/geometry/FourVector.inl @@ -10,159 +10,109 @@ #pragma once -#include <iostream> #include <type_traits> #include <corsika/framework/core/PhysicalUnits.hpp> #include <corsika/framework/geometry/Vector.hpp> namespace corsika { - template <typename TimeType, typename SpaceVecType> - TimeType FourVector<TimeType, SpaceVecType>::GetTimeLikeComponent() const { - return fTimeLike; + template <typename TTimeType, typename TSpaceVecType> + TTimeType FourVector<TTimeType, TSpaceVecType>::getTimeLikeComponent() const { + return timeLike_; } - template <typename TimeType, typename SpaceVecType> - SpaceVecType& FourVector<TimeType, SpaceVecType>::GetSpaceLikeComponents() { - return fSpaceLike; + template <typename TTimeType, typename TSpaceVecType> + TSpaceVecType& FourVector<TTimeType, TSpaceVecType>::spaceLikeComponents() { + return spaceLike_; } - template <typename TimeType, typename SpaceVecType> - const SpaceVecType& FourVector<TimeType, SpaceVecType>::GetSpaceLikeComponents() const { - return fSpaceLike; + template <typename TTimeType, typename TSpaceVecType> + const TSpaceVecType& FourVector<TTimeType, TSpaceVecType>::getSpaceLikeComponents() + const { + return spaceLike_; } - template <typename TimeType, typename SpaceVecType> - auto FourVector<TimeType, SpaceVecType>::GetNormSqr() const { - return GetTimeSquared() - fSpaceLike.squaredNorm(); + template <typename TTimeType, typename TSpaceVecType> + typename FourVector<TTimeType, TSpaceVecType>::norm_square_type + FourVector<TTimeType, TSpaceVecType>::getNormSqr() const { + return getTimeSquared() - spaceLike_.squaredNorm(); } - template <typename TimeType, typename SpaceVecType> - typename FourVector<TimeType, SpaceVecType>::SpaceType - FourVector<TimeType, SpaceVecType>::GetNorm() const { + template <typename TTimeType, typename TSpaceVecType> + typename FourVector<TTimeType, TSpaceVecType>::norm_type + FourVector<TTimeType, TSpaceVecType>::getNorm() const { - return sqrt(abs(GetNormSqr())); + return sqrt(abs(getNormSqr())); } - template <typename TimeType, typename SpaceVecType> - bool FourVector<TimeType, SpaceVecType>::IsTimelike() const { - return GetTimeSquared() < fSpaceLike.squaredNorm(); + template <typename TTimeType, typename TSpaceVecType> + bool FourVector<TTimeType, TSpaceVecType>::isTimelike() const { + return getTimeSquared() < spaceLike_.squaredNorm(); } - template <typename TimeType, typename SpaceVecType> - bool FourVector<TimeType, SpaceVecType>::IsSpacelike() const { - return GetTimeSquared() > fSpaceLike.squaredNorm(); + template <typename TTimeType, typename TSpaceVecType> + bool FourVector<TTimeType, TSpaceVecType>::isSpacelike() const { + return getTimeSquared() > spaceLike_.squaredNorm(); } - template <typename TimeType, typename SpaceVecType> - FourVector<TimeType, SpaceVecType>& FourVector<TimeType, SpaceVecType>::operator+=( + template <typename TTimeType, typename TSpaceVecType> + FourVector<TTimeType, TSpaceVecType>& FourVector<TTimeType, TSpaceVecType>::operator+=( const FourVector& b) { - fTimeLike += b.fTimeLike; - fSpaceLike += b.fSpaceLike; + timeLike_ += b.timeLike_; + spaceLike_ += b.spaceLike_; return *this; } - template <typename TimeType, typename SpaceVecType> - FourVector<TimeType, SpaceVecType>& FourVector<TimeType, SpaceVecType>::operator-=( + template <typename TTimeType, typename TSpaceVecType> + FourVector<TTimeType, TSpaceVecType>& FourVector<TTimeType, TSpaceVecType>::operator-=( const FourVector& b) { - fTimeLike -= b.fTimeLike; - fSpaceLike -= b.fSpaceLike; + timeLike_ -= b.timeLike_; + spaceLike_ -= b.spaceLike_; return *this; } - template <typename TimeType, typename SpaceVecType> - FourVector<TimeType, SpaceVecType>& FourVector<TimeType, SpaceVecType>::operator*=( + template <typename TTimeType, typename TSpaceVecType> + FourVector<TTimeType, TSpaceVecType>& FourVector<TTimeType, TSpaceVecType>::operator*=( const double b) { - fTimeLike *= b; - fSpaceLike *= b; + timeLike_ *= b; + spaceLike_ *= b; return *this; } - template <typename TimeType, typename SpaceVecType> - FourVector<TimeType, SpaceVecType>& FourVector<TimeType, SpaceVecType>::operator/=( + template <typename TTimeType, typename TSpaceVecType> + FourVector<TTimeType, TSpaceVecType>& FourVector<TTimeType, TSpaceVecType>::operator/=( const double b) { - fTimeLike /= b; - fSpaceLike.GetComponents() /= b; // TODO: WHY IS THIS?????? + timeLike_ /= b; + spaceLike_.GetComponents() /= b; return *this; } - template <typename TimeType, typename SpaceVecType> - FourVector<TimeType, SpaceVecType>& FourVector<TimeType, SpaceVecType>::operator/( + template <typename TTimeType, typename TSpaceVecType> + FourVector<TTimeType, TSpaceVecType>& FourVector<TTimeType, TSpaceVecType>::operator/( const double b) { *this /= b; return *this; } - template <typename TimeType, typename SpaceVecType> - typename FourVector<TimeType, SpaceVecType>::SpaceType - FourVector<TimeType, SpaceVecType>::operator*(const FourVector& b) { - if constexpr (std::is_same<typename std::decay<TimeType>::type, - decltype(std::declval<SpaceType>() / meter * - second)>::value) - return fTimeLike * b.fTimeLike * (constants::c * constants::c) - fSpaceLike.norm(); + template <typename TTimeType, typename TSpaceVecType> + typename FourVector<TTimeType, TSpaceVecType>::norm_type + FourVector<TTimeType, TSpaceVecType>::operator*(const FourVector& b) { + if constexpr (std::is_same<time_type, decltype(std::declval<space_type>() / meter * + second)>::value) + return timeLike_ * b.timeLike_ * constants::cSquared - spaceLike_.norm(); else - return fTimeLike * fTimeLike - fSpaceLike.norm(); + return timeLike_ * timeLike_ - spaceLike_.norm(); } - template <typename TimeType, typename SpaceVecType> - auto FourVector<TimeType, SpaceVecType>::GetTimeSquared() const { - if constexpr (std::is_same<typename std::decay<TimeType>::type, - decltype(std::declval<SpaceType>() / meter * - second)>::value) - return fTimeLike * fTimeLike * constants::cSquared; + template <typename TTimeType, typename TSpaceVecType> + typename FourVector<TTimeType, TSpaceVecType>::norm_square_type + FourVector<TTimeType, TSpaceVecType>::getTimeSquared() const { + if constexpr (std::is_same<time_type, decltype(std::declval<space_type>() / meter * + second)>::value) + return timeLike_ * timeLike_ * constants::cSquared; else - return fTimeLike * fTimeLike; - } - - /** - 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); - } - - /** - 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); - } - - /** - 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); - } - - /** - 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); + return timeLike_ * timeLike_; } } // namespace corsika diff --git a/corsika/framework/geometry/FourVector.hpp b/corsika/framework/geometry/FourVector.hpp index 3223914803e1f4a4339f3439a11d733c658308e3..c4f8cf1ac0916a1cd4ce10c4c9d5fd6fbb068c48 100644 --- a/corsika/framework/geometry/FourVector.hpp +++ b/corsika/framework/geometry/FourVector.hpp @@ -16,185 +16,179 @@ namespace corsika { /** - FourVector supports "full" units, e.g. E in [GeV/c] and p in [GeV], + FourVector fully supports 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]. + Thus, the input units of time-like and space-like coordinates + must either be idential (e.g. GeV) or scaled by "c" as in + [E/c]=[p]. + + 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. + When a FourVector is initialized with a lvalue references, + e.g. as `FourVector<TimeType&, Vector<length_d>&>`, references + are also used as internal data types, which should lead to + complete disappearance of the FourVector class during + optimization. */ - template <typename TimeType, typename SpaceVecType> + template <typename TTimeType, typename TSpaceVecType> class FourVector { public: - using SpaceType = typename std::decay<SpaceVecType>::type::Quantity; + using space_vec_type = typename std::decay<TSpaceVecType>::type; + using space_type = typename space_vec_type::Quantity; + using time_type = typename std::decay<TTimeType>::type; //! 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>() / meter * second)>::value, - "Units of time-like and space-like coordinates must either be idential " - "(e.g. GeV) or [E/c]=[p]"); + static_assert(std::is_same<time_type, space_type>::value || + std::is_same<time_type, decltype(std::declval<space_type>() / + meter * second)>::value, + "Units of time-like and space-like coordinates must either be idential " + "(e.g. GeV) or [E/c]=[p]"); + + using norm_type = space_type; + using norm_square_type = + decltype(std::declval<norm_type>() * std::declval<norm_type>()); public: FourVector() = default; - FourVector(const TimeType& eT, const SpaceVecType& eS) - : fTimeLike(eT) - , fSpaceLike(eS) {} + FourVector(TTimeType const& eT, TSpaceVecType const& eS) + : timeLike_(eT) + , spaceLike_(eS) {} /* * FIXME: These Getters are mis-leading and does not favor * locality. Adhere to Getter/Setter */ /** - * @brief * - * @return fTimeLike + * @return timeLike_ */ - TimeType GetTimeLikeComponent() const; + TTimeType getTimeLikeComponent() const; /** - * @brief * - * @return fSpaceLike + * @return spaceLike_ */ - SpaceVecType& GetSpaceLikeComponents(); + TSpaceVecType& spaceLikeComponents(); /** - * @brief * - * @return fSpaceLike; + * @return spaceLike_; */ - const SpaceVecType& GetSpaceLikeComponents() const; + TSpaceVecType const& getSpaceLikeComponents() const; /** - * @brief * - * @return + * @return $p_0^2 - \vec{p}^2$ */ - auto GetNormSqr() const; + norm_square_type getNormSqr() const; /** - * @brief * - * @return + * @return $sqrt(p_0^2 - \vec{p}^2)$ */ - SpaceType GetNorm() const; + norm_type getNorm() const; /* * FIXME: a better alternative would be to define an enumeration * enum { SpaceLike =-1, TimeLike, LightLike } V4R_Category; * and a method called V4R_Category GetCategory() const; + * + * RU: then you have to decide in the constructor which avoids "lazyness" */ /** - * @brief * * @return */ - bool IsTimelike() const; + bool isTimelike() const; /** - * @brief - * - * @return + * @defgroup math operators (class members) + * @{ */ - bool IsSpacelike() const; - - FourVector& operator+=(const FourVector& b); - - FourVector& operator-=(const FourVector& b); - - FourVector& operator*=(const double b); - - FourVector& operator/=(const double b); - - FourVector& operator/(const double b); - + bool isSpacelike() const; + + FourVector& operator+=(FourVector const&); + FourVector& operator-=(const FourVector&); + FourVector& operator*=(const double); + FourVector& operator/=(const double); + FourVector& operator/(const double); + /** + Scalar product of two FourVectors + 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); - + norm_type operator*(const FourVector& b); + + /** @} */ + 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>&); + TTimeType timeLike_; + TSpaceVecType spaceLike_; - 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>&); - - template <typename T, typename U> - friend FourVector<typename std::decay<T>::type, typename std::decay<U>::type> - operator*(const FourVector<T, U>&, const double); - - template <typename T, typename U> - friend FourVector<typename std::decay<T>::type, typename std::decay<U>::type> - operator/(const FourVector<T, U>&, const double); + /** + * @defgroup the friends: (free) math operators + * @{ + * + * We need to define them inline here since we do not want to + * implement them as template functions. They are valid only for + * the specific types as defined right here. + * + * + * Note, these are "free function" (event if they don't look as + * such). Thus, even if the input object uses internal references + * for storage, the free math operators, of course, must provide + * value-copies. + * + **/ + friend FourVector<time_type, space_vec_type> operator+(FourVector const& a, + FourVector const& b) { + return FourVector<time_type, space_vec_type>(a.timeLike_ + b.timeLike_, + a.spaceLike_ + b.spaceLike_); + } + + friend FourVector<time_type, space_vec_type> operator-(FourVector const& a, + FourVector const& b) { + return FourVector<time_type, space_vec_type>(a.timeLike_ - b.timeLike_, + a.spaceLike_ - b.spaceLike_); + } + + friend FourVector<time_type, space_vec_type> operator*(FourVector const& a, + const double b) { + return FourVector<time_type, space_vec_type>(a.timeLike_ * b, a.spaceLike_ * b); + } + + friend FourVector<time_type, space_vec_type> operator/(FourVector const& a, + double b) { + return FourVector<time_type, space_vec_type>(a.timeLike_ / b, a.spaceLike_ / b); + } + + /** @} */ private: /** - This function is automatically compiled to use of ignore the - extra factor of "c" for the time-like quantity + This function is there to automatically remove the eventual + extra factor of "c" for the time-like quantity. */ - auto GetTimeSquared() const; + norm_square_type getTimeSquared() const; }; - /** - 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); - - /** - 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); - - /** - The math operator* - FIXME: Add overload to deal with multiplication by a scalar and 3-vectors - */ - 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); - - /** - 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); - } // namespace corsika #include <corsika/detail/framework/geometry/FourVector.inl> diff --git a/tests/framework/testFourVector.cpp b/tests/framework/testFourVector.cpp index 06a96f3a375beddb7737da40ea64128fd80afb3c..e81277c6063764be1ad86fec0b25d942a35ade24 100644 --- a/tests/framework/testFourVector.cpp +++ b/tests/framework/testFourVector.cpp @@ -38,8 +38,8 @@ TEST_CASE("four vectors") { FourVector p0(E0, P0); - REQUIRE(p0.GetNormSqr() == -200_GeV * 1_GeV); - REQUIRE(p0.GetNorm() == sqrt(200_GeV * 1_GeV)); + REQUIRE(p0.getNormSqr() == -200_GeV * 1_GeV); + REQUIRE(p0.getNorm() == sqrt(200_GeV * 1_GeV)); } /* @@ -56,17 +56,14 @@ TEST_CASE("four vectors") { FourVector p1(E0, P1); FourVector p2(E0, P2); - CHECK(p0.IsSpacelike()); - CHECK(!p0.IsTimelike()); - // CHECK(!p0.IsPhotonlike()); + CHECK(p0.isSpacelike()); + CHECK(!p0.isTimelike()); - CHECK(!p1.IsSpacelike()); - CHECK(p1.IsTimelike()); - // CHECK(!p1.IsPhotonlike()); + CHECK(!p1.isSpacelike()); + CHECK(p1.isTimelike()); - CHECK(!p2.IsSpacelike()); - CHECK(!p2.IsTimelike()); - // CHECK(p2.IsPhotonlike()); + CHECK(!p2.isSpacelike()); + CHECK(!p2.isTimelike()); } /* @@ -82,8 +79,8 @@ TEST_CASE("four vectors") { 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))); + REQUIRE(p1.getNormSqr() / 1_GeV / 1_GeV == Approx(check)); + REQUIRE(p1.getNorm() / 1_GeV == Approx(sqrt(check))); } /** @@ -99,8 +96,8 @@ TEST_CASE("four vectors") { FourVector p2(T2, P2); - REQUIRE(p2.GetNormSqr() == check * 1_m * 1_m); - REQUIRE(p2.GetNorm() == sqrt(abs(check)) * 1_m); + REQUIRE(p2.getNormSqr() == check * 1_m * 1_m); + REQUIRE(p2.getNorm() == sqrt(abs(check)) * 1_m); } /** @@ -118,35 +115,35 @@ TEST_CASE("four vectors") { FourVector p1(E1, P1); const FourVector p2(E2, P2); - REQUIRE(p1.GetNorm() / 1_GeV == Approx(100.)); - REQUIRE(p2.GetNorm() / 1_GeV == Approx(10.)); + 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.))); + 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.)); + 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.))); + 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.)); + 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))); + 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.)); + REQUIRE(p3.getNorm() / 1_GeV == Approx(sqrt(100. * 100.))); + REQUIRE(p1.getNorm() / 1_GeV == Approx(100.)); + REQUIRE(p2.getNorm() / 1_GeV == Approx(10.)); } } @@ -177,8 +174,8 @@ TEST_CASE("four vectors") { // 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() / (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)))); + 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)))); } }