Feat: use double for intermediate computations

Author: ariedel-cern

PWGCF/Femto/Core/closePairRejection.h

@@ -214,8 +214,8 @@ class CloseTrackRejection
     mDeta = track1.eta() - track2.eta();
 
     for (size_t i = 0; i < TpcRadii.size(); i++) {
-      auto phistar1 = utils::dphistar(mMagField, TpcRadii[i], mChargeAbsTrack1 * track1.signedPt(), track1.phi());
-      auto phistar2 = utils::dphistar(mMagField, TpcRadii[i], mChargeAbsTrack2 * track2.signedPt(), track2.phi());
+      auto phistar1 = phistar(mMagField, TpcRadii[i], mChargeAbsTrack1 * track1.signedPt(), track1.phi());
+      auto phistar2 = phistar(mMagField, TpcRadii[i], mChargeAbsTrack2 * track2.signedPt(), track2.phi());
       if (phistar1 && phistar2) {
         mDphistar.at(i) = RecoDecay::constrainAngle(phistar1.value() - phistar2.value(), -o2::constants::math::PI); // constrain angular difference between -pi and pi
         mDphistarMask.at(i) = true;
@@ -309,6 +309,16 @@ class CloseTrackRejection
   bool isActivated() const { return mIsActivated; }
 
  private:
+  std::optional<float> phistar(float magfield, float radius, float signedPt, float phi)
+  {
+    double arg = 0.3 * (0.1 * magfield) * (0.01 * radius) / (2. * signedPt);
+    if (std::fabs(arg) <= 1.) {
+      double phistar = phi - std::asin(arg);
+      return static_cast<float>(RecoDecay::constrainAngle(phistar));
+    }
+    return std::nullopt;
+  }
+
   o2::framework::HistogramRegistry* mHistogramRegistry = nullptr;
   bool mPlotAllRadii = false;
   bool mPlotAverage = false;

PWGCF/Femto/Core/collisionBuilder.h

@@ -225,7 +225,7 @@ class CollisionSelection : public BaseSelection<float, o2::aod::femtodatatypes::
   template <typename T>
   void setSphericity(T tracks)
   {
-    mSphericity = utils::sphericity(tracks);
+    mSphericity = computeSphericity(tracks);
   }
 
   float getSphericity() const { return mSphericity; }
@@ -312,6 +312,39 @@ class CollisionSelection : public BaseSelection<float, o2::aod::femtodatatypes::
   };
 
  protected:
+  template <typename T>
+  float computeSphericity(T const& tracks)
+  {
+    int minNumberTracks = 2;
+    double maxSphericity = 2.f;
+    if (tracks.size() <= minNumberTracks) {
+      return maxSphericity;
+    }
+    // Initialize the transverse momentum tensor components
+    double sxx = 0.;
+    double syy = 0.;
+    double sxy = 0.;
+    double sumPt = 0.;
+    // Loop over the tracks to compute the tensor components
+    for (const auto& track : tracks) {
+      sxx += (track.px() * track.px()) / track.pt();
+      syy += (track.py() * track.py()) / track.pt();
+      sxy += (track.px() * track.py()) / track.pt();
+      sumPt += track.pt();
+    }
+    sxx /= sumPt;
+    syy /= sumPt;
+    sxy /= sumPt;
+    // Compute the eigenvalues (real values)
+    double lambda1 = ((sxx + syy) + std::sqrt((sxx + syy) * (sxx + syy) - 4. * (sxx * syy - sxy * sxy))) / 2.;
+    double lambda2 = ((sxx + syy) - std::sqrt((sxx + syy) * (sxx + syy) - 4. * (sxx * syy - sxy * sxy))) / 2.;
+    if (lambda1 <= 0. || lambda2 <= 0.) {
+      return maxSphericity;
+    }
+    // Compute sphericity
+    return static_cast<float>(2. * lambda2 / (lambda1 + lambda2));
+  }
+
   // filter cuts
   float mVtxZMin = -12.f;
   float mVtxZMax = -12.f;

PWGCF/Femto/Core/femtoUtils.h

@@ -16,11 +16,8 @@
 #ifndef PWGCF_FEMTO_CORE_FEMTOUTILS_H_
 #define PWGCF_FEMTO_CORE_FEMTOUTILS_H_
 
-#include "RecoDecay.h"
-
 #include "Common/Core/TableHelper.h"
 
-#include "CommonConstants/MathConstants.h"
 #include "CommonConstants/PhysicsConstants.h"
 #include "Framework/InitContext.h"
 
@@ -63,55 +60,16 @@ float itsSignal(T const& track)
   auto clSizeLayer6 = (clsizeflag >> (6 * 4)) & 0xf;
   int numLayers = 7;
   int sumClusterSizes = clSizeLayer0 + clSizeLayer1 + clSizeLayer2 + clSizeLayer3 + clSizeLayer4 + clSizeLayer5 + clSizeLayer6;
-  float cosLamnda = 1. / std::cosh(track.eta());
-  return (static_cast<float>(sumClusterSizes) / numLayers) * cosLamnda;
+  double cosLamnda = 1. / std::cosh(track.eta());
+  double signal = (static_cast<double>(sumClusterSizes) / numLayers) * cosLamnda;
+  return static_cast<float>(signal);
 };
 
