fixed merging but, improved printout

Framework/Utilities/COMBoost.cc

@@ -49,8 +49,8 @@ COMBoost<FourVector>::COMBoost(const FourVector& Pprojectile, const FourVector&
   //~ 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;
-  
+  std::cout << "COMBoost (1-beta)=" << 1 - beta << " gamma=" << coshEta << std::endl;
+
   fBoost << coshEta, sinhEta, sinhEta, coshEta;
 
   fInverseBoost << coshEta, -sinhEta, -sinhEta, coshEta;
@@ -80,8 +80,9 @@ FourVector COMBoost<FourVector>::fromCoM(const FourVector& p) const {
       (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;
-  
+            << " pcm=" << p.GetSpaceLikeComponents().GetComponents() / 1_GeV << "GeV"
+            << std::endl;
+
   auto const boostedZ = fInverseBoost * com;
   auto const E_lab = boostedZ(0) * 1_GeV;
 
@@ -90,7 +91,7 @@ FourVector COMBoost<FourVector>::fromCoM(const FourVector& p) const {
   pLab.eVector = fRotation.transpose() * pLab.eVector;
 
   std::cout << "COMBoost::fromCoM --> Elab=" << E_lab / 1_GeV << "GeV, "
-	    << " pcm=" << pLab / 1_GeV << "GeV" << std::endl;
+            << " pcm=" << pLab / 1_GeV << "GeV" << std::endl;
 
   return FourVector(E_lab, corsika::geometry::Vector(fCS, pLab));
 }

Framework/Utilities/testCOMBoost.cc

@@ -41,46 +41,130 @@ TEST_CASE("boosts") {
     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);
-  const FourVector PprojLab(eProjectileLab, 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(PprojLab, FourVector(targetMass, pTargetLab));
-
-  // 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() /
+  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, FourVector(targetMass, pTargetLab));
+
+    // 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));
+  }
+
+  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, FourVector(targetMass, pTargetLab));
+
+    // 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));
+  }
+
+  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, FourVector(targetMass, pTargetLab));
+
+    // 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));
+  }
+
+  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, FourVector(targetMass, pTargetLab));
+
+    // 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
+  }
 }

Processes/Sibyll/Interaction.h

@@ -168,9 +168,9 @@ namespace corsika::process::sibyll {
 
     /**
        In this function SIBYLL is called to produce one event. The
-       event is copied (and boosted) into the shower lab frame. 
+       event is copied (and boosted) into the shower lab frame.
      */
-    
+
     template <typename Particle, typename Stack>
     corsika::process::EProcessReturn DoInteraction(Particle& p, Stack& s) {
 
@@ -184,8 +184,7 @@ namespace corsika::process::sibyll {
       const auto corsikaBeamId = p.GetPID();
       cout << "ProcessSibyll: "
            << "DoInteraction: " << corsikaBeamId << " interaction? "
-           << process::sibyll::CanInteract(corsikaBeamId)
-	   << endl;
+           << process::sibyll::CanInteract(corsikaBeamId) << endl;
       if (process::sibyll::CanInteract(corsikaBeamId)) {
         const CoordinateSystem& rootCS =
             RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
@@ -198,30 +197,24 @@ namespace corsika::process::sibyll {
         // for Sibyll is always a single nucleon
         auto constexpr nucleon_mass = 0.5 * (corsika::particles::Proton::GetMass() +
                                              corsika::particles::Neutron::GetMass());
-<<<<<<< HEAD
         // 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);
-  
+        const FourVector PtargLab(eTargetLab, pTargetLab);
+
         // define projectile
         HEPEnergyType const eProjectileLab = p.GetEnergy();
         auto const pProjectileLab = p.GetMomentum();
-        const FourVector PprojLab(eProjectileLab, pProjectileLab);
 
         cout << "Interaction: ebeam lab: " << eProjectileLab / 1_GeV << endl
              << "Interaction: pbeam lab: " << pProjectileLab.GetComponents() / 1_GeV
              << endl;
         cout << "Interaction: etarget lab: " << eTargetLab / 1_GeV << endl
              << "Interaction: ptarget lab: " << pTargetLab.GetComponents() / 1_GeV
-=======
-        const auto Etarget = 0_GeV + nucleon_mass;
-        const auto pTarget = MomentumVector(rootCS, 0_GeV, 0_GeV, 0_GeV);
-        cout << "target momentum (GeV/c): " << pTarget.GetComponents() / 1_GeV << endl;
-        cout << "beam momentum (GeV/c): " << p.GetMomentum().GetComponents() / 1_GeV
->>>>>>> 4739165629caedbec096282d6834e264b573c3b2
              << 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
@@ -234,7 +227,6 @@ namespace corsika::process::sibyll {
         // boost target
         auto const PtargCoM = boost.toCoM(PtargLab);
 
-<<<<<<< HEAD
         cout << "Interaction: ebeam CoM: " << PprojCoM.GetTimeLikeComponent() / 1_GeV
              << endl
              << "Interaction: pbeam CoM: "
@@ -264,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(
-								      compVec.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];
@@ -292,30 +283,6 @@ namespace corsika::process::sibyll {
 
         // beam id for sibyll
         const int kBeam = process::sibyll::ConvertToSibyllRaw(corsikaBeamId);
-=======
-        auto const pProjectileLab = p.GetMomentum();
-        //{rootCS, {0_GeV / c, 1_PeV / c, 0_GeV / c}};
-        HEPEnergyType const eProjectileLab = p.GetEnergy();
-        // energy(projectileMass, pProjectileLab);
-
-        // define target kinematics in lab frame
-        HEPMassType const targetMass = nucleon_mass;
-        // define boost to com frame
-        COMBoost const boost(eProjectileLab, pProjectileLab, targetMass);
-
-        cout << "Interaction: new boost: ebeam lab: " << eProjectileLab / 1_GeV << endl
-             << "Interaction: new boost: pbeam lab: "
-             << pProjectileLab.GetComponents() / 1_GeV << 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 << endl;
-
-        int kBeam = process::sibyll::ConvertToSibyllRaw(p.GetPID());
->>>>>>> 4739165629caedbec096282d6834e264b573c3b2
 
         std::cout << "Interaction: "
                   << " DoInteraction: E(GeV):" << eProjectileLab / 1_GeV
@@ -323,7 +290,7 @@ namespace corsika::process::sibyll {
         if (eProjectileLab < 8.5_GeV || Ecm < 10_GeV) {
           std::cout << "Interaction: "
                     << " DoInteraction: should have dropped particle.. "
-		    << "THIS IS AN ERROR" << std::endl;
+                    << "THIS IS AN ERROR" << std::endl;
           // p.Delete(); delete later... different process
         } else {
           fCount++;