diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index a76ce24b4b2dcbb1f70a583a296b00372e90ff96..2cdcd21991fc83c1787beae36ae46dc8d59ab472 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -237,7 +237,7 @@ int main() {
// setup particle stack, and add primary particle
setup::Stack stack;
stack.Clear();
- const hep::EnergyType E0 = 100_TeV;
+ const hep::EnergyType E0 = 100_GeV;
double theta = 0.;
double phi = 0.;
{
diff --git a/Framework/Utilities/COMBoost.cc b/Framework/Utilities/COMBoost.cc
index 4e8a9688ea07cb6ccd8a09e3f49d3f525effc360..42d5995d2e54c01ef0ca5468dd98b1005d70fb3a 100644
--- a/Framework/Utilities/COMBoost.cc
+++ b/Framework/Utilities/COMBoost.cc
@@ -18,17 +18,23 @@ COMBoost::COMBoost(EnergyType eProjectile, COMBoost::MomentumVector const& pProj
: fRotation(Eigen::Matrix3d::Identity())
, fCS(pProjectile.GetCoordinateSystem()) {
// calculate matrix for rotating pProjectile to z-axis first
- // TODO: handle the case when pProjectile ~ (0, 0, -1)
auto const pProjNorm = pProjectile.norm();
auto const a = (pProjectile / pProjNorm).GetComponents().eVector;
- Eigen::Vector3d const b{0, 0, 1};
- auto const v = a.cross(b);
+ if (a(0) == 0 && a(1) == 0) {
+ // if pProjectile ~ (0, 0, -1), the standard formula for the rotation matrix breaks
+ // down but we can easily define a suitable rotation manually. We just need some SO(3)
+ // matrix that reverses the z-axis and I like this one:
+ fRotation << 1, 0, 0, 0, -1, 0, 0, 0, -1;
+ } else {
+ Eigen::Vector3d const b{0, 0, 1};
+ auto const v = a.cross(b);
- Eigen::Matrix3d vHat;
- vHat << 0, -v(2), v(1), v(2), 0, -v(0), -v(1), v(0), 0;
+ Eigen::Matrix3d vHat;
+ vHat << 0, -v(2), v(1), v(2), 0, -v(0), -v(1), v(0), 0;
- fRotation += vHat + vHat * vHat / (1 + a.dot(b));
+ fRotation += vHat + vHat * vHat / (1 + a.dot(b));
+ }
// calculate boost
double const x = pProjNorm * units::constants::c /
diff --git a/Framework/Utilities/testCOMBoost.cc b/Framework/Utilities/testCOMBoost.cc
index faac8c6d112e26448502910cf516eeaee949a674..3af2585630be330acd547c1087b55c174732fcc0 100644
--- a/Framework/Utilities/testCOMBoost.cc
+++ b/Framework/Utilities/testCOMBoost.cc
@@ -41,41 +41,46 @@ TEST_CASE("boosts") {
// define projectile kinematics in lab frame
MassType const projectileMass = 1._GeV / cSquared;
- Vector<momentum_d> pProjectileLab{rootCS, {0_GeV / c, 1_PeV / c, 0_GeV / c}};
- EnergyType const eProjectileLab = energy(projectileMass, pProjectileLab);
-
- // define target kinematics in lab frame
- MassType const targetMass = 1_GeV / cSquared;
- Vector<momentum_d> pTargetLab{rootCS, {0_Ns, 0_Ns, 0_Ns}};
- EnergyType 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 / c) == Approx(0).margin(absMargin));
-
- // mandelstam-s should be invariant under transformation
- CHECK(s(eProjectileLab + eTargetLab,
- pProjectileLab.GetComponents() + pTargetLab.GetComponents()) /
- (1_GeV / c) / (1_GeV / c) ==
- Approx(s(eProjectileCoM + eTargetCoM, pProjectileCoM + pTargetCoM) / (1_GeV / c) /
- (1_GeV / c)));
-
- // 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() ==
- Approx(0).margin(absMargin));
+ std::vector<Vector<momentum_d>> labProjectiles{
+ {rootCS, {0_GeV / c, 1_PeV / c, 0_GeV / c}}, // standard case
+ {rootCS, {0_GeV / c, 0_GeV / c, -1_GeV / c}}}; // "special" case
+
+ for (auto const& pProjectileLab : labProjectiles) {
+ EnergyType const eProjectileLab = energy(projectileMass, pProjectileLab);
+
+ // define target kinematics in lab frame
+ MassType const targetMass = 1_GeV / cSquared;
+ Vector<momentum_d> pTargetLab{rootCS, {0_Ns, 0_Ns, 0_Ns}};
+ EnergyType 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 / c) == Approx(0).margin(absMargin));
+
+ // mandelstam-s should be invariant under transformation
+ CHECK(s(eProjectileLab + eTargetLab,
+ pProjectileLab.GetComponents() + pTargetLab.GetComponents()) /
+ (1_GeV / c) / (1_GeV / c) ==
+ Approx(s(eProjectileCoM + eTargetCoM, pProjectileCoM + pTargetCoM) /
+ (1_GeV / c) / (1_GeV / c)));
+
+ // 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() ==
+ Approx(0).margin(absMargin));
+ }
}
diff --git a/Processes/Sibyll/CMakeLists.txt b/Processes/Sibyll/CMakeLists.txt
index 9224d9fe5f557aa8a9d2bd8e1973a52cea4c8830..b1635360e2125e3b5fa8da64bb82804585cf34a6 100644
--- a/Processes/Sibyll/CMakeLists.txt
+++ b/Processes/Sibyll/CMakeLists.txt
@@ -65,6 +65,7 @@ set_target_properties (
target_link_libraries (
ProcessSibyll
CORSIKAparticles
+ CORSIKAutilities
CORSIKAunits
CORSIKAthirdparty
CORSIKAgeometry
diff --git a/Processes/Sibyll/Interaction.h b/Processes/Sibyll/Interaction.h
index e9edec807e1ba1df784f5a6ecd15517ddbde7425..608d94d27142f1b61e907fd34b9afb37b0a75fc8 100644
--- a/Processes/Sibyll/Interaction.h
+++ b/Processes/Sibyll/Interaction.h
@@ -17,7 +17,7 @@
#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>
@@ -133,6 +133,7 @@ namespace corsika::process::sibyll {
corsika::process::EProcessReturn DoInteraction(Particle& p, Stack& s) {
using namespace corsika::units;
+ using namespace corsika::utl;
using namespace corsika::units::hep;
using namespace corsika::units::si;
using namespace corsika::geometry;
@@ -228,6 +229,26 @@ namespace corsika::process::sibyll {
std::cout << "Interaction: "
<< " DoDiscrete: gambet:" << gambet.GetComponents() << endl;
+ Vector<si::momentum_d> pProjectileLab = p.GetMomentum() / constants::c;
+ //{rootCS, {0_GeV / c, 1_PeV / c, 0_GeV / c}};
+ EnergyType const eProjectileLab = p.GetEnergy();
+ //energy(projectileMass, pProjectileLab);
+
+ // define target kinematics in lab frame
+ si::MassType const targetMass = nucleon_mass / constants::cSquared;
+ // define boost to com frame
+ COMBoost boost(eProjectileLab, pProjectileLab, targetMass);
+
+ cout << "Interaction: new boost: ebeam lab: " << eProjectileLab / 1_GeV << endl
+ << "Interaction: new boost: pbeam lab: " << pProjectileLab.GetComponents() / ( 1_GeV / constants::c ) << endl;
+
+ // 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 / constants::c ) << endl;
+
int kBeam = process::sibyll::ConvertToSibyllRaw(p.GetPID());
std::cout << "Interaction: "