From db64dd187b74c2451f4b251341041af6d2009d1a Mon Sep 17 00:00:00 2001
From: Felix Riehn <friehn@lip.pt>
Date: Wed, 31 Jul 2024 16:26:39 +0200
Subject: [PATCH] switch to nucleon-nucleon frame
---
.../detail/modules/pythia8/Interaction.inl | 25 ++++++++++---------
1 file changed, 13 insertions(+), 12 deletions(-)
diff --git a/corsika/detail/modules/pythia8/Interaction.inl b/corsika/detail/modules/pythia8/Interaction.inl
index 01e88c6d2..470cc66b2 100644
--- a/corsika/detail/modules/pythia8/Interaction.inl
+++ b/corsika/detail/modules/pythia8/Interaction.inl
@@ -46,11 +46,11 @@ namespace corsika::pythia8 {
pythia_.settings.mode("Beams:idA",
static_cast<PDGCodeIntType>(get_PDG(Code::Proton)));
pythia_.settings.mode("Beams:idB",
- static_cast<PDGCodeIntType>(get_PDG(Code::Nitrogen)));
+ static_cast<PDGCodeIntType>(get_PDG(Code::Oxygen)));
pythia_.readString("Beams:allowVariableEnergy = on");
pythia_.readString("Beams:frameType = 1"); // CoM
- pythia_.settings.parm("Beams:eCM", 10000.);
+ pythia_.settings.parm("Beams:eCM", 400_TeV/1_GeV);
// needed when regular Pythia (not Angantyr) takes over for pp
pythia_.readString("SoftQCD:all = on");
@@ -74,7 +74,7 @@ namespace corsika::pythia8 {
pythia_.readString("HadronLevel:Decay = off");
// Reduce printout and relax energy-momentum conservation.
- pythia_.readString("Print:quiet = on");
+ //pythia_.readString("Print:quiet = on");
pythia_.readString("Check:epTolErr = 0.1");
pythia_.readString("Check:epTolWarn = 0.0001");
pythia_.readString("Check:mTolErr = 0.01");
@@ -218,12 +218,12 @@ namespace corsika::pythia8 {
projectileId, corsika::pythia8::Interaction::canInteract(projectileId));
// define system (this is Nucleus-Nucleus!!)
- auto const sqrtS2 = (projectileP4 + targetP4).getNormSqr();
- HEPEnergyType const sqrtS = sqrt(sqrtS2);
+ auto const sqrtS2NN = (projectileP4 + targetP4/get_nucleus_A(targetId)).getNormSqr();
+ HEPEnergyType const sqrtSNN = sqrt(sqrtS2NN);
HEPEnergyType const eProjectileLab =
- calculate_lab_energy(sqrtS2, get_mass(projectileId), get_mass(targetId));
+ calculate_lab_energy(sqrtS2NN, get_mass(projectileId), get_mass(targetId)/get_nucleus_A(targetId));
- if (!isValid(projectileId, targetId, sqrtS)) {
+ if (!isValid(projectileId, targetId, sqrtSNN)) {
throw std::runtime_error("invalid target,projectile,energy combination.");
}
@@ -231,8 +231,8 @@ namespace corsika::pythia8 {
CORSIKA_LOG_DEBUG(
"Interaction: "
" doInteraction: E(GeV): {}"
- " Ecm(GeV): {}",
- eProjectileLab / 1_GeV, sqrtS / 1_GeV);
+ " Ecm_NN(GeV): {}",
+ eProjectileLab / 1_GeV, sqrtSNN / 1_GeV);
COMBoost const labFrameBoost{targetP4.getSpaceLikeComponents(), get_mass(targetId)};
// auto const proj4MomLab = labFrameBoost.toCoM(projectileP4);
@@ -241,9 +241,10 @@ namespace corsika::pythia8 {
// variables to talk to Pythia
double const eCM_pythia =
- sqrtS / 1_GeV; // CoM energy hadron-Nucleus or Nucleus-Nucleus
- double const idA_pythia = static_cast<int>(get_PDG(projectileId));
- double const idB_pythia = static_cast<int>(get_PDG(targetId));
+ sqrtSNN / 1_GeV; // CoM energy hadron-Nucleon
+ double const idA_pythia = static_cast<double>(get_PDG(projectileId));
+ double const idB_pythia = static_cast<double>(get_PDG(targetId));
+ CORSIKA_LOG_INFO("re-configuring pythia idA={}, idB={}, ecm={} GeV", idA_pythia, idB_pythia, eCM_pythia);
// configure event
pythia_.setBeamIDs(idA_pythia, idB_pythia);
pythia_.setKinematics(eCM_pythia);
--
GitLab