IAP GITLAB

Skip to content
Snippets Groups Projects
Commit 604198b8 authored by ralfulrich's avatar ralfulrich
Browse files

also allow reference members in FourVector

parent 8bf34d9f
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Framework/Geometry/FourVector.h
+ 173
− 190
View file @ 604198b8
......@@ -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
Framework/Geometry/testFourVector.cc
+ 27
− 22
View file @ 604198b8
......@@ -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))));
}
*/
}
Framework/Units/PhysicalUnits.h
+ 6
− 6
View file @ 604198b8
......@@ -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
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

data privacy - impressum