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/Utilities/COMBoost.cc b/Framework/Utilities/COMBoost.cc
index 48dd7b37d3a78a893a23d4b22434f265db077e30..effb02ebfabe972799668bf95172c6967833ee16 100644
--- a/Framework/Utilities/COMBoost.cc
+++ b/Framework/Utilities/COMBoost.cc
@@ -17,7 +17,7 @@ using namespace corsika::utl;
using namespace corsika::units::si;
template <typename FourVector>
-COMBoost<FourVector>::COMBoost(const FourVector& Pprojectile, const FourVector& Ptarget)
+COMBoost<FourVector>::COMBoost(const FourVector& Pprojectile, const HEPMassType massTarget)
: fRotation(Eigen::Matrix3d::Identity())
, fCS(Pprojectile.GetSpaceLikeComponents().GetCoordinateSystem()) {
// calculate matrix for rotating pProjectile to z-axis first
@@ -42,7 +42,7 @@ COMBoost<FourVector>::COMBoost(const FourVector& Pprojectile, const FourVector&
// calculate boost
double const beta =
- pProjNorm / (Pprojectile.GetTimeLikeComponent() + Ptarget.GetTimeLikeComponent());
+ pProjNorm / (Pprojectile.GetTimeLikeComponent() + massTarget);
/* Accurracy matters here, beta = 1 - epsilon for ultra-relativistic boosts */
double const coshEta = 1 / std::sqrt((1 + beta) * (1 - beta));
diff --git a/Framework/Utilities/COMBoost.h b/Framework/Utilities/COMBoost.h
index 08dfe75adb3f743cc3625e98dfe918b07f9fe90a..46e808c71aeb57d384fe09fc65d8ce898904e28b 100644
--- a/Framework/Utilities/COMBoost.h
+++ b/Framework/Utilities/COMBoost.h
@@ -12,6 +12,7 @@
#define _include_corsika_utilties_comboost_h_
#include <corsika/geometry/CoordinateSystem.h>
+#include <corsika/units/PhysicalUnits.h>
#include <Eigen/Dense>
@@ -29,8 +30,8 @@ namespace corsika::utl {
corsika::geometry::CoordinateSystem const& fCS;
public:
- //! construct a COMBoost given energy and momentum of projectile and mass of target
- COMBoost(const FourVector& Pprojectile, const FourVector& Ptarget);
+ //! 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
FourVector toCoM(const FourVector& p) const;
diff --git a/Framework/Utilities/testCOMBoost.cc b/Framework/Utilities/testCOMBoost.cc
index b134dcc65675e554e3f06d87b0c6f690bd14db1e..712c4996c89a986ac3449c5b4d24ccbaf38e316c 100644
--- a/Framework/Utilities/testCOMBoost.cc
+++ b/Framework/Utilities/testCOMBoost.cc
@@ -31,12 +31,13 @@ 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();
};
@@ -46,6 +47,10 @@ TEST_CASE("boosts") {
Vector<hepmomentum_d> pTargetLab{rootCS, {0_eV, 0_eV, 0_eV}};
HEPEnergyType const eTargetLab = energy(targetMass, pTargetLab);
+ /*
+ General tests check the interface and basic operation
+ */
+
SECTION("General tests") {
// define projectile kinematics in lab frame
@@ -55,7 +60,7 @@ TEST_CASE("boosts") {
const FourVector PprojLab(eProjectileLab, pProjectileLab);
// define boost to com frame
- COMBoost boost(PprojLab, FourVector(targetMass, pTargetLab));
+ COMBoost boost(PprojLab, targetMass);
// boost projecticle
auto const PprojCoM = boost.toCoM(PprojLab);
@@ -89,6 +94,10 @@ TEST_CASE("boosts") {
Approx(0).margin(absMargin));
}
+ /*
+ special case: projectile along -z
+ */
+
SECTION("Test boost along z-axis") {
// define projectile kinematics in lab frame
@@ -98,7 +107,7 @@ TEST_CASE("boosts") {
const FourVector PprojLab(eProjectileLab, pProjectileLab);
// define boost to com frame
- COMBoost boost(PprojLab, FourVector(targetMass, pTargetLab));
+ COMBoost boost(PprojLab, targetMass);
// boost projecticle
auto const PprojCoM = boost.toCoM(PprojLab);
@@ -112,6 +121,10 @@ TEST_CASE("boosts") {
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;
@@ -131,7 +144,7 @@ TEST_CASE("boosts") {
const FourVector PprojLab(eProjectileLab, pProjectileLab);
// define boost to com frame
- COMBoost boost(PprojLab, FourVector(targetMass, pTargetLab));
+ COMBoost boost(PprojLab, targetMass);
// boost projecticle
auto const PprojCoM = boost.toCoM(PprojLab);
@@ -145,6 +158,10 @@ TEST_CASE("boosts") {
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;
@@ -154,7 +171,7 @@ TEST_CASE("boosts") {
const FourVector PprojLab(eProjectileLab, pProjectileLab);
// define boost to com frame
- COMBoost boost(PprojLab, FourVector(targetMass, pTargetLab));
+ COMBoost boost(PprojLab, targetMass);
// boost projecticle
auto const PprojCoM = boost.toCoM(PprojLab);
diff --git a/Processes/Sibyll/Interaction.h b/Processes/Sibyll/Interaction.h
index 8d52bcf363d3d048a72fad5146b401c4ee2c2415..288e724fe086c05adaa280a4fe385be3fdb3e27f 100644
--- a/Processes/Sibyll/Interaction.h
+++ b/Processes/Sibyll/Interaction.h
@@ -218,7 +218,7 @@ namespace corsika::process::sibyll {
// define target kinematics in lab frame
// define boost to and from CoM frame
// CoM frame definition in Sibyll projectile: +z
- COMBoost const boost(PprojLab, PtargLab);
+ COMBoost const boost(PprojLab, nucleon_mass);
// just for show:
// boost projecticle