-template <typename T>
-float sphericity(T const& tracks)
-{
-
-  int minNumberTracks = 2;
-  float maxSphericity = 2.f;
-
-  if (tracks.size() <= minNumberTracks) {
-    return maxSphericity;
-  }
-
-  // Initialize the transverse momentum tensor components
-  float sxx = 0.f;
-  float syy = 0.f;
-  float sxy = 0.f;
-  float sumPt = 0.f;
-
-  // Loop over the tracks to compute the tensor components
-  for (const auto& track : tracks) {
-    sxx += (track.px() * track.px()) / track.pt();
-    syy += (track.py() * track.py()) / track.pt();
-    sxy += (track.px() * track.py()) / track.pt();
-    sumPt += track.pt();
-  }
-  sxx /= sumPt;
-  syy /= sumPt;
-  sxy /= sumPt;
-
-  // Compute the eigenvalues (real values)
-  float lambda1 = ((sxx + syy) + std::sqrt((sxx + syy) * (sxx + syy) - 4 * (sxx * syy - sxy * sxy))) / 2;
-  float lambda2 = ((sxx + syy) - std::sqrt((sxx + syy) * (sxx + syy) - 4 * (sxx * syy - sxy * sxy))) / 2;
-
-  if (lambda1 <= 0.f || lambda2 <= 0.f) {
-    return maxSphericity;
-  }
-
-  // Compute sphericity
-  return 2.f * lambda2 / (lambda1 + lambda2);
-}
-
-inline float getMass(int pdgCode)
+inline double getMass(int pdgCode)
 {
   // use this function instead of TDatabasePDG to return masses defined in the PhysicsConstants.h header
   // this approach saves a lot of memory and important partilces like deuteron are missing in TDatabasePDG anyway
-  float mass = 0.f;
+  double mass = 0.f;
   // add new particles if necessary here
   switch (std::abs(pdgCode)) {
     case kPiPlus:
@@ -175,15 +133,6 @@ float qn(T const& col)
   return qn;
 }
 
-inline std::optional<float> dphistar(float magfield, float radius, float signedPt, float phi)
-{
-  float arg = 0.3f * (0.1f * magfield) * (0.01 * radius) / (2.f * signedPt);
-  if (std::fabs(arg) <= 1.f) {
-    return RecoDecay::constrainAngle(phi - std::asin(arg));
-  }
-  return std::nullopt;
-}
-
 inline bool enableTable(const char* tableName, int userSetting, o2::framework::InitContext& initContext)
 {
   if (userSetting == 1) {

PWGCF/Femto/Core/pairHistManager.h

@@ -596,14 +596,14 @@ class PairHistManager
 
   float getKt(ROOT::Math::PtEtaPhiMVector const& part1, ROOT::Math::PtEtaPhiMVector const& part2)
   {
-    auto sum = part1 + part2;
-    return 0.5 * sum.Pt();
+    double kt = 0.5 * (part1 + part2).Pt();
+    return static_cast<float>(kt);
   }
 
   float getMt(ROOT::Math::PtEtaPhiMVector const& part1, ROOT::Math::PtEtaPhiMVector const& part2)
   {
     auto sum = part1 + part2;
-    float mt = 0;
+    double mt = 0;
     switch (mMtType) {
       case modes::TransverseMassType::kAveragePdgMass:
         mt = std::hypot(0.5 * sum.Pt(), mAverageMass);
@@ -617,7 +617,7 @@ class PairHistManager
       default:
         LOG(fatal) << "Invalid transverse mass type, breaking...";
     }
-    return mt;
+    return static_cast<float>(mt);
   }
 
   float getKstar(ROOT::Math::PtEtaPhiMVector const& part1, ROOT::Math::PtEtaPhiMVector const& part2)
@@ -630,16 +630,16 @@ class PairHistManager
     // get lorentz boost into pair rest frame
     ROOT::Math::Boost boostPrf(sum.BoostToCM());
     // boost particle 1 into pair rest frame and calculate its momentum, which has the same value as k*
-    return boostPrf(particle1Prf).P();
+    return static_cast<float>(boostPrf(particle1Prf).P());
   }
 
   o2::framework::HistogramRegistry* mHistogramRegistry = nullptr;
   float mPdgMass1 = 0.f;
   float mPdgMass2 = 0.f;
 
   modes::TransverseMassType mMtType = modes::TransverseMassType::kAveragePdgMass;
-  float mAverageMass = 0.f;
-  float mReducedMass = 0.f;
+  double mAverageMass = 0.f;
+  double mReducedMass = 0.f;
 
   int mAbsCharge1 = 1;
   int mAbsCharge2 = 1;