switch to nucleon-nucleon frame

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);