diff --git a/Framework/Geometry/FourVector.h b/Framework/Geometry/FourVector.h
index 76c5c7ecb2386f4926042b6a53c4caccc71de638..8d2b0945270c24d6acbcb90bd992121b3da87e08 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 8bc17488d2c347520c272731133bfa28df099c6e..b018830796ab1b7b33302dd7d8b4e97a00725599 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/ProcessSequence/ProcessSequence.h b/Framework/ProcessSequence/ProcessSequence.h
index c6a0cb5444c13bd798a56b7c92c3c86eef9df087..48eafee2839b1035bcff8019e9b94b366c4f3f21 100644
--- a/Framework/ProcessSequence/ProcessSequence.h
+++ b/Framework/ProcessSequence/ProcessSequence.h
@@ -138,7 +138,8 @@ namespace corsika::process {
template <typename Particle, typename Stack>
EProcessReturn SelectInteraction(
- Particle& p, Stack& s, [[maybe_unused]]corsika::units::si::InverseGrammageType lambda_select,
+ Particle& p, Stack& s,
+ [[maybe_unused]] corsika::units::si::InverseGrammageType lambda_select,
corsika::units::si::InverseGrammageType& lambda_inv_count) {
if constexpr (is_process_sequence<T1type>::value) {
@@ -204,9 +205,10 @@ namespace corsika::process {
// select decay process
template <typename Particle, typename Stack>
- EProcessReturn SelectDecay(Particle& p, Stack& s,
- [[maybe_unused]] corsika::units::si::InverseTimeType decay_select,
- corsika::units::si::InverseTimeType& decay_inv_count) {
+ EProcessReturn SelectDecay(
+ Particle& p, Stack& s,
+ [[maybe_unused]] corsika::units::si::InverseTimeType decay_select,
+ corsika::units::si::InverseTimeType& decay_inv_count) {
if constexpr (is_process_sequence<T1>::value) {
// if A is a process sequence --> check inside
const EProcessReturn ret = A.SelectDecay(p, s, decay_select, decay_inv_count);
diff --git a/Framework/Units/PhysicalUnits.h b/Framework/Units/PhysicalUnits.h
index 5861f9a45d24aae10bdc6640244e9f0692d5cd18..450dd9d02b75642de13d4124739570daaa0d7749 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
diff --git a/Processes/Sibyll/Decay.h b/Processes/Sibyll/Decay.h
index a798f318db53e994ba7a1edff828e9010dec979d..e7ed3165dc46bd50e193298f8a63956becf028f4 100644
--- a/Processes/Sibyll/Decay.h
+++ b/Processes/Sibyll/Decay.h
@@ -169,7 +169,7 @@ namespace corsika::process {
using corsika::geometry::Point;
using namespace corsika::units::si;
- // TODO: this should be done in a central, common place. Not here..
+ // TODO: this should be done in a central, common place. Not here..
#ifndef CORSIKA_OSX
feenableexcept(FE_INVALID);
#endif
@@ -185,7 +185,8 @@ namespace corsika::process {
pin.SetMomentum(p.GetMomentum());
// setting particle mass with Corsika values, may be inconsistent with sibyll
// internal values
- // TODO: #warning setting particle mass with Corsika values, may be inconsistent with sibyll internal values
+ // TODO: #warning setting particle mass with Corsika values, may be inconsistent
+ // with sibyll internal values
pin.SetMass(corsika::particles::GetMass(pCode));
// remember position
Point decayPoint = p.GetPosition();
@@ -219,7 +220,7 @@ namespace corsika::process {
// empty sibyll stack
ss.Clear();
- // TODO: this should be done in a central, common place. Not here..
+ // TODO: this should be done in a central, common place. Not here..
#ifndef CORSIKA_OSX
fedisableexcept(FE_INVALID);
#endif
diff --git a/Processes/Sibyll/Interaction.h b/Processes/Sibyll/Interaction.h
index fac0dd578b77cfea776b1b942b6a42b9e899f809..4375bb2f0d21fe71784e60a554d008cc6ab5b9df 100644
--- a/Processes/Sibyll/Interaction.h
+++ b/Processes/Sibyll/Interaction.h
@@ -14,13 +14,13 @@
#include <corsika/process/InteractionProcess.h>
+#include <corsika/particles/ParticleProperties.h>
#include <corsika/process/sibyll/ParticleConversion.h>
#include <corsika/process/sibyll/SibStack.h>
#include <corsika/process/sibyll/sibyll2.3c.h>
-#include <corsika/utl/COMBoost.h>
-#include <corsika/particles/ParticleProperties.h>
#include <corsika/random/RNGManager.h>
#include <corsika/units/PhysicalUnits.h>
+#include <corsika/utl/COMBoost.h>
namespace corsika::process::sibyll {
@@ -186,7 +186,7 @@ namespace corsika::process::sibyll {
// FOR NOW: target is always at rest
auto constexpr nucleon_mass = 0.5 * (corsika::particles::Proton::GetMass() +
- corsika::particles::Neutron::GetMass());
+ corsika::particles::Neutron::GetMass());
const auto Etarget = 0_GeV + nucleon_mass;
const auto pTarget = MomentumVector(rootCS, 0_GeV, 0_GeV, 0_GeV);
cout << "target momentum (GeV/c): " << pTarget.GetComponents() / 1_GeV << endl;
@@ -226,26 +226,27 @@ namespace corsika::process::sibyll {
std::cout << "Interaction: "
<< " DoDiscrete: gambet:" << gambet.GetComponents() << endl;
- auto const pProjectileLab = p.GetMomentum();
- //{rootCS, {0_GeV / c, 1_PeV / c, 0_GeV / c}};
- HEPEnergyType const eProjectileLab = p.GetEnergy();
- //energy(projectileMass, pProjectileLab);
+ auto const pProjectileLab = p.GetMomentum();
+ //{rootCS, {0_GeV / c, 1_PeV / c, 0_GeV / c}};
+ HEPEnergyType const eProjectileLab = p.GetEnergy();
+ // energy(projectileMass, pProjectileLab);
+
+ // define target kinematics in lab frame
+ HEPMassType const targetMass = nucleon_mass;
+ // define boost to com frame
+ COMBoost const boost(eProjectileLab, pProjectileLab, targetMass);
- // define target kinematics in lab frame
- HEPMassType const targetMass = nucleon_mass;
- // define boost to com frame
- COMBoost const boost(eProjectileLab, pProjectileLab, targetMass);
+ cout << "Interaction: new boost: ebeam lab: " << eProjectileLab / 1_GeV << endl
+ << "Interaction: new boost: pbeam lab: "
+ << pProjectileLab.GetComponents() / 1_GeV << endl;
- cout << "Interaction: new boost: ebeam lab: " << eProjectileLab / 1_GeV << endl
- << "Interaction: new boost: pbeam lab: " << pProjectileLab.GetComponents() / 1_GeV << endl;
+ // boost projecticle
+ auto const [eProjectileCoM, pProjectileCoM] =
+ boost.toCoM(eProjectileLab, pProjectileLab);
- // boost projecticle
- auto const [eProjectileCoM, pProjectileCoM] =
- boost.toCoM(eProjectileLab, pProjectileLab);
+ cout << "Interaction: new boost: ebeam com: " << eProjectileCoM / 1_GeV << endl
+ << "Interaction: new boost: pbeam com: " << pProjectileCoM / 1_GeV << endl;
- cout << "Interaction: new boost: ebeam com: " << eProjectileCoM / 1_GeV << endl
- << "Interaction: new boost: pbeam com: " << pProjectileCoM / 1_GeV << endl;
-
int kBeam = process::sibyll::ConvertToSibyllRaw(p.GetPID());
std::cout << "Interaction: "
diff --git a/Processes/Sibyll/testSibyll.cc b/Processes/Sibyll/testSibyll.cc
index 3d7d425d9c63aaacf5ffcef5723dfac730283b7e..844d4c3ac3e06b86ca0fe21ae11803d86c32b6a1 100644
--- a/Processes/Sibyll/testSibyll.cc
+++ b/Processes/Sibyll/testSibyll.cc
@@ -70,9 +70,9 @@ TEST_CASE("Sibyll", "[processes]") {
#include <corsika/units/PhysicalUnits.h>
+#include <corsika/particles/ParticleProperties.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
-#include <corsika/particles/ParticleProperties.h>
using namespace corsika::units::si;
using namespace corsika::units;
@@ -91,9 +91,9 @@ TEST_CASE("SibyllInterface", "[processes]") {
setup::Stack stack;
auto particle = stack.NewParticle();
-
+
Interaction model;
-
+
model.Init();
[[maybe_unused]] const process::EProcessReturn ret =
model.DoInteraction(particle, stack);
@@ -108,7 +108,8 @@ TEST_CASE("SibyllInterface", "[processes]") {
{
const HEPEnergyType E0 = 10_GeV;
particle.SetPID(particles::Code::Proton);
- HEPMomentumType P0 = sqrt(E0 * E0 - particles::Proton::GetMass() * particles::Proton::GetMass());
+ HEPMomentumType P0 =
+ sqrt(E0 * E0 - particles::Proton::GetMass() * particles::Proton::GetMass());
auto plab = stack::super_stupid::MomentumVector(cs, {0_GeV, 0_GeV, -P0});
particle.SetEnergy(E0);
particle.SetMomentum(plab);
@@ -116,9 +117,9 @@ TEST_CASE("SibyllInterface", "[processes]") {
geometry::Point p(cs, 0_m, 0_m, 0_m);
particle.SetPosition(p);
}
-
+
Decay model;
-
+
model.Init();
/*[[maybe_unused]] const process::EProcessReturn ret =*/model.DoDecay(particle,
stack);