Adding also estimations of the number of peaks, and a random variables to be used for random event rejection

Author: Ionut Cristian Arsene

PWGDQ/Core/VarManager.cxx

@@ -359,25 +359,24 @@ double VarManager::ComputePIDcalibration(int species, double nSigmaValue)
 }
 
 //__________________________________________________________________
-std::tuple<double, double, double, double, double> VarManager::BimodalityCoefficientUnbinned(const std::vector<double>& data)
+std::tuple<float, float, float, float, float> VarManager::BimodalityCoefficientUnbinned(const std::vector<float>& data)
 {
   // Bimodality coefficient = (skewness^2 + 1) / kurtosis
   // return a tuple including the coefficient, mean, RMS, skewness, and kurtosis
   size_t n = data.size();
   if (n < 3) {
     return std::make_tuple(-1.0, -1.0, -1.0, -1.0, -1.0);
   }
-  double mean = std::accumulate(data.begin(), data.end(), 0.0) / n;
+  float mean = std::accumulate(data.begin(), data.end(), 0.0) / n;
 
-  double m2 = 0.0, m3 = 0.0, m4 = 0.0;
-  for (double x : data) {
-    double diff = x - mean;
-    double diff2 = diff * diff;
-    double diff3 = diff2 * diff;
-    double diff4 = diff3 * diff;
+  float m2 = 0.0, m3 = 0.0, m4 = 0.0;
+  float diff, diff2;
+  for (float x : data) {
+    diff = x - mean;
+    diff2 = diff * diff;
     m2 += diff2;
-    m3 += diff3;
-    m4 += diff4;
+    m3 += diff2 * diff;
+    m4 += diff2 * diff2;
   }
 
   m2 /= n;
@@ -388,21 +387,23 @@ std::tuple<double, double, double, double, double> VarManager::BimodalityCoeffic
     return std::make_tuple(-1.0, -1.0, -1.0, -1.0, -1.0);
   }
 
-  double stddev = std::sqrt(m2);
-  double skewness = m3 / (stddev * stddev * stddev);
-  double kurtosis = m4 / (m2 * m2);
+  float stddev = std::sqrt(m2);
+  float skewness = m3 / (stddev * stddev * stddev);
+  float kurtosis = m4 / (m2 * m2);
 
   return std::make_tuple((skewness * skewness + 1.0) / kurtosis, mean, stddev, skewness, kurtosis);
 }
 
-std::tuple<double, double, double, double, double> VarManager::BimodalityCoefficient(const std::vector<double>& data, float binWidth, int trim, float min, float max)
+std::tuple<float, float, float, float, float, int> VarManager::BimodalityCoefficientAndNPeaks(const std::vector<float>& data, float binWidth, int trim, float min, float max)
 {
   // Bimodality coefficient = (skewness^2 + 1) / kurtosis
   // return a tuple including the coefficient, mean, RMS, skewness, and kurtosis
 
   // if the binWidth is < 0, use the unbinned calculation
   if (binWidth < 0) {
-    return BimodalityCoefficientUnbinned(data);
+    // get the tuple from the unbinned calculation
+    auto [bc, mean, stddev, skewness, kurtosis] = BimodalityCoefficientUnbinned(data);
+    return std::make_tuple(bc, mean, stddev, skewness, kurtosis, -1);
   }
 
   // bin the data and put it in a vector
@@ -428,46 +429,62 @@ std::tuple<double, double, double, double, double> VarManager::BimodalityCoeffic
     }
   }
 
+  // count the number of peaks
+  int nPeaks = 0;
+  bool inPeak = false;
+  for (int i = 0; i < nBins; ++i) {
+    if (counts[i] > 0) {
+      if (!inPeak) {
+        inPeak = true;
+        nPeaks++;
+      }
+    } else {
+      inPeak = false;
+    }
+  }
+
   // first compute the mean
-  double mean = 0.0;
-  double totalCounts = 0.0;
+  float mean = 0.0;
+  float totalCounts = 0.0;
   for (int i = 0; i < nBins; ++i) {
-    double binCenter = min + (i + 0.5) * binWidth;
+    float binCenter = min + (i + 0.5) * binWidth;
     mean += counts[i] * binCenter;
     totalCounts += counts[i];
   }
 
   if (totalCounts == 0) {
-    return std::make_tuple(-1.0, -1.0, -1.0, -1.0, -1.0);
+    return std::make_tuple(-1.0, -1.0, -1.0, -1.0, -1.0, -1);
   }
   mean /= totalCounts;
 
   // then compute the second, third, and fourth central moments
