diff --git a/corsika/detail/modules/epos/InteractionModel.inl b/corsika/detail/modules/epos/InteractionModel.inl
index f9e8925b7b2da6d94744a6d2afcad1c4fefdfa6a..dbaeba716a2aba70d7d3bf5a55dd617d2585609f 100644
--- a/corsika/detail/modules/epos/InteractionModel.inl
+++ b/corsika/detail/modules/epos/InteractionModel.inl
@@ -76,7 +76,7 @@ namespace corsika::epos {
::epos::aaset_(iarg);
// debug output settings
- ::epos::prnt1_.ish = 0; // debug level in epos, 0: off, 6: medium output
+ ::epos::prnt1_.ish = 0;
::epos::prnt3_.iwseed = 0; // 1: printout seeds, 0: off
::epos::files_.ifch = 6; // output unit, 6: screen
@@ -144,18 +144,18 @@ namespace corsika::epos {
inline void InteractionModel::initializeEventCoM(Code const idBeam, int const iBeamA,
int const iBeamZ, Code const idTarget,
int const iTargetA, int const iTargetZ,
- HEPEnergyType const Ecm) const {
+ HEPEnergyType const EcmNN) const {
CORSIKA_LOGGER_TRACE(logger_,
"initialize event in CoM frame!"
" Ecm={}",
- Ecm);
+ EcmNN);
::epos::lept1_.engy = -1.;
::epos::enrgy_.ecms = -1.;
::epos::enrgy_.elab = -1.;
::epos::enrgy_.ekin = -1.;
::epos::hadr1_.pnll = -1.;
- ::epos::enrgy_.ecms = Ecm / 1_GeV; // -> c.m.s. frame
+ ::epos::enrgy_.ecms = EcmNN / 1_GeV; // -> c.m.s. frame
CORSIKA_LOGGER_TRACE(logger_, "inside EPOS: Ecm={}, Elab={}", ::epos::enrgy_.ecms,
::epos::enrgy_.elab);
@@ -366,43 +366,48 @@ namespace corsika::epos {
count_ = count_ + 1;
- // define projectile
- // define projectile, in lab frame
- auto const sqrtS2 = (projectileP4 + targetP4).getNormSqr();
- HEPEnergyType const sqrtS = sqrt(sqrtS2);
- if (!isValid(projectileId, targetId, sqrtS)) {
+ // define nucleon-nucleon center-of-mass frame
+ auto const projectileP4NN =
+ projectileP4 / (is_nucleus(projectileId) ? get_nucleus_A(projectileId) : 1);
+ auto const targetP4NN =
+ targetP4 / (is_nucleus(targetId) ? get_nucleus_A(targetId) : 1);
+ auto const SNN = (projectileP4NN + targetP4NN).getNormSqr();
+ HEPEnergyType const sqrtSNN = sqrt(SNN);
+ if (!isValid(projectileId, targetId, sqrtSNN)) {
throw std::runtime_error("invalid projectile/target/energy combination.");
}
- HEPEnergyType const Elab = (sqrtS2 - static_pow<2>(get_mass(projectileId)) -
+ HEPEnergyType const Elab = (SNN - static_pow<2>(get_mass(projectileId)) -
static_pow<2>(get_mass(targetId))) /
(2 * get_mass(targetId));
// system of initial-state
- COMBoost const boost(projectileP4, targetP4);
+ COMBoost const boost(projectileP4NN, targetP4NN);
auto const& originalCS = boost.getOriginalCS();
auto const& csPrime =
boost.getRotatedCS(); // z is along the CM motion (projectile, in Cascade)
- HEPMomentumType const pLabMag =
- sqrt((Elab - get_mass(projectileId)) * (Elab + get_mass(projectileId)));
- MomentumVector pLab(csPrime, {0_eV, 0_eV, pLabMag});
-
CORSIKA_LOGGER_DEBUG(logger_,
"doInteraction: interaction, projectile id={}, E={}, p3={} ",
projectileId, projectileP4.getTimeLikeComponent(),
projectileP4.getSpaceLikeComponents());
+ CORSIKA_LOGGER_DEBUG(
+ logger_, "doInteraction: projectile per-nucleon ENN={}, p3NN={} ",
+ projectileP4NN.getTimeLikeComponent(), projectileP4NN.getSpaceLikeComponents());
CORSIKA_LOGGER_DEBUG(
logger_, "doInteraction: interaction, target id={}, E={}, p3={} ", targetId,
targetP4.getTimeLikeComponent(), targetP4.getSpaceLikeComponents());
- CORSIKA_LOGGER_DEBUG(logger_, "doInteraction: Elab={}, sqrtS={} ", Elab, sqrtS);
+ CORSIKA_LOGGER_DEBUG(logger_, "doInteraction: target per-nucleon ENN={}, p3NN={} ",
+ targetP4NN.getTimeLikeComponent(),
+ targetP4NN.getSpaceLikeComponents());
+ CORSIKA_LOGGER_DEBUG(logger_, "doInteraction: Elab={}, sqrtSNN={} ", Elab, sqrtSNN);
int beamA = 1;
int beamZ = 1;
if (is_nucleus(projectileId)) {
beamA = get_nucleus_A(projectileId);
beamZ = get_nucleus_Z(projectileId);
- CORSIKA_LOGGER_DEBUG(logger_, "A={}, Z={} ", beamA, beamZ);
+ CORSIKA_LOGGER_DEBUG(logger_, "projectile: A={}, Z={} ", beamA, beamZ);
}
// // from corsika7 interface
@@ -412,8 +417,9 @@ namespace corsika::epos {
if (is_nucleus(targetId)) {
targetA = get_nucleus_A(targetId);
targetZ = get_nucleus_Z(targetId);
+ CORSIKA_LOGGER_DEBUG(logger_, "target: A={}, Z={} ", beamA, beamZ);
}
- initializeEventCoM(projectileId, beamA, beamZ, targetId, targetA, targetZ, sqrtS);
+ initializeEventCoM(projectileId, beamA, beamZ, targetId, targetA, targetZ, sqrtSNN);
// create event
int iarg = 1;
diff --git a/corsika/detail/modules/qgsjetII/InteractionModel.inl b/corsika/detail/modules/qgsjetII/InteractionModel.inl
index 13ad2d0c814ff994a47c4bd00ec8d9632d85442d..649fb3ddb7e11192fa3ab2c0851efdd0626b17be 100644
--- a/corsika/detail/modules/qgsjetII/InteractionModel.inl
+++ b/corsika/detail/modules/qgsjetII/InteractionModel.inl
@@ -70,11 +70,14 @@ namespace corsika::qgsjetII {
}
// define projectile, in lab frame
- auto const sqrtS2 = (projectileP4 + targetP4).getNormSqr();
- auto const sqrtS = sqrt(sqrtS2);
- if (!isValid(projectileId, targetId, sqrtS)) { return CrossSectionType::zero(); }
+ auto const SNN =
+ (projectileP4 / (is_nucleus(projectileId) ? get_nucleus_A(projectileId) : 1) +
+ targetP4 / (is_nucleus(targetId) ? get_nucleus_A(targetId) : 1))
+ .getNormSqr();
+ auto const sqrtSNN = sqrt(SNN);
+ if (!isValid(projectileId, targetId, sqrtSNN)) { return CrossSectionType::zero(); }
HEPEnergyType const Elab =
- calculate_lab_energy(sqrtS2, get_mass(projectileId), get_mass(targetId));
+ calculate_lab_energy(SNN, get_mass(projectileId), get_mass(targetId));
int const iBeam = static_cast<QgsjetIIXSClassIntType>(
corsika::qgsjetII::getQgsjetIIXSCode(projectileId));
@@ -104,14 +107,25 @@ namespace corsika::qgsjetII {
projectileId, corsika::qgsjetII::canInteract(projectileId));
// define projectile, in lab frame
- auto const sqrtS2 = (projectileP4 + targetP4).getNormSqr();
- auto const sqrtS = sqrt(sqrtS2);
+ auto const SNN =
+ (projectileP4 / (is_nucleus(projectileId) ? get_nucleus_A(projectileId) : 1) +
+ targetP4 / (is_nucleus(targetId) ? get_nucleus_A(targetId) : 1))
+ .getNormSqr();
+ auto const sqrtS = sqrt(SNN);
if (!corsika::qgsjetII::canInteract(projectileId) ||
!isValid(projectileId, targetId, sqrtS)) {
throw std::runtime_error("invalid target/projectile/energy combination.");
}
- HEPEnergyType const Elab =
- calculate_lab_energy(sqrtS2, get_mass(projectileId), get_mass(targetId));
+ auto const projectileMass =
+ (is_nucleus(projectileId) ? constants::nucleonMass : get_mass(projectileId));
+ auto const targetMass =
+ (projectileId == Code::Proton
+ ? get_mass(Code::Proton)
+ : constants::nucleonMass); // qgsjet target is always proton or nucleon.
+ // always nucleon??
+
+ // lab energy/hadron
+ HEPEnergyType const Elab = calculate_lab_energy(SNN, projectileMass, targetMass);
int beamA = 0;
if (is_nucleus(projectileId)) { beamA = get_nucleus_A(projectileId); }
@@ -163,18 +177,18 @@ namespace corsika::qgsjetII {
// rest.
// system of initial-state
- COMBoost boost(projectileP4, targetP4);
+ COMBoost boost(projectileP4 / projectileMassNumber, targetP4 / targetMassNumber);
auto const& originalCS = boost.getOriginalCS();
auto const& csPrime =
boost.getRotatedCS(); // z is along the CM motion (projectile, in Cascade)
HEPMomentumType const pLabMag =
- sqrt((Elab - get_mass(projectileId)) * (Elab + get_mass(projectileId)));
+ sqrt((Elab - projectileMass) * (Elab + projectileMass));
MomentumVector pLab(csPrime, {0_eV, 0_eV, pLabMag});
- // internal QGSJetII system
- COMBoost boostInternal({Elab, pLab}, get_mass(targetId));
+ // internal QGSJetII system: hadron-nucleon lab. frame!
+ COMBoost boostInternal({Elab, pLab}, targetMass);
// fragments
QGSJetIIFragmentsStack qfs;
@@ -265,7 +279,7 @@ namespace corsika::qgsjetII {
Plab_final += pnew.getMomentum();
Elab_final += pnew.getEnergy();
}
- } // namespace corsika::qgsjetII
+ }
// secondaries
QGSJetIIStack qs;
diff --git a/tests/modules/testQGSJetII.cpp b/tests/modules/testQGSJetII.cpp
index 966928f9a06d8dfade6ddebca63c9639a7d2a90c..e9caaecf8c3341c53351e78c023a96838af2bd76 100644
--- a/tests/modules/testQGSJetII.cpp
+++ b/tests/modules/testQGSJetII.cpp
@@ -123,6 +123,8 @@ TEST_CASE("QgsjetIIInterface", "interaction,processes") {
logging::set_level(logging::level::info);
+ RNGManager<>::getInstance().registerRandomStream("qgsjet");
+
auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Oxygen);
auto const& cs = *csPtr;
[[maybe_unused]] auto const& env_dummy = env;
@@ -174,7 +176,7 @@ TEST_CASE("QgsjetIIInterface", "interaction,processes") {
corsika::qgsjetII::InteractionModel model;
model.doInteraction(view, pid, Code::Oxygen, projectileP4,
targetP4); // this also should produce some fragments
- CHECK(view.getSize() == Approx(300).margin(150)); // this is not physics validation
+ CHECK(view.getSize() == Approx(150).margin(150)); // this is not physics validation
int countFragments = 0;
for (auto const& sec : view) { countFragments += (is_nucleus(sec.getPID())); }
CHECK(countFragments == Approx(4).margin(2)); // this is not physics validation
@@ -216,7 +218,7 @@ TEST_CASE("QgsjetIIInterface", "interaction,processes") {
{sqrt(static_pow<2>(1_TeV) + static_pow<2>(Pi0::mass)),
MomentumVector{cs, 1_TeV, 0_GeV, 0_GeV}},
{Oxygen::mass, MomentumVector{cs, 0_eV, 0_eV, 0_eV}});
- CHECK(view.getSize() == Approx(10).margin(8)); // this is not physics validation
+ CHECK(view.getSize() == Approx(20).margin(20)); // this is not physics validation
}
{ // rho0 is internally converted into pi-/pi+
auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
@@ -228,7 +230,7 @@ TEST_CASE("QgsjetIIInterface", "interaction,processes") {
{sqrt(static_pow<2>(1_TeV) + static_pow<2>(Rho0::mass)),
MomentumVector{cs, 1_TeV, 0_GeV, 0_GeV}},
{Oxygen::mass, MomentumVector{cs, 0_eV, 0_eV, 0_eV}});
- CHECK(view.getSize() == Approx(25).margin(20)); // this is not physics validation
+ CHECK(view.getSize() == Approx(50).margin(50)); // this is not physics validation
}
{ // Lambda is internally converted into neutron
auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
@@ -241,7 +243,7 @@ TEST_CASE("QgsjetIIInterface", "interaction,processes") {
{sqrt(static_pow<2>(100_GeV) + static_pow<2>(Lambda0::mass)),
MomentumVector{cs, 100_GeV, 0_GeV, 0_GeV}},
{Oxygen::mass, MomentumVector{cs, 0_eV, 0_eV, 0_eV}});
- CHECK(view.getSize() == Approx(25).margin(20)); // this is not physics validation
+ CHECK(view.getSize() == Approx(50).margin(50)); // this is not physics validation
}
{ // AntiLambda is internally converted into anti neutron
auto [stackPtr, secViewPtr] = setup::testing::setup_stack(