[PWGLF] added second option for mass calculation (#11459)

Co-authored-by: ALICE Action Bot <alibuild@cern.ch>

Author: mherzer28

PWGLF/TableProducer/Nuspex/trHeAnalysis.cxx

@@ -35,6 +35,7 @@
 #include "PWGLF/DataModel/LFParticleIdentification.h"
 #include "ReconstructionDataFormats/PID.h"
 #include "ReconstructionDataFormats/Track.h"
+#include "PWGLF/DataModel/pidTOFGeneric.h"
 #include <TF1.h>
 
 namespace o2::aod
@@ -113,7 +114,7 @@ static const std::vector<int> particleCharge{1, 2};
 static const std::vector<float> particleChargeFactor{2.3, 2.55};
 static const std::vector<std::string> betheBlochParNames{
   "p0", "p1", "p2", "p3", "p4", "resolution"};
-constexpr float betheBlochDefault[nParticles][nBetheParams]{
+constexpr float BetheBlochDefault[nParticles][nBetheParams]{
   {0.248753, 3.58634, 0.0167065, 2.29194, 0.774344,
    0.07}, // triton
   {0.0274556, 18.3054, 3.99987e-05, 3.17219, 11.1775,
@@ -126,7 +127,8 @@ using TracksFull =
   soa::Join<aod::TracksIU, aod::TracksExtra, aod::TracksCovIU,
             o2::aod::TracksDCA, aod::pidTOFmass, aod::pidTOFbeta,
             aod::pidTPCLfFullTr, aod::pidTPCLfFullHe,
-            aod::TOFSignal, aod::TrackSelectionExtension>;
+            aod::TOFSignal, aod::TrackSelectionExtension,
+            o2::aod::EvTimeTOFFT0ForTrack>;
 
 class Particle
 {
@@ -186,7 +188,7 @@ struct TrHeAnalysis {
   } evselOptions;
 
   Configurable<bool> cfgTPCPidMethod{"cfgTPCPidMethod", false, "Using own or built in bethe parametrization"}; // false for built in
-
+  Configurable<int> cfgMassMethod{"cfgMassMethod", 0, "0: Using built in 1: mass calculated with beta 2: mass calculated with the event time"};
   // Set the multiplity event limits
   Configurable<float> cfgLowMultCut{"cfgLowMultCut", 0.0f, "Accepted multiplicity percentage lower limit"};
   Configurable<float> cfgHighMultCut{"cfgHighMultCut", 100.0f, "Accepted multiplicity percentage higher limit"};
@@ -226,7 +228,7 @@ struct TrHeAnalysis {
     Configurable<float> nsigmaTPCTr{"nsigmaTPCTr", 5.f, "Value of the Nsigma TPC cut for tritons"};
     Configurable<float> nsigmaTPCHe{"nsigmaTPCHe", 5.f, "Value of the Nsigma TPC cut for helium-3"};
   } nsigmaTPCvar;
-  Configurable<LabeledArray<float>> cfgBetheBlochParams{"cfgBetheBlochParams", {betheBlochDefault[0], nParticles, nBetheParams, particleNames, betheBlochParNames}, "TPC Bethe-Bloch parameterisation for light nuclei"};
+  Configurable<LabeledArray<float>> cfgBetheBlochParams{"cfgBetheBlochParams", {BetheBlochDefault[0], nParticles, nBetheParams, particleNames, betheBlochParNames}, "TPC Bethe-Bloch parameterisation for light nuclei"};
 
   void init(o2::framework::InitContext&)
   {
@@ -663,12 +665,14 @@ struct TrHeAnalysis {
     constexpr int NNumLayers = 8;
     constexpr int NBitsPerLayer = 4;
     constexpr int NBitMask = (1 << NBitsPerLayer) - 1;
+
     int sum = 0, n = 0;
     for (int i = 0; i < NNumLayers; i++) {
       int clsSize = (track.itsClusterSizes() >> (NBitsPerLayer * i)) & NBitMask;
       sum += clsSize;
-      if (clsSize)
+      if (clsSize) {
         n++;
+      }
     }
     return n > 0 ? static_cast<float>(sum) / n : 0.f;
   }
@@ -681,13 +685,34 @@ struct TrHeAnalysis {
     return hePID ? track.tpcInnerParam() / 2 : track.tpcInnerParam();
   }
   template <class T>
-  float getMass(T const& track)
+  float getMass(const T& track)
   {
-    const float beta = track.beta();
-    const float rigidity = getRigidity(track);
-    float gamma = 1 / std::sqrt(1 - beta * beta);
-    float mass = (rigidity / std::sqrt(gamma * gamma - 1));
-    return mass;
+    if (cfgMassMethod == 0) {
+      return track.mass();
+    }
+    if (cfgMassMethod == 1) {
+      const float beta = track.beta();
+      const float rigidity = getRigidity(track);
+      float gamma = 1 / std::sqrt(1 - beta * beta);
+      float mass = (rigidity / std::sqrt(gamma * gamma - 1.f));
+      return mass;
+    }
+    if (cfgMassMethod == 2) {
+      const float rigidity = getRigidity(track);
+      float tofStartTime = track.evTimeForTrack();
+      float tofTime = track.tofSignal();
+      constexpr float CInCmPs = 2.99792458e-2f;
+      float length = track.length();
+      float time = tofTime - tofStartTime;
+      if (time > 0.f && length > 0.f) {
+        float beta = length / (CInCmPs * time);
+        float gamma = 1 / std::sqrt(1 - beta * beta);
+        float mass = rigidity / std::sqrt(gamma * gamma - 1.f);
+        return mass;
+      }
+      return -1.f;
+    }
+    return -1.f;
   }
 };