-  double m2 = 0.0, m3 = 0.0, m4 = 0.0;
+  float m2 = 0.0, m3 = 0.0, m4 = 0.0;
+  float diff, diff2, binCenter;
   for (int i = 0; i < nBins; ++i) {
-    double binCenter = min + (i + 0.5) * binWidth;
-    double diff = binCenter - mean;
-    double diff2 = diff * diff;
-    double diff3 = diff2 * diff;
-    double diff4 = diff3 * diff;
+    if (counts[i] == 0) {
+      continue; // skip empty bins
+    }
+    binCenter = min + (i + 0.5) * binWidth;
+    diff = binCenter - mean;
+    diff2 = diff * diff;
     m2 += counts[i] * diff2;
-    m3 += counts[i] * diff3;
-    m4 += counts[i] * diff4;
+    m3 += counts[i] * diff2 * diff;
+    m4 += counts[i] * diff2 * diff2;
   }
 
   m2 /= totalCounts;
   m3 /= totalCounts;
   m4 /= totalCounts;
 
   if (m2 == 0.0) {
-    return std::make_tuple(-1.0, -1.0, -1.0, -1.0, -1.0);
+    return std::make_tuple(-1.0, -1.0, -1.0, -1.0, -1.0, -1);
   }
 
-  double stddev = std::sqrt(m2);
-  double skewness = m3 / (stddev * stddev * stddev);
-  double kurtosis = m4 / (m2 * m2); // Pearson's kurtosis, not excess
+  float stddev = std::sqrt(m2);
+  float skewness = m3 / (stddev * stddev * stddev);
+  float kurtosis = m4 / (m2 * m2); // Pearson's kurtosis, not excess
 
-  return std::make_tuple((skewness * skewness + 1.0) / kurtosis, mean, stddev, skewness, kurtosis);
+  return std::make_tuple((skewness * skewness + 1.0) / kurtosis, mean, stddev, skewness, kurtosis, nPeaks);
 }
 
 //__________________________________________________________________
@@ -501,6 +518,8 @@ void VarManager::SetDefaultVarNames()
   fgVariableUnits[kTimeFromSOR] = "min.";
   fgVariableNames[kBCOrbit] = "Bunch crossing";
   fgVariableUnits[kBCOrbit] = "";
+  fgVariableNames[kCollisionRandom] = "Random number (collision-level)";
+  fgVariableUnits[kCollisionRandom] = "";
   fgVariableNames[kIsPhysicsSelection] = "Physics selection";
   fgVariableUnits[kIsPhysicsSelection] = "";
   fgVariableNames[kVtxX] = "Vtx X ";
@@ -733,6 +752,14 @@ void VarManager::SetDefaultVarNames()
   fgVariableUnits[kDCAzFracAbove5mm] = "";
   fgVariableNames[kDCAzFracAbove10mm] = "Fraction of tracks with |DCAz| > 10 mm";
   fgVariableUnits[kDCAzFracAbove10mm] = "";
+  fgVariableNames[kDCAzNPeaks] = "Number of peaks in DCAz distribution";
+  fgVariableUnits[kDCAzNPeaks] = "";
+  fgVariableNames[kDCAzNPeaksTrimmed1] = "Number of peaks in binned DCAz distribution (trimmed 1)";
+  fgVariableUnits[kDCAzNPeaksTrimmed1] = "";
+  fgVariableNames[kDCAzNPeaksTrimmed2] = "Number of peaks in binned DCAz distribution (trimmed 2)";
+  fgVariableUnits[kDCAzNPeaksTrimmed2] = "";
+  fgVariableNames[kDCAzNPeaksTrimmed3] = "Number of peaks in binned DCAz distribution (trimmed 3)";
+  fgVariableUnits[kDCAzNPeaksTrimmed3] = "";
   fgVariableNames[kPt] = "p_{T}";
   fgVariableUnits[kPt] = "GeV/c";
   fgVariableNames[kPt1] = "p_{T1}";
@@ -1709,6 +1736,7 @@ void VarManager::SetDefaultVarNames()
   fgVarNamesMap["kCollisionTimeRes"] = kCollisionTimeRes;
   fgVarNamesMap["kBC"] = kBC;
   fgVarNamesMap["kBCOrbit"] = kBCOrbit;
+  fgVarNamesMap["kCollisionRandom"] = kCollisionRandom;
   fgVarNamesMap["kIsPhysicsSelection"] = kIsPhysicsSelection;
   fgVarNamesMap["kIsTVXTriggered"] = kIsTVXTriggered;
   fgVarNamesMap["kIsNoTFBorder"] = kIsNoTFBorder;
