diff --git a/Framework/Utilities/COMBoost.cc b/Framework/Utilities/COMBoost.cc
index 3e9aafb8a3dfbfe41d3f353930e6ff548d3bff7e..48dd7b37d3a78a893a23d4b22434f265db077e30 100644
--- a/Framework/Utilities/COMBoost.cc
+++ b/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));
}
diff --git a/Framework/Utilities/testCOMBoost.cc b/Framework/Utilities/testCOMBoost.cc
index bbc6c295ceaa856b770db7f607f3417d98c2f725..b134dcc65675e554e3f06d87b0c6f690bd14db1e 100644
--- a/Framework/Utilities/testCOMBoost.cc
+++ b/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
+ }
}
diff --git a/Processes/Sibyll/Interaction.h b/Processes/Sibyll/Interaction.h
index a1b52d5b87031ca2925e9b0a92ff735f94de91d1..8d52bcf363d3d048a72fad5146b401c4ee2c2415 100644
--- a/Processes/Sibyll/Interaction.h
+++ b/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++;