[PWGDQ] Adding a table mergeable with the event tables, with quantities for detecting bimodality in the DCAz distributions

PWGDQ/Core/VarManager.cxx

@@ -358,6 +358,135 @@ double VarManager::ComputePIDcalibration(int species, double nSigmaValue)
   }
 }
 
+//__________________________________________________________________
+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);
+  }
+  float mean = std::accumulate(data.begin(), data.end(), 0.0) / n;
+
+  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 += diff2 * diff;
+    m4 += diff2 * diff2;
+  }
+
+  m2 /= n;
+  m3 /= n;
+  m4 /= n;
+
+  if (m2 == 0.0) {
+    return std::make_tuple(-1.0, -1.0, -1.0, -1.0, -1.0);
+  }
+
+  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<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) {
+    // 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
+  int nBins = static_cast<int>((max - min) / binWidth);
+  std::vector<int> counts(nBins, 0.0);
+
+  for (float x : data) {
+    if (x < min || x >= max) {
+      continue; // skip out-of-range values
+    }
+    int bin = static_cast<int>((x - min) / binWidth);
+    if (bin >= 0 && bin < nBins) {
+      counts[bin]++;
+    }
+  }
+
+  // trim the distribution if requested, by requiring a minimum of "trim" counts in each bin
+  if (trim > 0) {
+    for (int i = 0; i < nBins; ++i) {
+      if (counts[i] < trim) {
+        counts[i] = 0;
+      }
+    }
+  }
+
+  // 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
+  float mean = 0.0;
+  float totalCounts = 0.0;
+  for (int i = 0; i < nBins; ++i) {
+    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, -1);
+  }
+  mean /= totalCounts;
+
+  // then compute the second, third, and fourth central moments
+  float m2 = 0.0, m3 = 0.0, m4 = 0.0;
+  float diff, diff2, binCenter;
+  for (int i = 0; i < nBins; ++i) {
+    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] * 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, -1);
+  }
+
+  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, nPeaks);
+}
+
 //__________________________________________________________________
 void VarManager::SetDefaultVarNames()
 {
@@ -389,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 ";
@@ -577,6 +708,58 @@ void VarManager::SetDefaultVarNames()
   fgVariableUnits[kNTPCmedianTimeShortA] = "#mu s";
   fgVariableNames[kNTPCmedianTimeShortC] = "# TPC-C pileup median time, short time range";
   fgVariableUnits[kNTPCmedianTimeShortC] = "#mu s";
+  fgVariableNames[kDCAzBimodalityCoefficient] = "Unbinned Bimodality Coeff of DCAz distribution";
+  fgVariableUnits[kDCAzBimodalityCoefficient] = "";
+  fgVariableNames[kDCAzBimodalityCoefficientBinned] = "Binned Bimodality Coeff of DCAz distribution";
+  fgVariableUnits[kDCAzBimodalityCoefficientBinned] = "";
+  fgVariableNames[kDCAzBimodalityCoefficientBinnedTrimmed1] = "Binned Bimodality Coeff of DCAz distribution (trimmed 1)";
+  fgVariableUnits[kDCAzBimodalityCoefficientBinnedTrimmed1] = "";
+  fgVariableNames[kDCAzBimodalityCoefficientBinnedTrimmed2] = "Binned Bimodality Coeff of DCAz distribution (trimmed 2)";
+  fgVariableUnits[kDCAzBimodalityCoefficientBinnedTrimmed2] = "";
+  fgVariableNames[kDCAzBimodalityCoefficientBinnedTrimmed3] = "Binned Bimodality Coeff of DCAz distribution (trimmed 3)";
+  fgVariableUnits[kDCAzBimodalityCoefficientBinnedTrimmed3] = "";
+  fgVariableNames[kDCAzMean] = "Mean of DCAz distribution";
+  fgVariableUnits[kDCAzMean] = "cm";
+  fgVariableNames[kDCAzMeanBinnedTrimmed1] = "Mean of DCAz distribution (trimmed 1)";
+  fgVariableUnits[kDCAzMeanBinnedTrimmed1] = "cm";
+  fgVariableNames[kDCAzMeanBinnedTrimmed2] = "Mean of DCAz distribution (trimmed 2)";
+  fgVariableUnits[kDCAzMeanBinnedTrimmed2] = "cm";
+  fgVariableNames[kDCAzMeanBinnedTrimmed3] = "Mean of DCAz distribution (trimmed 3)";
+  fgVariableUnits[kDCAzMeanBinnedTrimmed3] = "cm";
+  fgVariableNames[kDCAzRMS] = "RMS of DCAz distribution";
+  fgVariableUnits[kDCAzRMS] = "cm";
+  fgVariableNames[kDCAzRMSBinnedTrimmed1] = "RMS of DCAz distribution (trimmed 1)";
+  fgVariableUnits[kDCAzRMSBinnedTrimmed1] = "cm";
+  fgVariableNames[kDCAzRMSBinnedTrimmed2] = "RMS of DCAz distribution (trimmed 2)";
+  fgVariableUnits[kDCAzRMSBinnedTrimmed2] = "cm";
+  fgVariableNames[kDCAzRMSBinnedTrimmed3] = "RMS of DCAz distribution (trimmed 3)";
+  fgVariableUnits[kDCAzRMSBinnedTrimmed3] = "cm";
+  fgVariableNames[kDCAzSkewness] = "Skewness of DCAz distribution";
+  fgVariableUnits[kDCAzSkewness] = "";
+  fgVariableNames[kDCAzKurtosis] = "Kurtosis of DCAz distribution";
+  fgVariableUnits[kDCAzKurtosis] = "";
+  fgVariableNames[kDCAzFracAbove100um] = "Fraction of tracks with |DCAz| > 100 um";
+  fgVariableUnits[kDCAzFracAbove100um] = "";
+  fgVariableNames[kDCAzFracAbove200um] = "Fraction of tracks with |DCAz| > 200 um";
+  fgVariableUnits[kDCAzFracAbove200um] = "";
+  fgVariableNames[kDCAzFracAbove500um] = "Fraction of tracks with |DCAz| > 500 um";
+  fgVariableUnits[kDCAzFracAbove500um] = "";
+  fgVariableNames[kDCAzFracAbove1mm] = "Fraction of tracks with |DCAz| > 1 mm";
+  fgVariableUnits[kDCAzFracAbove1mm] = "";
+  fgVariableNames[kDCAzFracAbove2mm] = "Fraction of tracks with |DCAz| > 2 mm";
+  fgVariableUnits[kDCAzFracAbove2mm] = "";
+  fgVariableNames[kDCAzFracAbove5mm] = "Fraction of tracks with |DCAz| > 5 mm";
+  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}";
@@ -1553,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;
@@ -1638,6 +1822,32 @@ void VarManager::SetDefaultVarNames()
   fgVarNamesMap["kNTPCmeanTimeShortC"] = kNTPCmeanTimeShortC;
   fgVarNamesMap["kNTPCmedianTimeShortA"] = kNTPCmedianTimeShortA;
   fgVarNamesMap["kNTPCmedianTimeShortC"] = kNTPCmedianTimeShortC;