@@ -1816,6 +1844,10 @@ void VarManager::SetDefaultVarNames()
   fgVarNamesMap["kDCAzFracAbove2mm"] = kDCAzFracAbove2mm;
   fgVarNamesMap["kDCAzFracAbove5mm"] = kDCAzFracAbove5mm;
   fgVarNamesMap["kDCAzFracAbove10mm"] = kDCAzFracAbove10mm;
+  fgVarNamesMap["kDCAzNPeaks"] = kDCAzNPeaks;
+  fgVarNamesMap["kDCAzNPeaksTrimmed1"] = kDCAzNPeaksTrimmed1;
+  fgVarNamesMap["kDCAzNPeaksTrimmed2"] = kDCAzNPeaksTrimmed2;
+  fgVarNamesMap["kDCAzNPeaksTrimmed3"] = kDCAzNPeaksTrimmed3;
   fgVarNamesMap["kMCEventGeneratorId"] = kMCEventGeneratorId;
   fgVarNamesMap["kMCEventSubGeneratorId"] = kMCEventSubGeneratorId;
   fgVarNamesMap["kMCVtxX"] = kMCVtxX;

PWGDQ/Core/VarManager.h

@@ -206,6 +206,7 @@ class VarManager : public TObject
     kCollisionTimeRes,
     kBC,
     kBCOrbit,
+    kCollisionRandom,            // random number generated per collision (if required, can be used to perform random selections at the collision level)
     kIsPhysicsSelection,
     kIsTVXTriggered,             // Is trigger TVX
     kIsNoTFBorder,               // No time frame border
@@ -315,6 +316,10 @@ class VarManager : public TObject
     kDCAzFracAbove2mm,
     kDCAzFracAbove5mm,
     kDCAzFracAbove10mm,
+    kDCAzNPeaks,
+    kDCAzNPeaksTrimmed1,
+    kDCAzNPeaksTrimmed2,
+    kDCAzNPeaksTrimmed3,
     kMCEventGeneratorId,
     kMCEventSubGeneratorId,
     kMCVtxX,
@@ -1278,8 +1283,8 @@ class VarManager : public TObject
   template <typename T, typename T1>
   static o2::dataformats::VertexBase RecalculatePrimaryVertex(T const& track0, T const& track1, const T1& collision);
 
-  static std::tuple<double, double, double, double, double> BimodalityCoefficientUnbinned(const std::vector<double>& data);
-  static std::tuple<double, double, double, double, double> BimodalityCoefficient(const std::vector<double>& data, float binWidth, int trim = 0, float min = -15.0, float max = 15.0);
+  static std::tuple<float, float, float, float, float> BimodalityCoefficientUnbinned(const std::vector<float>& data);
+  static std::tuple<float, float, float, float, float, int> BimodalityCoefficientAndNPeaks(const std::vector<float>& data, float binWidth, int trim = 0, float min = -15.0, float max = 15.0);
 
   template <typename T, typename C>
   static o2::track::TrackParCovFwd FwdToTrackPar(const T& track, const C& cov);
@@ -1867,6 +1872,10 @@ void VarManager::FillEvent(T const& event, float* values)
     values[kTimestamp] = event.timestamp();
   }
 
+  if (fgUsedVars[kCollisionRandom]) {
+    values[kCollisionRandom] = gRandom->Rndm();
+  }
+
   if constexpr ((fillMap & Collision) > 0) {
     // TODO: trigger info from the event selection requires a separate flag
     //       so that it can be switched off independently of the rest of Collision variables (e.g. if event selection is not available)
@@ -2384,58 +2393,71 @@ void VarManager::FillEventTracks(T const& tracks, float* values)
   }
 
   // compute the unbinned bimodality coefficient and related statistics
-  auto [bimodality, mean, stddev, skewness, kurtosis] = BimodalityCoefficient(dcazValues, -1.0);
+  auto [bimodality, mean, stddev, skewness, kurtosis, nPeaks] = BimodalityCoefficient(dcazValues, -1.0);
   if (stddev > -1.0) {
     values[kDCAzBimodalityCoefficient] = bimodality;
     values[kDCAzMean] = mean;
     values[kDCAzRMS] = stddev;
     values[kDCAzSkewness] = skewness;
     values[kDCAzKurtosis] = kurtosis;
+    values[kDCAzNPeaks] = -9999.0;
   } else {
     values[kDCAzBimodalityCoefficient] = -9999.0;
     values[kDCAzMean] = -9999.0;
     values[kDCAzRMS] = -9999.0;
     values[kDCAzSkewness] = -9999.0;
     values[kDCAzKurtosis] = -9999.0;
-  }
+    values[kDCAzNPeaks] = -9999.0;
   // compute the binned bimodality coefficient and related statistics with a bin width of 50um
