Merge branch '115-central-boost-routines-produce-NaN-in-sibyll-interface' into 'master'

Resolve "central boost routines produce NaN in sibyll interface"

Closes #115

See merge request !47

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.;
   {

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 /

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));
+  }
 }

Processes/Sibyll/CMakeLists.txt

@@ -65,6 +65,7 @@ set_target_properties (
 target_link_libraries (
   ProcessSibyll
   CORSIKAparticles
+  CORSIKAutilities
   CORSIKAunits
   CORSIKAthirdparty
   CORSIKAgeometry

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: "