diff --git a/Processes/QGSJetII/Interaction.cc b/Processes/QGSJetII/Interaction.cc index d1b84e5c1a11b73b2d67f6474c45a1d4d791a6ca..c7a2ffc1846e96fdeb4df54bd044ddad33ce834a 100644 --- a/Processes/QGSJetII/Interaction.cc +++ b/Processes/QGSJetII/Interaction.cc @@ -201,12 +201,10 @@ namespace corsika::process::qgsjetII { HEPEnergyType const projectileEnergyLab = vP.GetEnergy(); auto const projectileMomentumLab = vP.GetMomentum(); - int beamA = 0; + int beamA = 1; if (particles::IsNucleus(corsikaBeamId)) beamA = vP.GetNuclearA(); - const HEPEnergyType projectileEnergyLabPerNucleon = - particles::IsNucleus(corsikaBeamId) ? projectileEnergyLab / beamA - : projectileEnergyLab; + const HEPEnergyType projectileEnergyLabPerNucleon = projectileEnergyLab / beamA; cout << "Interaction: ebeam lab: " << projectileEnergyLab / 1_GeV << endl << "Interaction: pbeam lab: " << projectileMomentumLab.GetComponents() / 1_GeV @@ -315,15 +313,15 @@ namespace corsika::process::qgsjetII { Z = 0; } - const HEPMassType projectileMass = particles::GetMass(idFragm); + const HEPMassType nucleonMass = particles::GetMass(idFragm); auto momentum = geometry::Vector( zAxisFrame, corsika::geometry::QuantityVector<hepmomentum_d>{ 0.0_GeV, 0.0_GeV, - sqrt((projectileEnergyLab + projectileMass) * - (projectileEnergyLab - projectileMass))}); + sqrt((projectileEnergyLabPerNucleon + nucleonMass) * + (projectileEnergyLabPerNucleon - nucleonMass))}); - auto const energy = sqrt(momentum.squaredNorm() + square(projectileMass)); + auto const energy = sqrt(momentum.squaredNorm() + square(nucleonMass)); momentum.rebase(originalCS); // transform back into standard lab frame std::cout << "secondary fragment> id=" << idFragm << " p=" << momentum.GetComponents() << std::endl; @@ -349,7 +347,7 @@ namespace corsika::process::qgsjetII { } } - if (idFragm == particles::Code::Nucleus) { + if (idFragm == particles::Code::Nucleus) { // thus: not p or n const HEPMassType nucleusMass = particles::Proton::GetMass() * Z + particles::Neutron::GetMass() * (A - Z); auto momentum = geometry::Vector(