-  auto [bimodalityBin, meanBin, stddevBin, skewnessBin, kurtosisBin] = BimodalityCoefficient(dcazValues, 0.005);
+  auto [bimodalityBin, meanBin, stddevBin, skewnessBin, kurtosisBin, nPeaksBinned] = BimodalityCoefficient(dcazValues, 0.005);
   if (stddevBin > -1.0) {
     values[kDCAzBimodalityCoefficientBinned] = bimodalityBin;
+    values[kDCAzNPeaks] = nPeaksBinned;
   } else {
     values[kDCAzBimodalityCoefficientBinned] = -9999.0;
+    values[kDCAzNPeaks] = -9999.0;
   }
+  cout << "Bimodality coefficient: " << bimodality << ", mean: " << mean << ", stddev: " << stddev << ", skewness: " << skewness << ", kurtosis: " << kurtosis << ", nPeaks: " << nPeaks << endl;
   // compute the binned bimodality coefficient and related statistics with a bin width of 50um and trimming of 1, 2, 3 entries per bin
-  auto [bimodalityBinTrimm1, meanBinTrimm1, stddevBinTrimm1, skewnessBinTrimm1, kurtosisBinTrimm1] = BimodalityCoefficient(dcazValues, 0.005, 2);
+  auto [bimodalityBinTrimm1, meanBinTrimm1, stddevBinTrimm1, skewnessBinTrimm1, kurtosisBinTrimm1, nPeaksTrimm1] = BimodalityCoefficient(dcazValues, 0.005, 2);
   if (stddevBinTrimm1 > -1.0) {
     values[kDCAzBimodalityCoefficientBinnedTrimmed1] = bimodalityBinTrimm1;
     values[kDCAzMeanBinnedTrimmed1] = meanBinTrimm1;
     values[kDCAzRMSBinnedTrimmed1] = stddevBinTrimm1;
+    values[kDCAzNPeaksTrimmed1] = nPeaksTrimm1;
   } else {
     values[kDCAzBimodalityCoefficientBinnedTrimmed1] = -9999.0;
     values[kDCAzMeanBinnedTrimmed1] = -9999.0;
     values[kDCAzRMSBinnedTrimmed1] = -9999.0;
+    values[kDCAzNPeaksTrimmed1] = -9999.0;
   }
-  auto [bimodalityBinTrimm2, meanBinTrimm2, stddevBinTrimm2, skewnessBinTrimm2, kurtosisBinTrimm2] = BimodalityCoefficient(dcazValues, 0.005, 3);
+  cout << "Bimodality coefficient (trimmed 1): " << bimodalityBinTrimm1 << ", mean: " << meanBinTrimm1 << ", stddev: " << stddevBinTrimm1 << ", skewness: " << skewnessBinTrimm1 << ", kurtosis: " << kurtosisBinTrimm1 << ", nPeaks: " << nPeaksTrimm1 << endl;
+  auto [bimodalityBinTrimm2, meanBinTrimm2, stddevBinTrimm2, skewnessBinTrimm2, kurtosisBinTrimm2, nPeaksTrimm2] = BimodalityCoefficient(dcazValues, 0.005, 3);
   if (stddevBinTrimm2 > -1.0) {
     values[kDCAzBimodalityCoefficientBinnedTrimmed2] = bimodalityBinTrimm2;
     values[kDCAzMeanBinnedTrimmed2] = meanBinTrimm2;
     values[kDCAzRMSBinnedTrimmed2] = stddevBinTrimm2;
+    values[kDCAzNPeaksTrimmed2] = nPeaksTrimm2;
   } else {
     values[kDCAzBimodalityCoefficientBinnedTrimmed2] = -9999.0;
     values[kDCAzMeanBinnedTrimmed2] = -9999.0;
     values[kDCAzRMSBinnedTrimmed2] = -9999.0;
+    values[kDCAzNPeaksTrimmed2] = -9999.0;
   }
