Fix: fix Q3 computation

Author: ariedel-cern

PWGCF/Femto/Core/tripletHistManager.h

@@ -27,6 +27,7 @@
 #include <Math/GenVector/Boost.h>
 #include <Math/Vector4D.h>
 
+#include <algorithm>
 #include <array>
 #include <cmath>
 #include <map>
@@ -273,9 +274,9 @@ class TripletHistManager
     mParticle3 = ROOT::Math::PtEtaPhiMVector(mAbsCharge3 * particle3.pt(), particle3.eta(), particle3.phi(), mPdgMass3);
 
     // set mT
-    mMt = getMt(mParticle1, mPdgMass1, mParticle2, mPdgMass2, mParticle3, mPdgMass3);
+    mMt = getMt(mParticle1, mParticle2, mParticle3);
     // set Q3
-    mQ3 = getQ3(mParticle1, mPdgMass1, mParticle2, mPdgMass2, mParticle3, mPdgMass3);
+    mQ3 = getQ3(mParticle1, mParticle2, mParticle3);
 
     // if one of the particles has a mass getter, we cache the value for the filling later
     if constexpr (modes::hasMass(particleType1)) {
@@ -323,9 +324,9 @@ class TripletHistManager
     mTrueParticle3 = ROOT::Math::PtEtaPhiMVector(mAbsCharge3 * mcParticle3.pt(), mcParticle3.eta(), mcParticle3.phi(), mPdgMass3);
 
     // set true mT
-    mTrueMt = getMt(mTrueParticle1, mPdgMass1, mTrueParticle2, mPdgMass2, mTrueParticle3, mPdgMass3);
+    mTrueMt = getMt(mTrueParticle1, mTrueParticle2, mTrueParticle3);
     // set true Q3
-    mTrueQ3 = getQ3(mTrueParticle1, mPdgMass1, mTrueParticle2, mPdgMass2, mTrueParticle3, mPdgMass3);
+    mTrueQ3 = getQ3(mTrueParticle1, mTrueParticle2, mTrueParticle3);
   }
 
   template <typename T1, typename T2, typename T3, typename T4, typename T5, typename T6>
@@ -382,26 +383,13 @@ class TripletHistManager
   float getQ3() const { return mQ3; }
 
  private:
-  template <typename T>
-  T getqij(const T& vecparti, const T& vecpartj)
+  float getQ3(ROOT::Math::PtEtaPhiMVector const& part1, ROOT::Math::PtEtaPhiMVector const& part2, ROOT::Math::PtEtaPhiMVector const& part3)
   {
-    T trackSum = vecparti + vecpartj;
-    T trackDifference = vecparti - vecpartj;
-    float scaling = trackDifference.Dot(trackSum) / trackSum.Dot(trackSum);
-    return trackDifference - scaling * trackSum;
-  }
-
-  template <typename T>
-  float getQ3(const T& part1, const float mass1, const T& part2, const float mass2, const T& part3, const float mass3)
-  {
-    ROOT::Math::PtEtaPhiMVector vecpart1(part1.pt(), part1.eta(), part1.phi(), mass1);
-    ROOT::Math::PtEtaPhiMVector vecpart2(part2.pt(), part2.eta(), part2.phi(), mass2);
-    ROOT::Math::PtEtaPhiMVector vecpart3(part3.pt(), part3.eta(), part3.phi(), mass3);
-    auto q12 = getqij(vecpart1, vecpart2);
-    auto q23 = getqij(vecpart2, vecpart3);
-    auto q31 = getqij(vecpart3, vecpart1);
-    float q = q12.M2() + q23.M2() + q31.M2();
-    return std::sqrt(-q);
+    double q12 = -(part1 - part2).M2();
+    double q23 = -(part2 - part3).M2();
+    double q31 = -(part3 - part1).M2();
+    double q = q12 + q23 + q31;
+    return static_cast<float>(std::sqrt(std::max(0., q))); // protection if rounding error give a value below 0
   }
 
   void initAnalysis(std::map<TripletHist, std::vector<o2::framework::AxisSpec>> const& Specs)
@@ -497,22 +485,19 @@ class TripletHistManager
     }
   }
 
-  float getMt(ROOT::Math::PtEtaPhiMVector const& part1,
-              float pdgMass1,
-              ROOT::Math::PtEtaPhiMVector const& part2,
-              float pdgMass2)
+  double getMt(ROOT::Math::PtEtaPhiMVector const& part1, ROOT::Math::PtEtaPhiMVector const& part2)
   {
     auto sum = part1 + part2;
-    float mt = 0;
-    float averageMass = 0;
-    float reducedMass = 0;
+    double mt = 0.;
+    double averageMass = 0.;
+    double reducedMass = 0.;
     switch (mMtType) {
       case modes::TransverseMassType::kAveragePdgMass:
-        averageMass = 0.5 * (pdgMass1 + pdgMass2);
+        averageMass = 0.5 * (part1.M() + part2.M());
         mt = std::hypot(0.5 * sum.Pt(), averageMass);
         break;
       case modes::TransverseMassType::kReducedPdgMass:
-        reducedMass = 2.f * (pdgMass1 * pdgMass2 / (pdgMass1 + pdgMass2));
+        reducedMass = 2. * (part1.M() * part2.M() / (part1.M() + part2.M()));
         mt = std::hypot(0.5 * sum.Pt(), reducedMass);
         break;
       case modes::TransverseMassType::kMt4Vector:
@@ -524,20 +509,15 @@ class TripletHistManager
     return mt;
   }
 
-  float getMt(ROOT::Math::PtEtaPhiMVector const& part1,
-              float mass1,
-              ROOT::Math::PtEtaPhiMVector const& part2,
-              float mass2,
-              ROOT::Math::PtEtaPhiMVector const& part3,
-              float mass3)
+  float getMt(ROOT::Math::PtEtaPhiMVector const& part1, ROOT::Math::PtEtaPhiMVector const& part2, ROOT::Math::PtEtaPhiMVector const& part3)
   {
-    return (getMt(part1, mass1, part2, mass2) + getMt(part2, mass2, part3, mass3) + getMt(part1, mass1, part3, mass3)) / 3.;
+    return static_cast<float>((getMt(part1, part2) + getMt(part2, part3) + getMt(part1, part3)) / 3.);
   }
 
   o2::framework::HistogramRegistry* mHistogramRegistry = nullptr;
-  float mPdgMass1 = 0.f;
-  float mPdgMass2 = 0.f;
-  float mPdgMass3 = 0.f;
+  double mPdgMass1 = 0.f;
+  double mPdgMass2 = 0.f;
+  double mPdgMass3 = 0.f;
 
   modes::TransverseMassType mMtType = modes::TransverseMassType::kAveragePdgMass;