From 8aa5eb7920fbd83183dd0c9f488f49275afc8409 Mon Sep 17 00:00:00 2001
From: nburmaso <nazar.burmasov@cern.ch>
Date: Thu, 21 Nov 2024 16:37:18 +0300
Subject: [PATCH] Translate Pythia6 status codes to Pythia8-like format

---
 include/UpcGenerator.h   |  4 ++--
 src/UpcGenerator.cpp     | 10 +++++-----
 src/UpcPythia6Helper.cpp | 19 ++++++++++++++++++-
 3 files changed, 25 insertions(+), 8 deletions(-)

diff --git a/include/UpcGenerator.h b/include/UpcGenerator.h
index 96dcedd..5b29d37 100644
--- a/include/UpcGenerator.h
+++ b/include/UpcGenerator.h
@@ -195,7 +195,7 @@ class UpcGenerator
 
   // helper struct for ROOT file output
   struct outParticle {
-    int eventNumber;
+    long int eventNumber;
     int pdgCode;
     int particleID;
     int statusID;
@@ -211,7 +211,7 @@ class UpcGenerator
   outParticle particle{};
   std::vector<TParticle> genParticles;
 
-  void writeEvent(int evt,
+  void writeEvent(long int evt,
                   const std::vector<int>& pdgs,
                   const std::vector<int>& statuses,
                   const std::vector<int>& mothers,
diff --git a/src/UpcGenerator.cpp b/src/UpcGenerator.cpp
index c75d74c..08651f2 100644
--- a/src/UpcGenerator.cpp
+++ b/src/UpcGenerator.cpp
@@ -442,7 +442,7 @@ void UpcGenerator::processInPythia(vector<int>& pdgs,
     }
     int pdg = part->GetPdgCode();
     int status = part->GetStatusCode();
-    int mother = part->GetFirstMother();
+    int mother = abs(status) == 23 ? 0 : part->GetFirstMother() + 1;
     int fDaughter = part->GetFirstDaughter();
     int lDaughter = part->GetLastDaughter();
     part->Momentum(tlVector);
@@ -518,7 +518,7 @@ bool UpcGenerator::checkKinCuts(std::vector<TLorentzVector>& particles)
   return pass;
 }
 
-void UpcGenerator::writeEvent(int evt,
+void UpcGenerator::writeEvent(long int evt,
                               const vector<int>& pdgs,
                               const vector<int>& statuses,
                               const vector<int>& mothers,
@@ -541,7 +541,7 @@ void UpcGenerator::writeEvent(int evt,
 
   if (useHepMCOut) {
     writerHepMC->writeEventInfo(evt, static_cast<int>(particles.size()), 1);
-    for (int i = 0; i < particles.size(); i++) {
+    for (int i = 0; i < particles.size(); ++i) {
       writerHepMC->writeParticleInfo(i + 1, mothers[i], pdgs[i],
                                      particles[i].Px(), particles[i].Py(), particles[i].Pz(), particles[i].E(), particles[i].M(),
                                      statuses[i]);
@@ -702,13 +702,13 @@ long int UpcGenerator::generateEvent(vector<int>& pdgs,
 
   // uniform pion decays into photon pairs
   if (procID == 111) {
-    twoPartDecayUniform(pdgs, statuses, mothers, particles, 0, 0., 22);
     twoPartDecayUniform(pdgs, statuses, mothers, particles, 1, 0., 22);
+    twoPartDecayUniform(pdgs, statuses, mothers, particles, 2, 0., 22);
   }
 
   // uniform ALP decay into two photons
   if (procID == 51) {
-    twoPartDecayUniform(pdgs, statuses, mothers, particles, 0, 0., 22);
+    twoPartDecayUniform(pdgs, statuses, mothers, particles, 1, 0., 22);
   }
 
   // fill genParticles
diff --git a/src/UpcPythia6Helper.cpp b/src/UpcPythia6Helper.cpp
index 2a3109c..2de1cff 100644
--- a/src/UpcPythia6Helper.cpp
+++ b/src/UpcPythia6Helper.cpp
@@ -45,6 +45,9 @@ void UpcPythia6Helper::process(std::vector<int>& pdgs, std::vector<int>& statuse
   }
 }
 
+// status codes:
+// https://www.pythia.org/download/pythia6/pythia6200.pdf
+// https://pythia.org/download/pdf/worksheet8153.pdf
 int UpcPythia6Helper::import(TClonesArray* particles)
 {
   int nParts = 0;
@@ -53,7 +56,21 @@ int UpcPythia6Helper::import(TClonesArray* particles)
   for (int i = 0; i < mPartHolder.size(); i++) {
     for (int j = 0; j < mPartHolder[i].GetEntriesFast(); j++) {
       auto* part = (TParticle*)mPartHolder[i].At(j);
-      int mother = part->GetFirstMother() == 0 ? 0 : part->GetFirstMother() + shift;
+      bool isPrimary = part->GetFirstMother() == 0;
+      bool isFinal = part->GetFirstDaughter() == 0;
+      int statusP6 = part->GetStatusCode();
+      int statusP8 = -1;
+      if (statusP6 == 11 && isPrimary) { // primary, decayed
+        statusP8 = -23;
+      }
+      if (statusP6 == 11 && !isFinal && !isPrimary) { // secondary, decayed
+        statusP8 = -91;
+      }
+      if (statusP6 == 1 && isFinal && !isPrimary) { // secondary, "final state"
+        statusP8 = 91;
+      }
+      part->SetStatusCode(statusP8);
+      int mother = part->GetFirstMother() == 0 ? -1 : part->GetFirstMother() + shift - 1;
       part->SetFirstMother(mother);
       new (clonesParticles[nParts]) TParticle(*part);
       nParts++;