[PWGCF] FemtoUniverse: Add quantum weight in MC Truth same-event pairs, fix bug

PWGCF/FemtoUniverse/Core/FemtoUniverseMath.h

@@ -21,7 +21,7 @@
 
 #include "Math/Boost.h"
 #include "Math/Vector4D.h"
 #include "TLorentzVector.h"
 #include "TMath.h"
 
 #include <algorithm>
@@ -177,7 +177,7 @@
   /// \param mass2 Mass of particle 2
   /// \param isiden Identical or non-identical particle pair
   template <typename T>
-  static std::vector<double> newpairfunc(const T& part1, const float mass1, const T& part2, const float mass2, bool isiden)
+  static std::vector<double> newpairfunc(const T& part1, const float mass1, const T& part2, const float mass2, bool isiden, bool isWeight = 0)
   {
     const double e1 = std::sqrt(std::pow(part1.px(), 2) + std::pow(part1.py(), 2) + std::pow(part1.pz(), 2) + std::pow(mass1, 2));
     const double e2 = std::sqrt(std::pow(part2.px(), 2) + std::pow(part2.py(), 2) + std::pow(part2.pz(), 2) + std::pow(mass2, 2));
@@ -186,6 +186,9 @@
     const ROOT::Math::PxPyPzEVector vecpart2(part2.px(), part2.py(), part2.pz(), e2);
     const ROOT::Math::PxPyPzEVector trackSum = vecpart1 + vecpart2;
 
+    const ROOT::Math::PtEtaPhiMVector vecspace_part1(part1.pt(), part1.eta(), part1.phi(), mass1);
+    const ROOT::Math::PtEtaPhiMVector vecspace_part2(part2.pt(), part2.eta(), part2.phi(), mass2);
+
     std::vector<double> vect;
 
     const double tPx = trackSum.px();
@@ -222,6 +225,7 @@
     const double fDKOutLCMS = px1LCMS - px2LCMS;
     const double fDKSideLCMS = py1LCMS - py2LCMS;
     const double fDKLongLCMS = pz1LCMS - pz2LCMS;
+    const double mDE = pE1LCMS - pE2LCMS;
 
     // Boost to PRF
 
@@ -249,6 +253,27 @@
       vect.push_back(fDKSide);
       vect.push_back(fDKLong);
     }
+
+    if (isiden && isWeight) {
+
+      const double x1_lcms = (vecspace_part1.x() * tPx + vecspace_part1.y() * tPy) / tPt;
+      const double y1_lcms = (-vecspace_part1.x() * tPy + vecspace_part1.y() * tPx) / tPt;
+      const double x2_lcms = (vecspace_part2.x() * tPx + vecspace_part2.y() * tPy) / tPt;
+      const double y2_lcms = (-vecspace_part2.x() * tPy + vecspace_part2.y() * tPx) / tPt;
+
+      const double z1_lcms = gamma * (vecspace_part1.z() - beta * vecspace_part1.t());
+      const double t1_lcms = gamma * (vecspace_part1.t() - beta * vecspace_part1.z());
+      const double z2_lcms = gamma * (vecspace_part2.z() - beta * vecspace_part2.t());
+      const double t2_lcms = gamma * (vecspace_part2.t() - beta * vecspace_part2.z());
+
+      const double mRO = (x1_lcms - x2_lcms) / 0.197327;
+      const double mRS = (y1_lcms - y2_lcms) / 0.197327;
+      const double mRL = (z1_lcms - z2_lcms) / 0.197327;
+      const double mDT = (t1_lcms - t2_lcms) / 0.197327;
+
+      const double quantumweight = 1.0 + TMath::Cos(-fDKOutLCMS * mRO - fDKSideLCMS * mRS - fDKLongLCMS * mRL + mDE * mDT);
+      vect.push_back(quantumweight);
+    }
     return vect;
   }
 };