From d46f3d3ade7db4dcd743c64a77e821a0181fd8ae Mon Sep 17 00:00:00 2001
From: Joe Latessa <latessa@users.noreply.github.com>
Date: Mon, 9 Oct 2023 21:08:07 -0400
Subject: [PATCH] Fix rapidities in pythia from JETSCAPE 3.6 #160

Fix rapidities in pythia IS modules and make
constants consistent #160
---
 src/framework/JetScapeConstants.h |  7 +++----
 src/initialstate/PGun.cc          | 13 ++++++++++---
 src/initialstate/PythiaGun.cc     | 13 +++++++------
 3 files changed, 20 insertions(+), 13 deletions(-)

diff --git a/src/framework/JetScapeConstants.h b/src/framework/JetScapeConstants.h
index 5f440d53..33306282 100644
--- a/src/framework/JetScapeConstants.h
+++ b/src/framework/JetScapeConstants.h
@@ -2,7 +2,7 @@
  * Copyright (c) The JETSCAPE Collaboration, 2018
  *
  * Modular, task-based framework for simulating all aspects of heavy-ion collisions
- * 
+ *
  * For the list of contributors see AUTHORS.
  *
  * Report issues at https://github.com/JETSCAPE/JETSCAPE/issues
@@ -34,11 +34,10 @@ static double Nc = 3.0;
 static double Lambda_QCD = 0.2;
 // 0.4 is the value chosen in JETSET
 
-static double fmToGeVinv = 5.0;
-/// < should be 1/0.197, but 5 helps in debugging.
-
 static const double hbarC = 0.197327053;
 
+static const double fmToGeVinv = 1.0 / hbarC;
+
 static double zeta3 = 1.20206;
 
 static double mu = 0.722;
diff --git a/src/initialstate/PGun.cc b/src/initialstate/PGun.cc
index 912f3f6a..8c1b7b66 100644
--- a/src/initialstate/PGun.cc
+++ b/src/initialstate/PGun.cc
@@ -58,7 +58,7 @@ void PGun::ExecuteTask() {
 
   double p[4], xLoc[4];
 
-  double pT, rapidity, phi;
+  double pT, rapidity, phi, pseudorapidity;
   double eta_cut = 1.0;
   double tempRand;
   const double maxN = 1.0 * RAND_MAX;
@@ -91,6 +91,13 @@ void PGun::ExecuteTask() {
   p[3] = sqrt(pT * pT + mass * mass) * sinh(rapidity);
   p[0] = sqrt(pT * pT + mass * mass) * cosh(rapidity);
 
+  double p_abs = std::sqrt(pT*pT + p[3]*p[3]);
+  if(std::abs(p_abs - p[3]) > rounding_error){
+    pseudorapidity = 0.5 * std::log((p_abs + p[3]) / (p_abs - p[3]));
+  } else {
+    JSWARN << "Particle in PGun has infinite pseudorapidity.";
+  }
+
   // Roll for a starting point
   // See: https://stackoverflow.com/questions/15039688/random-generator-from-vector-with-probability-distribution-in-c
   for (int i = 0; i <= 3; i++) {
@@ -115,9 +122,9 @@ void PGun::ExecuteTask() {
   xLoc[3] = 0.0;
   
   if (flag_useHybridHad != 1) {
-    AddParton(make_shared<Parton>(0, parID, 0, pT, rapidity, phi, p[0], xLoc));
+    AddParton(make_shared<Parton>(0, parID, 0, pT, pseudorapidity, phi, p[0], xLoc));
   } else {
-    auto ptn = make_shared<Parton>(0, parID, 0, pT, rapidity, phi, p[0], xLoc);
+    auto ptn = make_shared<Parton>(0, parID, 0, pT, pseudorapidity, phi, p[0], xLoc);
     ptn->set_color((parID > 0) ? 100 : 0);
     ptn->set_anti_color(((parID > 0) || (parID == 21)) ? 0 : 101);
     ptn->set_max_color(102);
diff --git a/src/initialstate/PythiaGun.cc b/src/initialstate/PythiaGun.cc
index 0b58756b..0f5d62fd 100644
--- a/src/initialstate/PythiaGun.cc
+++ b/src/initialstate/PythiaGun.cc
@@ -170,10 +170,10 @@ void PythiaGun::ExecuteTask() {
       if (!printer.empty()){
             std::ofstream sigma_printer;
             sigma_printer.open(printer, std::ios::out | std::ios::app);
-            
+
             sigma_printer << "sigma = " << GetSigmaGen() << " Err =  " << GetSigmaErr() << endl ;
             //sigma_printer.close();
-      
+
 //      JSINFO << BOLDYELLOW << " sigma = " << GetSigmaGen() << " sigma err = " << GetSigmaErr() << " printer = " << printer << " is " << sigma_printer.is_open() ;
     };
 
@@ -275,12 +275,13 @@ void PythiaGun::ExecuteTask() {
     // {
     if (flag_useHybridHad != 1) {
       AddParton(make_shared<Parton>(0, particle.id(), 0, particle.pT(),
-                                    particle.y(), particle.phi(), particle.e(),
-                                    xLoc));
+                                    particle.eta(), particle.phi(),
+                                    particle.e(), xLoc));
     } else {
       auto ptn =
-          make_shared<Parton>(0, particle.id(), 0, particle.pT(), particle.y(),
-                              particle.phi(), particle.e(), xLoc);
+          make_shared<Parton>(0, particle.id(), 0, particle.pT(),
+                              particle.eta(),particle.phi(),
+                              particle.e(), xLoc);
       ptn->set_color(particle.col());
       ptn->set_anti_color(particle.acol());
       ptn->set_max_color(1000 * (np + 1));