diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 09936820642852b0f9b0949dd6abb40faa94a3bc..50a951202666c2b29a2a9c863f96a2841b9f6356 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -18,7 +18,9 @@ build:
- cmake ..
- cmake --build .
- ctest -j4 -V >& test.log
- - gzip -9 test.log
+ - gzip -9 -S .gz test.log
+ - ls
+ - pwd
artifacts:
paths:
- build/test.log.gz
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index 8e21990c4557a9b4bdcfecba6f746e0fe0497774..9e44b146b0958fbfddf2b297b82121deab04f47a 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -230,7 +230,7 @@ int main() {
// setup processes, decays and interactions
tracking_line::TrackingLine<setup::Stack> tracking(env);
stack_inspector::StackInspector<setup::Stack> p0(true);
-
+
corsika::random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
corsika::process::sibyll::Interaction sibyll(env);
corsika::process::sibyll::Decay decay;
@@ -246,7 +246,7 @@ int main() {
// setup particle stack, and add primary particle
setup::Stack stack;
stack.Clear();
- const HEPEnergyType E0 = 100_TeV;
+ const HEPEnergyType E0 = 100_TeV; // 1_PeV crashes with bad COMboost in second interaction (crash later)
double theta = 0.;
double phi = 0.;
{
diff --git a/Environment/Environment.h b/Environment/Environment.h
index 1ab22e03877300d00c3270d17b01a4a378671db5..1fc1e443b26c617acb56b5d4b9406dacaed87514 100644
--- a/Environment/Environment.h
+++ b/Environment/Environment.h
@@ -21,7 +21,7 @@
namespace corsika::environment {
using BaseNodeType = VolumeTreeNode<corsika::setup::IEnvironmentModel>;
-
+
struct Universe : public corsika::geometry::Sphere {
Universe(corsika::geometry::CoordinateSystem const& pCS)
: corsika::geometry::Sphere(
diff --git a/Environment/IMediumModel.h b/Environment/IMediumModel.h
index 4278a79a3afc6a5b22be9c1dbc787add9126fe15..cdfb30cd65a5e8d2a1cfc7f4157d9a88d564cb3b 100644
--- a/Environment/IMediumModel.h
+++ b/Environment/IMediumModel.h
@@ -16,8 +16,8 @@
#include <corsika/geometry/Line.h>
#include <corsika/geometry/Point.h>
#include <corsika/geometry/Trajectory.h>
-#include <random>
#include <corsika/units/PhysicalUnits.h>
+#include <random>
namespace corsika::environment {
@@ -39,18 +39,19 @@ namespace corsika::environment {
corsika::units::si::GrammageType) const = 0;
virtual NuclearComposition const& GetNuclearComposition() const = 0;
-
+
template <class TRNG>
corsika::particles::Code SampleTarget(
- std::vector<corsika::units::si::CrossSectionType> const& sigma, TRNG& randomStream) const {
+ std::vector<corsika::units::si::CrossSectionType> const& sigma,
+ TRNG& randomStream) const {
using namespace corsika::units::si;
-
+
auto const& nuclComp = GetNuclearComposition();
auto const& fractions = nuclComp.GetFractions();
assert(sigma.size() == fractions.size());
-
+
std::vector<float> weights(fractions.size());
-
+
for (size_t i = 0; i < fractions.size(); ++i) {
std::cout << "HomogeneousMedium: fraction: " << fractions[i] << std::endl;
weights[i] = fractions[i] * sigma[i].magnitude();
diff --git a/Environment/NuclearComposition.h b/Environment/NuclearComposition.h
index 3bc1101a4246dbe46a53b0ddb1fffb3c912171cd..3678b1305bb6b0e2bb29dae4ac18917afd0fa5fc 100644
--- a/Environment/NuclearComposition.h
+++ b/Environment/NuclearComposition.h
@@ -13,9 +13,9 @@
#define _include_NuclearComposition_h
#include <corsika/particles/ParticleProperties.h>
+#include <cassert>
#include <numeric>
#include <stdexcept>
-#include <cassert>
#include <vector>
namespace corsika::environment {
@@ -29,7 +29,7 @@ namespace corsika::environment {
std::vector<float> pFractions)
: fNumberFractions(pFractions)
, fComponents(pComponents) {
-
+
assert(pComponents.size() == pFractions.size());
auto const sumFractions =
std::accumulate(pFractions.cbegin(), pFractions.cend(), 0.f);
diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h
index 999894f70efbaa015c5b7b82ed4cb234dc8f55fc..024044aea5d6699dd311352f24ec91459a180845 100644
--- a/Framework/Cascade/Cascade.h
+++ b/Framework/Cascade/Cascade.h
@@ -56,7 +56,7 @@ namespace corsika::cascade {
}
}
-private:
+ private:
void Step(Particle& particle) {
using namespace corsika::units::si;
diff --git a/Framework/Cascade/testCascade.cc b/Framework/Cascade/testCascade.cc
index 12885054b18cafdb1312818418f0bf46e95c9a88..33914fe527128826155a449452f0c5af0f8e8ac4 100644
--- a/Framework/Cascade/testCascade.cc
+++ b/Framework/Cascade/testCascade.cc
@@ -68,11 +68,14 @@ corsika::environment::Environment MakeDummyEnv() {
return env;
}
-static int fCount = 0;
-
class ProcessSplit : public corsika::process::ContinuousProcess<ProcessSplit> {
+
+ int fCount = 0;
+ int fCalls = 0;
+ HEPEnergyType fEcrit;
+
public:
- ProcessSplit() {}
+ ProcessSplit(HEPEnergyType e) : fEcrit(e) {}
template <typename Particle, typename T>
LengthType MaxStepLength(Particle&, T&) const {
@@ -80,9 +83,10 @@ public:
}
template <typename Particle, typename T, typename Stack>
- EProcessReturn DoContinuous(Particle& p, T&, Stack& s) const {
+ EProcessReturn DoContinuous(Particle& p, T&, Stack& s) {
+ fCalls++;
HEPEnergyType E = p.GetEnergy();
- if (E < 85_MeV) {
+ if (E < fEcrit) {
p.Delete();
fCount++;
} else {
@@ -98,9 +102,10 @@ public:
return EProcessReturn::eOk;
}
- void Init() { fCount = 0; }
+ void Init() { fCount = 0; fCalls =0; }
int GetCount() const { return fCount; }
+ int GetCalls() const { return fCalls; }
private:
};
@@ -113,7 +118,9 @@ TEST_CASE("Cascade", "[Cascade]") {
tracking_line::TrackingLine<setup::Stack> tracking(env);
stack_inspector::StackInspector<setup::Stack> p0(true);
- ProcessSplit p1;
+
+ const HEPEnergyType Ecrit = 85_MeV;
+ ProcessSplit p1(Ecrit);
auto sequence = p0 << p1;
setup::Stack stack;
@@ -133,6 +140,9 @@ TEST_CASE("Cascade", "[Cascade]") {
EAS.Init();
EAS.Run();
+ CHECK( p1.GetCount() == 2048 );
+ CHECK( p1.GetCalls() == 4095 );
+
/*
SECTION("sectionTwo") {
for (int i = 0; i < 0; ++i) {
diff --git a/Framework/Geometry/CMakeLists.txt b/Framework/Geometry/CMakeLists.txt
index 4d708ff109e6d259a352e82536cece0f0bc4c53a..bb96c2ea6bf5cf6e228b23fb7471c73c1b282923 100644
--- a/Framework/Geometry/CMakeLists.txt
+++ b/Framework/Geometry/CMakeLists.txt
@@ -17,7 +17,7 @@ set (
BaseVector.h
QuantityVector.h
Trajectory.h
- # BaseTrajectory.h
+ FourVector.h
)
set (
@@ -57,16 +57,22 @@ install (
PUBLIC_HEADER DESTINATION include/${GEOMETRY_NAMESPACE}
)
-
# --------------------
# code unit testing
add_executable (testGeometry testGeometry.cc)
-
target_link_libraries (
testGeometry
CORSIKAgeometry
CORSIKAunits
CORSIKAthirdparty # for catch2
)
-
CORSIKA_ADD_TEST(testGeometry)
+
+add_executable (testFourVector testFourVector.cc)
+target_link_libraries (
+ testFourVector
+ CORSIKAgeometry
+ CORSIKAunits
+ CORSIKAthirdparty # for catch2
+ )
+CORSIKA_ADD_TEST(testFourVector)
diff --git a/Framework/Geometry/FourVector.h b/Framework/Geometry/FourVector.h
new file mode 100644
index 0000000000000000000000000000000000000000..d3fe31f1838e64e052587f474baa2fdc259870bf
--- /dev/null
+++ b/Framework/Geometry/FourVector.h
@@ -0,0 +1,210 @@
+#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 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 GetTimeLikeComponent() const { return fTimeLike; }
+ SpaceVecType& GetSpaceLikeComponents() { return fSpaceLike; }
+ const SpaceVecType& GetSpaceLikeComponents() const { return fSpaceLike; }
+
+ 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
+
+ /* this is not numerically stable
+ 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>&);
+
+ 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);
+ };
+
+ /**
+ 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);
+ }
+
+} // namespace corsika::geometry
+
+#endif
diff --git a/Framework/Geometry/QuantityVector.h b/Framework/Geometry/QuantityVector.h
index bb498f086d7f0a20fbfd1aea0682f182c7dcaf80..e0893e52706312ee6ed8b58478ee4b6802ff6cf1 100644
--- a/Framework/Geometry/QuantityVector.h
+++ b/Framework/Geometry/QuantityVector.h
@@ -29,8 +29,7 @@ namespace corsika::geometry {
template <typename dim>
class QuantityVector {
- protected:
- // todo: check if we need to move "quantity" into namespace corsika::units
+ public:
using Quantity = phys::units::quantity<dim, double>; //< the phys::units::quantity
// corresponding to the dimension
diff --git a/Framework/Geometry/Vector.h b/Framework/Geometry/Vector.h
index 94394b90dd859761d4f15d0a2389f173726ee34d..712eece61380c75bd36dea77b17a5650690b79b5 100644
--- a/Framework/Geometry/Vector.h
+++ b/Framework/Geometry/Vector.h
@@ -31,6 +31,7 @@ namespace corsika::geometry {
template <typename dim>
class Vector : public BaseVector<dim> {
+ public:
using Quantity = phys::units::quantity<dim, double>;
public:
diff --git a/Framework/Geometry/testFourVector.cc b/Framework/Geometry/testFourVector.cc
new file mode 100644
index 0000000000000000000000000000000000000000..301fb97e93d71dccc79fe0eda844376b7ecc133d
--- /dev/null
+++ b/Framework/Geometry/testFourVector.cc
@@ -0,0 +1,193 @@
+
+/**
+ * (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.));
+ }
+ }
+
+ /**
+ 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 T = 10_m / corsika::units::constants::c;
+ Vector<length_d> P(rootCS, {10_m, 5_m, 5_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() / (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 6af26b0abdbc9cac073d1c833e085c4c8295ce78..48eafee2839b1035bcff8019e9b94b366c4f3f21 100644
--- a/Framework/ProcessSequence/ProcessSequence.h
+++ b/Framework/ProcessSequence/ProcessSequence.h
@@ -261,39 +261,6 @@ namespace corsika::process {
return ProcessSequence<P1, P2>(A.GetRef(), B.GetRef());
}
- /* #define OPSEQ(C1, C2) \ */
- /* template < \ */
- /* typename P1, typename P2, \ */
- /* typename std::enable_if<is_process<typename std::decay<P1>::type>::value &&
- * \ */
- /* is_process<typename
- * std::decay<P2>::type>::value>::type...> \ */
- /* inline auto operator+(P1&& A, P2&& B)->ProcessSequence<P1, P2> { \ */
- /* return ProcessSequence<P1, P2>(A.GetRef(), B.GetRef()); \ */
- /* } */
-
- /* /\*template <typename T1, typename T2> \ */
- /* inline ProcessSequence<T1, T2> operator%(C1<T1>& A, C2<T2>& B) { \ */
- /* return ProcessSequence<T1, T2>(A.GetRef(), B.GetRef()); \ */
- /* }*\/ */
-
- /* OPSEQ(BaseProcess, BaseProcess) */
- /* OPSEQ(BaseProcess, InteractionProcess) */
- /* OPSEQ(BaseProcess, ContinuousProcess) */
- /* OPSEQ(BaseProcess, DecayProcess) */
- /* OPSEQ(ContinuousProcess, BaseProcess) */
- /* OPSEQ(ContinuousProcess, InteractionProcess) */
- /* OPSEQ(ContinuousProcess, ContinuousProcess) */
- /* OPSEQ(ContinuousProcess, DecayProcess) */
- /* OPSEQ(InteractionProcess, BaseProcess) */
- /* OPSEQ(InteractionProcess, InteractionProcess) */
- /* OPSEQ(InteractionProcess, ContinuousProcess) */
- /* OPSEQ(InteractionProcess, DecayProcess) */
- /* OPSEQ(DecayProcess, BaseProcess) */
- /* OPSEQ(DecayProcess, InteractionProcess) */
- /* OPSEQ(DecayProcess, ContinuousProcess) */
- /* OPSEQ(DecayProcess, DecayProcess) */
-
/// marker to identify objectas ProcessSequence
template <typename A, typename B>
struct is_process_sequence<corsika::process::ProcessSequence<A, B> > {
diff --git a/Framework/Units/PhysicalUnits.h b/Framework/Units/PhysicalUnits.h
index 262e55d1c52070e36add9404f0d2ba67bae347bc..eff077af38644180888c02fa0c5d6eafaf1cf455 100644
--- a/Framework/Units/PhysicalUnits.h
+++ b/Framework/Units/PhysicalUnits.h
@@ -79,21 +79,22 @@ namespace corsika::units::si {
/**
* @file PhysicalUnits
*
- * Define _XeV literals, alowing 10_GeV in the code.
- * Define _barn literal
*/
namespace phys {
namespace units {
namespace literals {
+
+ /**
+ * Define new _XeV literals, alowing 10_GeV in the code.
+ * Define new _barn literal
+ */
+
QUANTITY_DEFINE_SCALING_LITERALS(eV, hepenergy_d, 1)
QUANTITY_DEFINE_SCALING_LITERALS(barn, corsika::units::si::sigma_d,
magnitude(corsika::units::constants::barn))
- // QUANTITY_DEFINE_SCALING_LITERALS(Ns, corsika::units::si::momentum_d,
- // magnitude(1_m * 1_kg / 1_s))
-
} // namespace literals
} // namespace units
} // namespace phys
diff --git a/Framework/Utilities/COMBoost.cc b/Framework/Utilities/COMBoost.cc
index 1bf8f839f1818fb12aa5c813ee166504e3351bd1..effb02ebfabe972799668bf95172c6967833ee16 100644
--- a/Framework/Utilities/COMBoost.cc
+++ b/Framework/Utilities/COMBoost.cc
@@ -8,16 +8,20 @@
* the license.
*/
+#include <corsika/geometry/CoordinateSystem.h>
+#include <corsika/geometry/Vector.h>
+#include <corsika/units/PhysicalUnits.h>
#include <corsika/utl/COMBoost.h>
using namespace corsika::utl;
using namespace corsika::units::si;
-COMBoost::COMBoost(HEPEnergyType eProjectile, COMBoost::MomentumVector const& pProjectile,
- HEPMassType mTarget)
+template <typename FourVector>
+COMBoost<FourVector>::COMBoost(const FourVector& Pprojectile, const HEPMassType massTarget)
: fRotation(Eigen::Matrix3d::Identity())
- , fCS(pProjectile.GetCoordinateSystem()) {
+ , fCS(Pprojectile.GetSpaceLikeComponents().GetCoordinateSystem()) {
// calculate matrix for rotating pProjectile to z-axis first
+ auto const pProjectile = Pprojectile.GetSpaceLikeComponents();
auto const pProjNorm = pProjectile.norm();
auto const a = (pProjectile / pProjNorm).GetComponents().eVector;
@@ -37,47 +41,77 @@ COMBoost::COMBoost(HEPEnergyType eProjectile, COMBoost::MomentumVector const& pP
}
// calculate boost
- double const x = pProjNorm / (eProjectile + mTarget);
+ double const beta =
+ pProjNorm / (Pprojectile.GetTimeLikeComponent() + massTarget);
- /* Accurracy matters here, x = 1 - epsilon for ultra-relativistic boosts */
- double const coshEta = 1 / std::sqrt((1 + x) * (1 - x));
- //~ double const coshEta = 1 / std::sqrt((1-x*x));
- double const sinhEta = -x * coshEta;
+ /* Accurracy matters here, beta = 1 - epsilon for ultra-relativistic boosts */
+ double const coshEta = 1 / std::sqrt((1 + beta) * (1 - beta));
+ //~ double const coshEta = 1 / std::sqrt((1-beta*beta));
+ double const sinhEta = -beta * coshEta;
+
+ std::cout << "COMBoost (1-beta)=" << 1 - beta << " gamma=" << coshEta << std::endl;
fBoost << coshEta, sinhEta, sinhEta, coshEta;
fInverseBoost << coshEta, -sinhEta, -sinhEta, coshEta;
}
-std::tuple<HEPEnergyType, corsika::geometry::QuantityVector<hepmomentum_d>>
-COMBoost::toCoM(HEPEnergyType E, COMBoost::MomentumVector p) const {
- corsika::geometry::QuantityVector<hepmomentum_d> pComponents = p.GetComponents(fCS);
+template <typename FourVector>
+FourVector COMBoost<FourVector>::toCoM(const FourVector& p) const {
+ auto pComponents = p.GetSpaceLikeComponents().GetComponents(fCS);
Eigen::Vector3d eVecRotated = fRotation * pComponents.eVector;
Eigen::Vector2d lab;
- lab << (E * (1 / 1_GeV)), (eVecRotated(2) * (1 / 1_GeV).magnitude());
+ lab << (p.GetTimeLikeComponent() * (1 / 1_GeV)),
+ (eVecRotated(2) * (1 / 1_GeV).magnitude());
auto const boostedZ = fBoost * lab;
auto const E_CoM = boostedZ(0) * 1_GeV;
eVecRotated(2) = boostedZ(1) * (1_GeV).magnitude();
- return std::make_tuple(E_CoM,
- corsika::geometry::QuantityVector<hepmomentum_d>{eVecRotated});
+ return FourVector(E_CoM, corsika::geometry::Vector<hepmomentum_d>(fCS, eVecRotated));
}
-std::tuple<HEPEnergyType, COMBoost::MomentumVector> COMBoost::fromCoM(
- HEPEnergyType E,
- corsika::geometry::QuantityVector<units::si::hepmomentum_d> pCoM) const {
+template <typename FourVector>
+FourVector COMBoost<FourVector>::fromCoM(const FourVector& p) const {
Eigen::Vector2d com;
- com << (E * (1 / 1_GeV)), (pCoM.eVector(2) * (1 / 1_GeV).magnitude());
+ com << (p.GetTimeLikeComponent() * (1 / 1_GeV)),
+ (p.GetSpaceLikeComponents().GetComponents().eVector(2) * (1 / 1_GeV).magnitude());
+
+ std::cout << "COMBoost::fromCoM Ecm=" << p.GetTimeLikeComponent() / 1_GeV << "GeV, "
+ << " pcm=" << p.GetSpaceLikeComponents().GetComponents() / 1_GeV << "GeV"
+ << std::endl;
auto const boostedZ = fInverseBoost * com;
- auto const E_CoM = boostedZ(0) * 1_GeV;
+ auto const E_lab = boostedZ(0) * 1_GeV;
- auto pLab = pCoM;
+ auto pLab = p.GetSpaceLikeComponents().GetComponents();
pLab.eVector(2) = boostedZ(1) * (1_GeV).magnitude();
pLab.eVector = fRotation.transpose() * pLab.eVector;
- return std::make_tuple(E_CoM, MomentumVector(fCS, pLab));
+ std::cout << "COMBoost::fromCoM --> Elab=" << E_lab / 1_GeV << "GeV, "
+ << " pcm=" << pLab / 1_GeV << "GeV" << std::endl;
+
+ return FourVector(E_lab, corsika::geometry::Vector(fCS, pLab));
}
+
+/*
+ Here we instantiate all physically meaningful versions of COMBoost
+ */
+
+#include <corsika/geometry/FourVector.h>
+
+namespace corsika::utl {
+
+ using corsika::geometry::FourVector;
+ using corsika::geometry::Vector;
+ using namespace corsika::units;
+
+ // for normal HEP energy/momentum units in [GeV]
+ template class COMBoost<FourVector<HEPEnergyType, Vector<hepmomentum_d>>>;
+ // template class COMBoost<FourVector<HEPEnergyType&, Vector<hepmomentum_d>&>>;
+
+ // template class COMBoost<FourVector<TimeType, Vector<length_d>>>;
+ // template class COMBoost<FourVector<TimeType&, Vector<length_d>&>>;
+} // namespace corsika::utl
diff --git a/Framework/Utilities/COMBoost.h b/Framework/Utilities/COMBoost.h
index fa57d8f3fa56c6c34432ed2dbdbbec3db14c5a14..46e808c71aeb57d384fe09fc65d8ce898904e28b 100644
--- a/Framework/Utilities/COMBoost.h
+++ b/Framework/Utilities/COMBoost.h
@@ -12,33 +12,32 @@
#define _include_corsika_utilties_comboost_h_
#include <corsika/geometry/CoordinateSystem.h>
-#include <corsika/geometry/QuantityVector.h>
-#include <corsika/geometry/Vector.h>
#include <corsika/units/PhysicalUnits.h>
+
#include <Eigen/Dense>
-#include <tuple>
namespace corsika::utl {
+
+ /**
+ This utility class handles Lorentz boost between different
+ referenence frames, using FourVectors.
+ */
+
+ template <typename FourVector>
class COMBoost {
Eigen::Matrix3d fRotation;
Eigen::Matrix2d fBoost, fInverseBoost;
corsika::geometry::CoordinateSystem const& fCS;
- using MomentumVector = corsika::geometry::Vector<units::si::hepmomentum_d>;
public:
- //! construct a COMBoost given energy and momentum of projectile and mass of target
- COMBoost(units::si::HEPEnergyType eProjectile, MomentumVector const& pProjectile,
- units::si::HEPMassType mTarget);
+ //! construct a COMBoost given four-vector of prjectile and mass of target
+ COMBoost(const FourVector& Pprojectile, const corsika::units::si::HEPEnergyType massTarget);
//! transforms a 4-momentum from lab frame to the center-of-mass frame
- std::tuple<units::si::HEPEnergyType,
- geometry::QuantityVector<units::si::hepmomentum_d>>
- toCoM(units::si::HEPEnergyType E, MomentumVector p) const;
+ FourVector toCoM(const FourVector& p) const;
//! transforms a 4-momentum from the center-of-mass frame back to lab frame
- std::tuple<units::si::HEPEnergyType, MomentumVector> fromCoM(
- units::si::HEPEnergyType E,
- geometry::QuantityVector<units::si::hepmomentum_d> pCoM) const;
+ FourVector fromCoM(const FourVector& pCoM) const;
};
} // namespace corsika::utl
diff --git a/Framework/Utilities/testCOMBoost.cc b/Framework/Utilities/testCOMBoost.cc
index 6b3b950518bd5a488910bbbd3538f93f1cda8721..712c4996c89a986ac3449c5b4d24ccbaf38e316c 100644
--- a/Framework/Utilities/testCOMBoost.cc
+++ b/Framework/Utilities/testCOMBoost.cc
@@ -12,7 +12,9 @@
// cpp file
#include <catch2/catch.hpp>
+#include <corsika/geometry/FourVector.h>
#include <corsika/geometry/RootCoordinateSystem.h>
+#include <corsika/geometry/Vector.h>
#include <corsika/units/PhysicalUnits.h>
#include <corsika/utl/COMBoost.h>
#include <iostream>
@@ -29,53 +31,157 @@ TEST_CASE("boosts") {
CoordinateSystem& rootCS =
RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
+ // helper function for energy-momentum
// relativistic energy
auto energy = [](HEPMassType m, Vector<hepmomentum_d> const& p) {
return sqrt(m * m + p.squaredNorm());
};
-
- // mandelstam-s
+
+ // helper function for mandelstam-s
auto s = [](HEPEnergyType E, QuantityVector<hepmomentum_d> const& p) {
return E * E - p.squaredNorm();
};
- // define projectile kinematics in lab frame
- HEPMassType const projectileMass = 1._GeV;
- Vector<hepmomentum_d> pProjectileLab{rootCS, {0_GeV, 1_PeV, 0_GeV}};
- HEPEnergyType const eProjectileLab = energy(projectileMass, pProjectileLab);
-
// define target kinematics in lab frame
HEPMassType const targetMass = 1_GeV;
Vector<hepmomentum_d> pTargetLab{rootCS, {0_eV, 0_eV, 0_eV}};
HEPEnergyType const eTargetLab = energy(targetMass, pTargetLab);
- // define boost to com frame
- COMBoost boost(eProjectileLab, pProjectileLab, targetMass);
-
- // boost projecticle
- auto const [eProjectileCoM, pProjectileCoM] =
- boost.toCoM(eProjectileLab, pProjectileLab);
-
- // boost target
- auto const [eTargetCoM, pTargetCoM] = boost.toCoM(eTargetLab, pTargetLab);
-
- // sum of momenta in CoM, should be 0
- auto const sumPCoM = pProjectileCoM + pTargetCoM;
- CHECK(sumPCoM[2] / 1_GeV == Approx(0).margin(absMargin));
-
- // mandelstam-s should be invariant under transformation
- CHECK(s(eProjectileLab + eTargetLab,
- pProjectileLab.GetComponents() + pTargetLab.GetComponents()) /
- 1_GeV / 1_GeV ==
- Approx(s(eProjectileCoM + eTargetCoM, pProjectileCoM + pTargetCoM) / 1_GeV /
- 1_GeV));
-
- // boost back...
- auto const [eProjectileBack, pProjectileBack] =
- boost.fromCoM(eProjectileCoM, pProjectileCoM);
-
- // ...should yield original values before the boosts
- CHECK(eProjectileBack / eProjectileLab == Approx(1));
- CHECK((pProjectileBack - pProjectileLab).norm() / pProjectileLab.norm() ==
+ /*
+ General tests check the interface and basic operation
+ */
+
+ SECTION("General tests") {
+
+ // define projectile kinematics in lab frame
+ HEPMassType const projectileMass = 1._GeV;
+ Vector<hepmomentum_d> pProjectileLab{rootCS, {0_GeV, 1_PeV, 0_GeV}};
+ HEPEnergyType const eProjectileLab = energy(projectileMass, pProjectileLab);
+ const FourVector PprojLab(eProjectileLab, pProjectileLab);
+
+ // define boost to com frame
+ COMBoost boost(PprojLab, targetMass);
+
+ // boost projecticle
+ auto const PprojCoM = boost.toCoM(PprojLab);
+
+ // boost target
+ auto const PtargCoM = boost.toCoM(FourVector(targetMass, pTargetLab));
+
+ // sum of momenta in CoM, should be 0
+ auto const sumPCoM =
+ PprojCoM.GetSpaceLikeComponents() + PtargCoM.GetSpaceLikeComponents();
+ CHECK(sumPCoM.norm() / 1_GeV == Approx(0).margin(absMargin));
+
+ // mandelstam-s should be invariant under transformation
+ CHECK(s(eProjectileLab + eTargetLab,
+ pProjectileLab.GetComponents() + pTargetLab.GetComponents()) /
+ 1_GeV / 1_GeV ==
+ Approx(s(PprojCoM.GetTimeLikeComponent() + PtargCoM.GetTimeLikeComponent(),
+ PprojCoM.GetSpaceLikeComponents().GetComponents() +
+ PtargCoM.GetSpaceLikeComponents().GetComponents()) /
+ 1_GeV / 1_GeV));
+
+ // boost back...
+ auto const PprojBack = boost.fromCoM(PprojCoM);
+
+ // ...should yield original values before the boosts
+ CHECK(PprojBack.GetTimeLikeComponent() / PprojLab.GetTimeLikeComponent() ==
+ Approx(1));
+ CHECK(
+ (PprojBack.GetSpaceLikeComponents() - PprojLab.GetSpaceLikeComponents()).norm() /
+ PprojLab.GetSpaceLikeComponents().norm() ==
Approx(0).margin(absMargin));
+ }
+
+ /*
+ special case: projectile along -z
+ */
+
+ SECTION("Test boost along z-axis") {
+
+ // define projectile kinematics in lab frame
+ HEPMassType const projectileMass = 1_GeV;
+ Vector<hepmomentum_d> pProjectileLab{rootCS, {0_GeV, 0_PeV, -1_PeV}};
+ HEPEnergyType const eProjectileLab = energy(projectileMass, pProjectileLab);
+ const FourVector PprojLab(eProjectileLab, pProjectileLab);
+
+ // define boost to com frame
+ COMBoost boost(PprojLab, targetMass);
+
+ // boost projecticle
+ auto const PprojCoM = boost.toCoM(PprojLab);
+
+ // boost target
+ auto const PtargCoM = boost.toCoM(FourVector(targetMass, pTargetLab));
+
+ // sum of momenta in CoM, should be 0
+ auto const sumPCoM =
+ PprojCoM.GetSpaceLikeComponents() + PtargCoM.GetSpaceLikeComponents();
+ CHECK(sumPCoM.norm() / 1_GeV == Approx(0).margin(absMargin));
+ }
+
+ /*
+ special case: projectile with arbitrary direction
+ */
+
+ SECTION("Test boost along tilted axis") {
+
+ const HEPMomentumType P0 = 1_PeV;
+ double theta = 33.;
+ double phi = -10.;
+ auto momentumComponents = [](double theta, double phi, HEPMomentumType ptot) {
+ return std::make_tuple(ptot * sin(theta) * cos(phi), ptot * sin(theta) * sin(phi),
+ -ptot * cos(theta));
+ };
+ auto const [px, py, pz] =
+ momentumComponents(theta / 180. * M_PI, phi / 180. * M_PI, P0);
+
+ // define projectile kinematics in lab frame
+ HEPMassType const projectileMass = 1_GeV;
+ Vector<hepmomentum_d> pProjectileLab(rootCS, {px, py, pz});
+ HEPEnergyType const eProjectileLab = energy(projectileMass, pProjectileLab);
+ const FourVector PprojLab(eProjectileLab, pProjectileLab);
+
+ // define boost to com frame
+ COMBoost boost(PprojLab, targetMass);
+
+ // boost projecticle
+ auto const PprojCoM = boost.toCoM(PprojLab);
+
+ // boost target
+ auto const PtargCoM = boost.toCoM(FourVector(targetMass, pTargetLab));
+
+ // sum of momenta in CoM, should be 0
+ auto const sumPCoM =
+ PprojCoM.GetSpaceLikeComponents() + PtargCoM.GetSpaceLikeComponents();
+ CHECK(sumPCoM.norm() / 1_GeV == Approx(0).margin(absMargin));
+ }
+
+ /*
+ test the ultra-high energy behaviour: E=ZeV
+ */
+
+ SECTION("High energy") {
+ // define projectile kinematics in lab frame
+ HEPMassType const projectileMass = 1_GeV;
+ HEPMomentumType P0 = 1_ZeV;
+ Vector<hepmomentum_d> pProjectileLab{rootCS, {0_GeV, 0_PeV, -P0}};
+ HEPEnergyType const eProjectileLab = energy(projectileMass, pProjectileLab);
+ const FourVector PprojLab(eProjectileLab, pProjectileLab);
+
+ // define boost to com frame
+ COMBoost boost(PprojLab, targetMass);
+
+ // boost projecticle
+ auto const PprojCoM = boost.toCoM(PprojLab);
+
+ // boost target
+ auto const PtargCoM = boost.toCoM(FourVector(targetMass, pTargetLab));
+
+ // sum of momenta in CoM, should be 0
+ auto const sumPCoM =
+ PprojCoM.GetSpaceLikeComponents() + PtargCoM.GetSpaceLikeComponents();
+ CHECK(sumPCoM.norm() / P0 == Approx(0).margin(absMargin)); // MAKE RELATIVE CHECK
+ }
}
diff --git a/Processes/Sibyll/Interaction.h b/Processes/Sibyll/Interaction.h
index f35b244f854dd8436f97b9e0791b7f64fbc5ecbc..288e724fe086c05adaa280a4fe385be3fdb3e27f 100644
--- a/Processes/Sibyll/Interaction.h
+++ b/Processes/Sibyll/Interaction.h
@@ -16,6 +16,7 @@
#include <corsika/environment/Environment.h>
#include <corsika/environment/NuclearComposition.h>
+#include <corsika/geometry/FourVector.h>
#include <corsika/particles/ParticleProperties.h>
#include <corsika/process/sibyll/ParticleConversion.h>
#include <corsika/process/sibyll/SibStack.h>
@@ -97,11 +98,11 @@ namespace corsika::process::sibyll {
corsika::particles::Neutron::GetMass());
// FOR NOW: assume target is at rest
- MomentumVector pTarget(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
+ MomentumVector pTarget(rootCS, {0_GeV, 0_GeV, 0_GeV});
// total momentum and energy
HEPEnergyType Elab = p.GetEnergy() + nucleon_mass;
- MomentumVector pTotLab(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
+ MomentumVector pTotLab(rootCS, {0_GeV, 0_GeV, 0_GeV});
pTotLab += p.GetMomentum();
pTotLab += pTarget;
auto const pTotLabNorm = pTotLab.norm();
@@ -165,6 +166,11 @@ namespace corsika::process::sibyll {
return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
}
+ /**
+ In this function SIBYLL is called to produce one event. The
+ event is copied (and boosted) into the shower lab frame.
+ */
+
template <typename Particle, typename Stack>
corsika::process::EProcessReturn DoInteraction(Particle& p, Stack& s) {
@@ -194,16 +200,11 @@ namespace corsika::process::sibyll {
// FOR NOW: target is always at rest
const auto eTargetLab = 0_GeV + nucleon_mass;
const auto pTargetLab = MomentumVector(rootCS, 0_GeV, 0_GeV, 0_GeV);
+ const FourVector PtargLab(eTargetLab, pTargetLab);
// define projectile
- auto const pProjectileLab = p.GetMomentum();
HEPEnergyType const eProjectileLab = p.GetEnergy();
-
- // define target kinematics in lab frame
- HEPMassType const targetMass = nucleon_mass;
- // define boost to and from CoM frame
- // CoM frame definition in Sibyll projectile: +z
- COMBoost const boost(eProjectileLab, pProjectileLab, targetMass);
+ auto const pProjectileLab = p.GetMomentum();
cout << "Interaction: ebeam lab: " << eProjectileLab / 1_GeV << endl
<< "Interaction: pbeam lab: " << pProjectileLab.GetComponents() / 1_GeV
@@ -212,18 +213,28 @@ namespace corsika::process::sibyll {
<< "Interaction: ptarget lab: " << pTargetLab.GetComponents() / 1_GeV
<< endl;
+ const FourVector PprojLab(eProjectileLab, pProjectileLab);
+
+ // define target kinematics in lab frame
+ // define boost to and from CoM frame
+ // CoM frame definition in Sibyll projectile: +z
+ COMBoost const boost(PprojLab, nucleon_mass);
+
// just for show:
// boost projecticle
- auto const [eProjectileCoM, pProjectileCoM] =
- boost.toCoM(eProjectileLab, pProjectileLab);
+ auto const PprojCoM = boost.toCoM(PprojLab);
// boost target
- auto const [eTargetCoM, pTargetCoM] = boost.toCoM(eTargetLab, pTargetLab);
+ auto const PtargCoM = boost.toCoM(PtargLab);
- cout << "Interaction: ebeam CoM: " << eProjectileCoM / 1_GeV << endl
- << "Interaction: pbeam CoM: " << pProjectileCoM / 1_GeV << endl;
- cout << "Interaction: etarget CoM: " << eTargetCoM / 1_GeV << endl
- << "Interaction: ptarget CoM: " << pTargetCoM / 1_GeV << endl;
+ cout << "Interaction: ebeam CoM: " << PprojCoM.GetTimeLikeComponent() / 1_GeV
+ << endl
+ << "Interaction: pbeam CoM: "
+ << PprojCoM.GetSpaceLikeComponents().GetComponents() / 1_GeV << endl;
+ cout << "Interaction: etarget CoM: " << PtargCoM.GetTimeLikeComponent() / 1_GeV
+ << endl
+ << "Interaction: ptarget CoM: "
+ << PtargCoM.GetSpaceLikeComponents().GetComponents() / 1_GeV << endl;
cout << "Interaction: position of interaction: " << pOrig.GetCoordinates()
<< endl;
@@ -245,8 +256,7 @@ namespace corsika::process::sibyll {
*/
#warning reading interaction cross section again, should not be necessary
auto const& compVec = mediumComposition.GetComponents();
- std::vector<si::CrossSectionType> cross_section_of_components(
- mediumComposition.GetComponents().size());
+ std::vector<si::CrossSectionType> cross_section_of_components(compVec.size());
for (size_t i = 0; i < compVec.size(); ++i) {
auto const targetId = compVec[i];
@@ -279,7 +289,8 @@ namespace corsika::process::sibyll {
<< " Ecm(GeV): " << Ecm / 1_GeV << std::endl;
if (eProjectileLab < 8.5_GeV || Ecm < 10_GeV) {
std::cout << "Interaction: "
- << " DoInteraction: should have dropped particle.." << std::endl;
+ << " DoInteraction: should have dropped particle.. "
+ << "THIS IS AN ERROR" << std::endl;
// p.Delete(); delete later... different process
} else {
fCount++;
@@ -310,15 +321,15 @@ namespace corsika::process::sibyll {
if (psib.HasDecayed()) continue;
// transform energy to lab. frame
- auto const pCoM = psib.GetMomentum().GetComponents();
+ auto const pCoM = psib.GetMomentum();
HEPEnergyType const eCoM = psib.GetEnergy();
- auto const [eLab, pLab] = boost.fromCoM(eCoM, pCoM);
+ auto const Plab = boost.fromCoM(FourVector(eCoM, pCoM));
// add to corsika stack
auto pnew = s.NewParticle();
pnew.SetPID(process::sibyll::ConvertFromSibyll(psib.GetPID()));
- pnew.SetEnergy(eLab);
- pnew.SetMomentum(pLab);
+ pnew.SetEnergy(Plab.GetTimeLikeComponent());
+ pnew.SetMomentum(Plab.GetSpaceLikeComponents());
pnew.SetPosition(pOrig);
pnew.SetTime(tOrig);
diff --git a/ThirdParty/phys/units/quantity.hpp b/ThirdParty/phys/units/quantity.hpp
index 890061b950a0bb84903a271115eed9e89016b3e8..bec55318685ffb0e69bfeeb19260ef20e6d1eed1 100644
--- a/ThirdParty/phys/units/quantity.hpp
+++ b/ThirdParty/phys/units/quantity.hpp
@@ -438,7 +438,7 @@ namespace phys {
typedef dimensions<0, 0, 0, 0, 1, 0, 0, 0> thermodynamic_temperature_d;
typedef dimensions<0, 0, 0, 0, 0, 1, 0, 0> amount_of_substance_d;
typedef dimensions<0, 0, 0, 0, 0, 0, 1, 0> luminous_intensity_d;
- typedef dimensions<0, 0, 0, 0, 0, 0, 0, 1> hepenergy_d;
+ typedef dimensions<0, 0, 0, 0, 0, 0, 0, 1> hepenergy_d; // this is not an SI unit !
// Addition operators