+  fgVarNamesMap["kDCAzBimodalityCoefficient"] = kDCAzBimodalityCoefficient;
+  fgVarNamesMap["kDCAzBimodalityCoefficientBinned"] = kDCAzBimodalityCoefficientBinned;
+  fgVarNamesMap["kDCAzBimodalityCoefficientBinnedTrimmed1"] = kDCAzBimodalityCoefficientBinnedTrimmed1;
+  fgVarNamesMap["kDCAzBimodalityCoefficientBinnedTrimmed2"] = kDCAzBimodalityCoefficientBinnedTrimmed2;
+  fgVarNamesMap["kDCAzBimodalityCoefficientBinnedTrimmed3"] = kDCAzBimodalityCoefficientBinnedTrimmed3;
+  fgVarNamesMap["kDCAzMean"] = kDCAzMean;
+  fgVarNamesMap["kDCAzMeanBinnedTrimmed1"] = kDCAzMeanBinnedTrimmed1;
+  fgVarNamesMap["kDCAzMeanBinnedTrimmed2"] = kDCAzMeanBinnedTrimmed2;
+  fgVarNamesMap["kDCAzMeanBinnedTrimmed3"] = kDCAzMeanBinnedTrimmed3;
+  fgVarNamesMap["kDCAzRMS"] = kDCAzRMS;
+  fgVarNamesMap["kDCAzRMSBinnedTrimmed1"] = kDCAzRMSBinnedTrimmed1;
+  fgVarNamesMap["kDCAzRMSBinnedTrimmed2"] = kDCAzRMSBinnedTrimmed2;
+  fgVarNamesMap["kDCAzRMSBinnedTrimmed3"] = kDCAzRMSBinnedTrimmed3;
+  fgVarNamesMap["kDCAzSkewness"] = kDCAzSkewness;
+  fgVarNamesMap["kDCAzKurtosis"] = kDCAzKurtosis;
+  fgVarNamesMap["kDCAzFracAbove100um"] = kDCAzFracAbove100um;
+  fgVarNamesMap["kDCAzFracAbove200um"] = kDCAzFracAbove200um;
+  fgVarNamesMap["kDCAzFracAbove500um"] = kDCAzFracAbove500um;
+  fgVarNamesMap["kDCAzFracAbove1mm"] = kDCAzFracAbove1mm;
+  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;