improved testCascade

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) {

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

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;

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

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