-  auto [bimodalityBinTrimm3, meanBinTrimm3, stddevBinTrimm3, skewnessBinTrimm3, kurtosisBinTrimm3] = BimodalityCoefficient(dcazValues, 0.005, 4);
+  cout << "Bimodality coefficient (trimmed 2): " << bimodalityBinTrimm2 << ", mean: " << meanBinTrimm2 << ", stddev: " << stddevBinTrimm2 << ", skewness: " << skewnessBinTrimm2 << ", kurtosis: " << kurtosisBinTrimm2 << ", nPeaks: " << nPeaksTrimm2 << endl;
+  auto [bimodalityBinTrimm3, meanBinTrimm3, stddevBinTrimm3, skewnessBinTrimm3, kurtosisBinTrimm3, nPeaksTrimm3] = BimodalityCoefficient(dcazValues, 0.005, 4);
   if (stddevBinTrimm3 > -1.0) {
     values[kDCAzBimodalityCoefficientBinnedTrimmed3] = bimodalityBinTrimm3;
     values[kDCAzMeanBinnedTrimmed3] = meanBinTrimm3;
     values[kDCAzRMSBinnedTrimmed3] = stddevBinTrimm3;
+    values[kDCAzNPeaksTrimmed3] = nPeaksTrimm3;
   } else {
     values[kDCAzBimodalityCoefficientBinnedTrimmed3] = -9999.0;
     values[kDCAzMeanBinnedTrimmed3] = -9999.0;
     values[kDCAzRMSBinnedTrimmed3] = -9999.0;
+    values[kDCAzNPeaksTrimmed3] = -9999.0;
   }
+  cout << "Bimodality coefficient (trimmed 3): " << bimodalityBinTrimm3 << ", mean: " << meanBinTrimm3 << ", stddev: " << stddevBinTrimm3 << ", skewness: " << skewnessBinTrimm3 << ", kurtosis: " << kurtosisBinTrimm3 << ", nPeaks: " << nPeaksTrimm3 << endl;
 
   // compute fraction of tracks with |DCAz| > 100um, 200um, 500um, 1mm, 2mm, 5mm, 10mm
   // make a loop over the DCAz values and count how many are above each threshold

PWGDQ/DataModel/ReducedInfoTables.h

@@ -88,6 +88,10 @@ DECLARE_SOA_COLUMN(DCAzFracAbove1mm, dcazFracAbove1mm, float);
 DECLARE_SOA_COLUMN(DCAzFracAbove2mm, dcazFracAbove2mm, float);                   //!  Fraction of tracks in the event with |DCAz| > 2mm
 DECLARE_SOA_COLUMN(DCAzFracAbove5mm, dcazFracAbove5mm, float);                   //!  Fraction of tracks in the event with |DCAz| > 5mm
 DECLARE_SOA_COLUMN(DCAzFracAbove10mm, dcazFracAbove10mm, float);                 //!  Fraction of tracks in the event with |DCAz| > 10mm
+DECLARE_SOA_COLUMN(DCAzNPeaks, dcazNPeaks, int);                                 //!  Number of peaks in the DCAz distribution of the tracks in the event
+DECLARE_SOA_COLUMN(DCAzNPeaksTrimmed1, dcazNPeaksTrimmed1, int);                 //!  Number of peaks in the binned DCAz distribution (trimmed 1)
+DECLARE_SOA_COLUMN(DCAzNPeaksTrimmed2, dcazNPeaksTrimmed2, int);                 //!  Number of peaks in the binned DCAz distribution (trimmed 2)
+DECLARE_SOA_COLUMN(DCAzNPeaksTrimmed3, dcazNPeaksTrimmed3, int);                 //!  Number of peaks in the binned DCAz distribution (trimmed 3)
 
 // Columns declared to guarantee the backward compatibility of the tables
 DECLARE_SOA_COLUMN(QvecBPosRe, qvecBPosRe, float);
@@ -226,7 +230,8 @@ DECLARE_SOA_TABLE(ReducedEventsMergingTable, "AOD", "REMERGE", //!   Collision m
                   reducedevent::DCAzRMS, reducedevent::DCAzRMSBinnedTrimmed1, reducedevent::DCAzRMSBinnedTrimmed2, reducedevent::DCAzRMSBinnedTrimmed3,
                   reducedevent::DCAzSkewness, reducedevent::DCAzKurtosis,
                   reducedevent::DCAzFracAbove100um, reducedevent::DCAzFracAbove200um, reducedevent::DCAzFracAbove500um,
-                  reducedevent::DCAzFracAbove1mm, reducedevent::DCAzFracAbove2mm, reducedevent::DCAzFracAbove5mm, reducedevent::DCAzFracAbove10mm);
+                  reducedevent::DCAzFracAbove1mm, reducedevent::DCAzFracAbove2mm, reducedevent::DCAzFracAbove5mm, reducedevent::DCAzFracAbove10mm,
+                  reducedevent::DCAzNPeaks, reducedevent::DCAzNPeaksTrimmed1, reducedevent::DCAzNPeaksTrimmed2, reducedevent::DCAzNPeaksTrimmed3);
 
 // TODO and NOTE: This table is just an extension of the ReducedEvents table
 //       There is no explicit accounting for MC events which were not reconstructed!!!