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