diff --git a/PWGDQ/Core/VarManager.cxx b/PWGDQ/Core/VarManager.cxx index 1c385185f9a..af3bdcdefea 100644 --- a/PWGDQ/Core/VarManager.cxx +++ b/PWGDQ/Core/VarManager.cxx @@ -358,6 +358,135 @@ double VarManager::ComputePIDcalibration(int species, double nSigmaValue) } } +//__________________________________________________________________ +std::tuple VarManager::BimodalityCoefficientUnbinned(const std::vector& 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 VarManager::BimodalityCoefficientAndNPeaks(const std::vector& 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 result = BimodalityCoefficientUnbinned(data); + return std::make_tuple(std::get<0>(result), std::get<1>(result), std::get<2>(result), std::get<3>(result), std::get<4>(result), -1); + } + + // bin the data and put it in a vector + int nBins = static_cast((max - min) / binWidth); + std::vector counts(nBins, 0.0); + + for (float x : data) { + if (x < min || x >= max) { + continue; // skip out-of-range values + } + int bin = static_cast((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; diff --git a/PWGDQ/Core/VarManager.h b/PWGDQ/Core/VarManager.h index 49bdf4eef5b..cd3d31d7f05 100644 --- a/PWGDQ/Core/VarManager.h +++ b/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 @@ -293,6 +294,32 @@ class VarManager : public TObject kNTPCmeanTimeShortC, kNTPCmedianTimeShortA, kNTPCmedianTimeShortC, + kDCAzBimodalityCoefficient, + kDCAzMean, + kDCAzRMS, + kDCAzSkewness, + kDCAzKurtosis, + kDCAzBimodalityCoefficientBinned, + kDCAzBimodalityCoefficientBinnedTrimmed1, + kDCAzBimodalityCoefficientBinnedTrimmed2, + kDCAzBimodalityCoefficientBinnedTrimmed3, + kDCAzMeanBinnedTrimmed1, + kDCAzMeanBinnedTrimmed2, + kDCAzMeanBinnedTrimmed3, + kDCAzRMSBinnedTrimmed1, + kDCAzRMSBinnedTrimmed2, + kDCAzRMSBinnedTrimmed3, + kDCAzFracAbove100um, + kDCAzFracAbove200um, + kDCAzFracAbove500um, + kDCAzFracAbove1mm, + kDCAzFracAbove2mm, + kDCAzFracAbove5mm, + kDCAzFracAbove10mm, + kDCAzNPeaks, + kDCAzNPeaksTrimmed1, + kDCAzNPeaksTrimmed2, + kDCAzNPeaksTrimmed3, kMCEventGeneratorId, kMCEventSubGeneratorId, kMCVtxX, @@ -1255,6 +1282,10 @@ class VarManager : public TObject } template static o2::dataformats::VertexBase RecalculatePrimaryVertex(T const& track0, T const& track1, const T1& collision); + + static std::tuple BimodalityCoefficientUnbinned(const std::vector& data); + static std::tuple BimodalityCoefficientAndNPeaks(const std::vector& data, float binWidth, int trim = 0, float min = -15.0, float max = 15.0); + template static o2::track::TrackParCovFwd FwdToTrackPar(const T& track, const C& cov); template @@ -1270,6 +1301,8 @@ class VarManager : public TObject template static void FillEvent(T const& event, float* values = nullptr); template + static void FillEventTracks(T const& tracks, float* values = nullptr); + template static void FillTimeFrame(T const& tfTable, float* values = nullptr); template static void FillEventFlowResoFactor(T const& hs_sp, T const& hs_ep, float* values = nullptr); @@ -1486,7 +1519,7 @@ class VarManager : public TObject VarManager& operator=(const VarManager& c); VarManager(const VarManager& c); - ClassDef(VarManager, 5); + ClassDef(VarManager, 6); }; template @@ -1618,6 +1651,7 @@ o2::dataformats::VertexBase VarManager::RecalculatePrimaryVertex(T const& track0 o2::dataformats::VertexBase primaryVertexRec = {std::move(vtxXYZ), std::move(vtxCov)}; return primaryVertexRec; } + template o2::dataformats::GlobalFwdTrack VarManager::PropagateMuon(const T& muon, const C& collision, const int endPoint) { @@ -1839,6 +1873,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) @@ -2329,6 +2367,143 @@ void VarManager::FillEvent(T const& event, float* values) // FillEventDerived(values); } +template +void VarManager::FillEventTracks(T const& tracks, float* values) +{ + if (!values) { + values = fgValues; + } + + // compute event properties based on DCAz of the tracks + std::vector dcazValues; + for (const auto& track : tracks) { + if (!track.hasITS()) { + continue; // skip tracks without ITS information + } + if (!track.hasTPC()) { + continue; // skip tracks without TPC information + } + if (track.dcaZ() > 998) { + continue; // skip tracks without valid DCAz + } + dcazValues.push_back(track.dcaZ()); + } + + if (dcazValues.empty()) { + return; + } + + // compute the unbinned bimodality coefficient and related statistics + auto [bimodalityCoefficient, mean, stddev, skewness, kurtosis, nPeaks] = BimodalityCoefficientAndNPeaks(dcazValues, -1.0); + if (stddev > -1.0) { + values[kDCAzBimodalityCoefficient] = bimodalityCoefficient; + values[kDCAzMean] = mean; + values[kDCAzRMS] = stddev; + values[kDCAzSkewness] = skewness; + values[kDCAzKurtosis] = kurtosis; + values[kDCAzNPeaks] = nPeaks; + } 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 [bimodalityCoefficientBin, meanBin, stddevBin, skewnessBin, kurtosisBin, nPeaksBin] = BimodalityCoefficientAndNPeaks(dcazValues, 0.005); + if (stddevBin > -1.0) { + values[kDCAzBimodalityCoefficientBinned] = bimodalityCoefficientBin; + values[kDCAzNPeaks] = nPeaksBin; + } else { + values[kDCAzBimodalityCoefficientBinned] = -9999.0; + values[kDCAzNPeaks] = -9999.0; + } + cout << "Bimodality coefficient binned: " << bimodalityCoefficientBin << ", mean: " << mean << ", stddev: " << stddev << ", skewness: " << skewness << ", kurtosis: " << kurtosis << ", nPeaks: " << nPeaksBin << 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 [bimodalityCoefficientBinTrimmed1, meanBinTrimmed1, stddevBinTrimmed1, skewnessBinTrimmed1, kurtosisBinTrimmed1, nPeaksBinTrimmed1] = BimodalityCoefficientAndNPeaks(dcazValues, 0.005, 2); + if (stddevBinTrimmed1 > -1.0) { + values[kDCAzBimodalityCoefficientBinnedTrimmed1] = bimodalityCoefficientBinTrimmed1; + values[kDCAzMeanBinnedTrimmed1] = meanBinTrimmed1; + values[kDCAzRMSBinnedTrimmed1] = stddevBinTrimmed1; + values[kDCAzNPeaksTrimmed1] = nPeaksBinTrimmed1; + } else { + values[kDCAzBimodalityCoefficientBinnedTrimmed1] = -9999.0; + values[kDCAzMeanBinnedTrimmed1] = -9999.0; + values[kDCAzRMSBinnedTrimmed1] = -9999.0; + values[kDCAzNPeaksTrimmed1] = -9999.0; + } + cout << "Bimodality coefficient (trimmed 1): " << bimodalityCoefficientBinTrimmed1 << ", mean: " << meanBinTrimmed1 << ", stddev: " << stddevBinTrimmed1 << ", skewness: " << skewnessBinTrimmed1 << ", kurtosis: " << kurtosisBinTrimmed1 << ", nPeaks: " << nPeaksBinTrimmed1 << endl; + auto [bimodalityCoefficientBinTrimmed2, meanBinTrimmed2, stddevBinTrimmed2, skewnessBinTrimmed2, kurtosisBinTrimmed2, nPeaksBinTrimmed2] = BimodalityCoefficientAndNPeaks(dcazValues, 0.005, 3); + if (stddevBinTrimmed2 > -1.0) { + values[kDCAzBimodalityCoefficientBinnedTrimmed2] = bimodalityCoefficientBinTrimmed2; + values[kDCAzMeanBinnedTrimmed2] = meanBinTrimmed2; + values[kDCAzRMSBinnedTrimmed2] = stddevBinTrimmed2; + values[kDCAzNPeaksTrimmed2] = nPeaksBinTrimmed2; + } else { + values[kDCAzBimodalityCoefficientBinnedTrimmed2] = -9999.0; + values[kDCAzMeanBinnedTrimmed2] = -9999.0; + values[kDCAzRMSBinnedTrimmed2] = -9999.0; + values[kDCAzNPeaksTrimmed2] = -9999.0; + } + cout << "Bimodality coefficient (trimmed 2): " << bimodalityCoefficientBinTrimmed2 << ", mean: " << meanBinTrimmed2 << ", stddev: " << stddevBinTrimmed2 << ", skewness: " << skewnessBinTrimmed2 << ", kurtosis: " << kurtosisBinTrimmed2 << ", nPeaks: " << nPeaksBinTrimmed2 << endl; + auto [bimodalityCoefficientBinTrimmed3, meanBinTrimmed3, stddevBinTrimmed3, skewnessBinTrimmed3, kurtosisBinTrimmed3, nPeaksBinTrimmed3] = BimodalityCoefficientAndNPeaks(dcazValues, 0.005, 4); + if (stddevBinTrimmed3 > -1.0) { + values[kDCAzBimodalityCoefficientBinnedTrimmed3] = bimodalityCoefficientBinTrimmed3; + values[kDCAzMeanBinnedTrimmed3] = meanBinTrimmed3; + values[kDCAzRMSBinnedTrimmed3] = stddevBinTrimmed3; + values[kDCAzNPeaksTrimmed3] = nPeaksBinTrimmed3; + } else { + values[kDCAzBimodalityCoefficientBinnedTrimmed3] = -9999.0; + values[kDCAzMeanBinnedTrimmed3] = -9999.0; + values[kDCAzRMSBinnedTrimmed3] = -9999.0; + values[kDCAzNPeaksTrimmed3] = -9999.0; + } + cout << "Bimodality coefficient (trimmed 3): " << bimodalityCoefficientBinTrimmed3 << ", mean: " << meanBinTrimmed3 << ", stddev: " << stddevBinTrimmed3 << ", skewness: " << skewnessBinTrimmed3 << ", kurtosis: " << kurtosisBinTrimmed3 << ", nPeaks: " << nPeaksBinTrimmed3 << 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 + int counter100um = 0; + int counter200um = 0; + int counter500um = 0; + int counter1mm = 0; + int counter2mm = 0; + int counter5mm = 0; + int counter10mm = 0; + for (auto& d : dcazValues) { + double absD = std::abs(d); + if (absD > 0.01) { + counter100um++; + if (absD > 0.02) { + counter200um++; + if (absD > 0.05) { + counter500um++; + if (absD > 0.1) { + counter1mm++; + if (absD > 0.2) { + counter2mm++; + if (absD > 0.5) { + counter5mm++; + if (absD > 1.0) { + counter10mm++; + } + } + } + } + } + } + } + } + int totalTracks = static_cast(dcazValues.size()); + values[kDCAzFracAbove100um] = static_cast(counter100um) / totalTracks; + values[kDCAzFracAbove200um] = static_cast(counter200um) / totalTracks; + values[kDCAzFracAbove500um] = static_cast(counter500um) / totalTracks; + values[kDCAzFracAbove1mm] = static_cast(counter1mm) / totalTracks; + values[kDCAzFracAbove2mm] = static_cast(counter2mm) / totalTracks; + values[kDCAzFracAbove5mm] = static_cast(counter5mm) / totalTracks; + values[kDCAzFracAbove10mm] = static_cast(counter10mm) / totalTracks; +} + template void VarManager::FillEventFlowResoFactor(T const& hs_sp, T const& hs_ep, float* values) { diff --git a/PWGDQ/DataModel/ReducedInfoTables.h b/PWGDQ/DataModel/ReducedInfoTables.h index cb607ff02a3..94c22011d6e 100644 --- a/PWGDQ/DataModel/ReducedInfoTables.h +++ b/PWGDQ/DataModel/ReducedInfoTables.h @@ -66,6 +66,32 @@ DECLARE_SOA_COLUMN(NTPCoccupMeanTimeShortA, nTPCoccupMeanTimeShortA, float); DECLARE_SOA_COLUMN(NTPCoccupMeanTimeShortC, nTPCoccupMeanTimeShortC, float); //! TPC pileup mean time on C side (short time range) DECLARE_SOA_COLUMN(NTPCoccupMedianTimeShortA, nTPCoccupMedianTimeShortA, float); //! TPC pileup median time on A side (short time range) DECLARE_SOA_COLUMN(NTPCoccupMedianTimeShortC, nTPCoccupMedianTimeShortC, float); //! TPC pileup median time on C side (short time range) +DECLARE_SOA_COLUMN(DCAzBimodalityCoefficient, dcazBimodalityCoefficient, float); //! Bimodality coefficient of the DCAz distribution of the tracks in the event +DECLARE_SOA_COLUMN(DCAzBimodalityCoefficientBinned, dcazBimodalityCoefficientBinned, float); //! Bimodality coefficient of the DCAz distribution of the tracks in the event, binned +DECLARE_SOA_COLUMN(DCAzBimodalityCoefficientBinnedTrimmed1, dcazBimodalityCoefficientBinnedTrimmed1, float); //! Bimodality coefficient of the DCAz distribution of the tracks in the event, binned and trimmed 1 +DECLARE_SOA_COLUMN(DCAzBimodalityCoefficientBinnedTrimmed2, dcazBimodalityCoefficientBinnedTrimmed2, float); //! Bimodality coefficient of the DCAz distribution of the tracks in the event, binned and trimmed 2 +DECLARE_SOA_COLUMN(DCAzBimodalityCoefficientBinnedTrimmed3, dcazBimodalityCoefficientBinnedTrimmed3, float); //! Bimodality coefficient of the DCAz distribution of the tracks in the event, binned and trimmed 3 +DECLARE_SOA_COLUMN(DCAzMean, dcazMean, float); //! Mean of the DCAz distribution of the tracks in the event +DECLARE_SOA_COLUMN(DCAzMeanBinnedTrimmed1, dcazMeanBinnedTrimmed1, float); //! Mean of the DCAz distribution of the tracks in the event, binned and trimmed 1 +DECLARE_SOA_COLUMN(DCAzMeanBinnedTrimmed2, dcazMeanBinnedTrimmed2, float); //! Mean of the DCAz distribution of the tracks in the event, binned and trimmed 2 +DECLARE_SOA_COLUMN(DCAzMeanBinnedTrimmed3, dcazMeanBinnedTrimmed3, float); //! Mean of the DCAz distribution of the tracks in the event, binned and trimmed 3 +DECLARE_SOA_COLUMN(DCAzRMS, dcazRMS, float); //! RMS of the DCAz distribution of the tracks in the event +DECLARE_SOA_COLUMN(DCAzRMSBinnedTrimmed1, dcazRMSBinnedTrimmed1, float); //! RMS of the DCAz distribution of the tracks in the event, binned and trimmed 1 +DECLARE_SOA_COLUMN(DCAzRMSBinnedTrimmed2, dcazRMSBinnedTrimmed2, float); //! RMS of the DCAz distribution of the tracks in the event, binned and trimmed 2 +DECLARE_SOA_COLUMN(DCAzRMSBinnedTrimmed3, dcazRMSBinnedTrimmed3, float); //! RMS of the DCAz distribution of the tracks in the event, binned and trimmed 3 +DECLARE_SOA_COLUMN(DCAzSkewness, dcazSkewness, float); //! Skewness of the DCAz distribution of the tracks in the event +DECLARE_SOA_COLUMN(DCAzKurtosis, dcazKurtosis, float); //! Kurtosis of the DCAz distribution of the tracks in the event +DECLARE_SOA_COLUMN(DCAzFracAbove100um, dcazFracAbove100um, float); //! Fraction of tracks in the event with |DCAz| > 100um +DECLARE_SOA_COLUMN(DCAzFracAbove200um, dcazFracAbove200um, float); //! Fraction of tracks in the event with |DCAz| > 200um +DECLARE_SOA_COLUMN(DCAzFracAbove500um, dcazFracAbove500um, float); //! Fraction of tracks in the event with |DCAz| > 500um +DECLARE_SOA_COLUMN(DCAzFracAbove1mm, dcazFracAbove1mm, float); //! Fraction of tracks in the event with |DCAz| > 1mm +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); @@ -197,6 +223,16 @@ DECLARE_SOA_TABLE(ReducedEventsQvectorZN, "AOD", "REQVECTORZN", //! Event Q-v DECLARE_SOA_TABLE(ReducedEventsInfo, "AOD", "REDUCEVENTINFO", //! Main event index table reducedevent::CollisionId); +DECLARE_SOA_TABLE(ReducedEventsMergingTable, "AOD", "REMERGE", //! Collision merging quatities + reducedevent::DCAzBimodalityCoefficient, reducedevent::DCAzBimodalityCoefficientBinned, + reducedevent::DCAzBimodalityCoefficientBinnedTrimmed1, reducedevent::DCAzBimodalityCoefficientBinnedTrimmed2, reducedevent::DCAzBimodalityCoefficientBinnedTrimmed3, + reducedevent::DCAzMean, reducedevent::DCAzMeanBinnedTrimmed1, reducedevent::DCAzMeanBinnedTrimmed2, reducedevent::DCAzMeanBinnedTrimmed3, + 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::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!!! // However, for analysis which will require these events, a special skimming process function diff --git a/PWGDQ/TableProducer/tableMakerMC_withAssoc.cxx b/PWGDQ/TableProducer/tableMakerMC_withAssoc.cxx index 441287b2de0..49b4ca8d84a 100644 --- a/PWGDQ/TableProducer/tableMakerMC_withAssoc.cxx +++ b/PWGDQ/TableProducer/tableMakerMC_withAssoc.cxx @@ -143,6 +143,7 @@ struct TableMakerMC { Produces multPV; Produces multAll; Produces eventMClabels; + Produces mergingTable; Produces trackBarrelInfo; Produces trackBasic; @@ -273,6 +274,36 @@ struct TableMakerMC { // RCT flag checker RCTFlagsChecker rctChecker{"CBT"}; + // variables to store quantities needed for tagging collision merging candidates + struct { + std::map bimodalityCoeffDCAz; // Bimodality coefficient of the DCAz distribution of tracks associated to a collision + std::map bimodalityCoeffDCAzBinned; // Bimodality coefficient of the DCAz distribution of tracks associated to a collision, binned + std::map bimodalityCoeffDCAzBinnedTrimmed1; // Bimodality coefficient of the DCAz distribution of tracks associated to a collision, binned and trimmed 1 + std::map bimodalityCoeffDCAzBinnedTrimmed2; // Bimodality coefficient of the DCAz distribution of tracks associated to a collision, binned and trimmed 2 + std::map bimodalityCoeffDCAzBinnedTrimmed3; // Bimodality coefficient of the DCAz distribution of tracks associated to a collision, binned and trimmed 3 + std::map meanDCAz; + std::map meanDCAzBinnedTrimmed1; + std::map meanDCAzBinnedTrimmed2; + std::map meanDCAzBinnedTrimmed3; + std::map rmsDCAz; + std::map rmsDCAzBinnedTrimmed1; + std::map rmsDCAzBinnedTrimmed2; + std::map rmsDCAzBinnedTrimmed3; + std::map skewnessDCAz; + std::map kurtosisDCAz; + std::map fraction100umDCAz; // fraction of tracks with |DCAz|>100um + std::map fraction200umDCAz; // fraction of tracks with |DCAz|>200um + std::map fraction500umDCAz; // fraction of tracks with |DCAz|>500um + std::map fraction1mmDCAz; // fraction of tracks with |DCAz|>1mm + std::map fraction2mmDCAz; // fraction of tracks with |DCAz|>2mm + std::map fraction5mmDCAz; // fraction of tracks with |DCAz|>5mm + std::map fraction10mmDCAz; // fraction of tracks with |DCAz|>10mm + std::map nPeaksDCAz; // number of peaks in the DCAz distribution of tracks associated to a collision + std::map nPeaksDCAzTrimmed1; // number of peaks in the binned DCAz distribution (trimmed 1) + std::map nPeaksDCAzTrimmed2; // number of peaks in the binned DCAz distribution (trimmed 2) + std::map nPeaksDCAzTrimmed3; // number of peaks in the binned DCAz distribution (trimmed 3) + } fCollMergingTag; + void init(o2::framework::InitContext& context) { // Check whether barrel or muon are enabled @@ -472,6 +503,8 @@ struct TableMakerMC { Preslice fwdtrackIndicesPerCollision = aod::track_association::collisionId; Preslice mfttrackIndicesPerCollision = aod::track_association::collisionId; + Preslice presliceBarrelTracksWithCov = aod::track_association::collisionId; + void skimMCCollisions(MyEventsMcWithMults const& mcCollisions) { // skim MC collisions @@ -582,6 +615,72 @@ struct TableMakerMC { } // end loop over mc stack } + template + void computeCollMergingTag(TEvents const& collisions, TTracks const& tracks, Preslice& preslice) + { + // this function uses the standard track - collision association to compute several quantities which could be related to collision merging + // clear the maps for this time frame + fCollMergingTag.bimodalityCoeffDCAz.clear(); + fCollMergingTag.bimodalityCoeffDCAzBinned.clear(); + fCollMergingTag.bimodalityCoeffDCAzBinnedTrimmed1.clear(); + fCollMergingTag.bimodalityCoeffDCAzBinnedTrimmed2.clear(); + fCollMergingTag.bimodalityCoeffDCAzBinnedTrimmed3.clear(); + fCollMergingTag.meanDCAz.clear(); + fCollMergingTag.meanDCAzBinnedTrimmed1.clear(); + fCollMergingTag.meanDCAzBinnedTrimmed2.clear(); + fCollMergingTag.meanDCAzBinnedTrimmed3.clear(); + fCollMergingTag.rmsDCAz.clear(); + fCollMergingTag.rmsDCAzBinnedTrimmed1.clear(); + fCollMergingTag.rmsDCAzBinnedTrimmed2.clear(); + fCollMergingTag.rmsDCAzBinnedTrimmed3.clear(); + fCollMergingTag.skewnessDCAz.clear(); + fCollMergingTag.kurtosisDCAz.clear(); + fCollMergingTag.fraction100umDCAz.clear(); + fCollMergingTag.fraction200umDCAz.clear(); + fCollMergingTag.fraction500umDCAz.clear(); + fCollMergingTag.fraction1mmDCAz.clear(); + fCollMergingTag.fraction2mmDCAz.clear(); + fCollMergingTag.fraction5mmDCAz.clear(); + fCollMergingTag.fraction10mmDCAz.clear(); + fCollMergingTag.nPeaksDCAz.clear(); + fCollMergingTag.nPeaksDCAzTrimmed1.clear(); + fCollMergingTag.nPeaksDCAzTrimmed2.clear(); + fCollMergingTag.nPeaksDCAzTrimmed3.clear(); + + for (const auto& collision : collisions) { + // make a slice for this collision and compute the DCAz based event quantities + auto thisCollTracks = tracks.sliceBy(preslice, collision.globalIndex()); + VarManager::FillEventTracks(thisCollTracks); // fill the VarManager arrays with the information of the tracks associated to this collision, needed for the cuts and histograms + // add the computed variables to the maps with the collision index as key + fCollMergingTag.bimodalityCoeffDCAz[collision.globalIndex()] = VarManager::fgValues[VarManager::kDCAzBimodalityCoefficient]; + fCollMergingTag.bimodalityCoeffDCAzBinned[collision.globalIndex()] = VarManager::fgValues[VarManager::kDCAzBimodalityCoefficientBinned]; + fCollMergingTag.bimodalityCoeffDCAzBinnedTrimmed1[collision.globalIndex()] = VarManager::fgValues[VarManager::kDCAzBimodalityCoefficientBinnedTrimmed1]; + fCollMergingTag.bimodalityCoeffDCAzBinnedTrimmed2[collision.globalIndex()] = VarManager::fgValues[VarManager::kDCAzBimodalityCoefficientBinnedTrimmed2]; + fCollMergingTag.bimodalityCoeffDCAzBinnedTrimmed3[collision.globalIndex()] = VarManager::fgValues[VarManager::kDCAzBimodalityCoefficientBinnedTrimmed3]; + fCollMergingTag.meanDCAz[collision.globalIndex()] = VarManager::fgValues[VarManager::kDCAzMean]; + fCollMergingTag.meanDCAzBinnedTrimmed1[collision.globalIndex()] = VarManager::fgValues[VarManager::kDCAzMeanBinnedTrimmed1]; + fCollMergingTag.meanDCAzBinnedTrimmed2[collision.globalIndex()] = VarManager::fgValues[VarManager::kDCAzMeanBinnedTrimmed2]; + fCollMergingTag.meanDCAzBinnedTrimmed3[collision.globalIndex()] = VarManager::fgValues[VarManager::kDCAzMeanBinnedTrimmed3]; + fCollMergingTag.rmsDCAz[collision.globalIndex()] = VarManager::fgValues[VarManager::kDCAzRMS]; + fCollMergingTag.rmsDCAzBinnedTrimmed1[collision.globalIndex()] = VarManager::fgValues[VarManager::kDCAzRMSBinnedTrimmed1]; + fCollMergingTag.rmsDCAzBinnedTrimmed2[collision.globalIndex()] = VarManager::fgValues[VarManager::kDCAzRMSBinnedTrimmed2]; + fCollMergingTag.rmsDCAzBinnedTrimmed3[collision.globalIndex()] = VarManager::fgValues[VarManager::kDCAzRMSBinnedTrimmed3]; + fCollMergingTag.skewnessDCAz[collision.globalIndex()] = VarManager::fgValues[VarManager::kDCAzSkewness]; + fCollMergingTag.kurtosisDCAz[collision.globalIndex()] = VarManager::fgValues[VarManager::kDCAzKurtosis]; + fCollMergingTag.fraction100umDCAz[collision.globalIndex()] = VarManager::fgValues[VarManager::kDCAzFracAbove100um]; + fCollMergingTag.fraction200umDCAz[collision.globalIndex()] = VarManager::fgValues[VarManager::kDCAzFracAbove200um]; + fCollMergingTag.fraction500umDCAz[collision.globalIndex()] = VarManager::fgValues[VarManager::kDCAzFracAbove500um]; + fCollMergingTag.fraction1mmDCAz[collision.globalIndex()] = VarManager::fgValues[VarManager::kDCAzFracAbove1mm]; + fCollMergingTag.fraction2mmDCAz[collision.globalIndex()] = VarManager::fgValues[VarManager::kDCAzFracAbove2mm]; + fCollMergingTag.fraction5mmDCAz[collision.globalIndex()] = VarManager::fgValues[VarManager::kDCAzFracAbove5mm]; + fCollMergingTag.fraction10mmDCAz[collision.globalIndex()] = VarManager::fgValues[VarManager::kDCAzFracAbove10mm]; + fCollMergingTag.nPeaksDCAz[collision.globalIndex()] = static_cast(VarManager::fgValues[VarManager::kDCAzNPeaks]); + fCollMergingTag.nPeaksDCAzTrimmed1[collision.globalIndex()] = static_cast(VarManager::fgValues[VarManager::kDCAzNPeaksTrimmed1]); + fCollMergingTag.nPeaksDCAzTrimmed2[collision.globalIndex()] = static_cast(VarManager::fgValues[VarManager::kDCAzNPeaksTrimmed2]); + fCollMergingTag.nPeaksDCAzTrimmed3[collision.globalIndex()] = static_cast(VarManager::fgValues[VarManager::kDCAzNPeaksTrimmed3]); + } + } + template void skimCollisions(TEvents const& collisions, BCsWithTimestamps const& /*bcs*/) { @@ -698,6 +797,17 @@ struct TableMakerMC { 0, 0, 0.0, 0.0, 0.0, 0.0, 0, 0, 0.0, 0.0, 0.0, 0.0); } + mergingTable(fCollMergingTag.bimodalityCoeffDCAz[collision.globalIndex()], fCollMergingTag.bimodalityCoeffDCAzBinned[collision.globalIndex()], + fCollMergingTag.bimodalityCoeffDCAzBinnedTrimmed1[collision.globalIndex()], fCollMergingTag.bimodalityCoeffDCAzBinnedTrimmed2[collision.globalIndex()], fCollMergingTag.bimodalityCoeffDCAzBinnedTrimmed3[collision.globalIndex()], + fCollMergingTag.meanDCAz[collision.globalIndex()], fCollMergingTag.meanDCAzBinnedTrimmed1[collision.globalIndex()], fCollMergingTag.meanDCAzBinnedTrimmed2[collision.globalIndex()], fCollMergingTag.meanDCAzBinnedTrimmed3[collision.globalIndex()], + fCollMergingTag.rmsDCAz[collision.globalIndex()], fCollMergingTag.rmsDCAzBinnedTrimmed1[collision.globalIndex()], fCollMergingTag.rmsDCAzBinnedTrimmed2[collision.globalIndex()], fCollMergingTag.rmsDCAzBinnedTrimmed3[collision.globalIndex()], + fCollMergingTag.skewnessDCAz[collision.globalIndex()], fCollMergingTag.kurtosisDCAz[collision.globalIndex()], + fCollMergingTag.fraction100umDCAz[collision.globalIndex()], fCollMergingTag.fraction200umDCAz[collision.globalIndex()], + fCollMergingTag.fraction500umDCAz[collision.globalIndex()], fCollMergingTag.fraction1mmDCAz[collision.globalIndex()], + fCollMergingTag.fraction2mmDCAz[collision.globalIndex()], fCollMergingTag.fraction5mmDCAz[collision.globalIndex()], + fCollMergingTag.fraction10mmDCAz[collision.globalIndex()], fCollMergingTag.nPeaksDCAz[collision.globalIndex()], fCollMergingTag.nPeaksDCAzTrimmed1[collision.globalIndex()], + fCollMergingTag.nPeaksDCAzTrimmed2[collision.globalIndex()], fCollMergingTag.nPeaksDCAzTrimmed3[collision.globalIndex()]); + // add an element for this collision into the map fCollIndexMap[collision.globalIndex()] = event.lastIndex(); } @@ -1518,6 +1628,7 @@ struct TableMakerMC { MyBarrelTracksWithCov const& tracksBarrel, aod::TrackAssoc const& trackAssocs, MyEventsMcWithMults const& mcCollisions, aod::McParticles const& mcParticles) { + computeCollMergingTag(collisions, tracksBarrel, presliceBarrelTracksWithCov); fullSkimming(collisions, bcs, tracksBarrel, nullptr, nullptr, trackAssocs, nullptr, nullptr, mcCollisions, mcParticles, nullptr); } @@ -1525,6 +1636,7 @@ struct TableMakerMC { MyBarrelTracksWithCov const& tracksBarrel, aod::TrackAssoc const& trackAssocs, MyEventsMcWithMults const& mcCollisions, aod::McParticles const& mcParticles) { + computeCollMergingTag(collisions, tracksBarrel, presliceBarrelTracksWithCov); fullSkimming(collisions, bcs, tracksBarrel, nullptr, nullptr, trackAssocs, nullptr, nullptr, mcCollisions, mcParticles, nullptr); } diff --git a/PWGDQ/TableProducer/tableMaker_withAssoc.cxx b/PWGDQ/TableProducer/tableMaker_withAssoc.cxx index dcdc55e61c4..f6df99aa368 100644 --- a/PWGDQ/TableProducer/tableMaker_withAssoc.cxx +++ b/PWGDQ/TableProducer/tableMaker_withAssoc.cxx @@ -175,6 +175,7 @@ struct TableMaker { Produces fit; Produces multPV; Produces multAll; + Produces mergingTable; Produces trackBarrelInfo; Produces trackBasic; Produces trackBarrel; @@ -373,6 +374,36 @@ struct TableMaker { std::map oContribLongC; } fOccup; + // variables to store quantities needed for tagging collision merging candidates + struct { + std::map bimodalityCoeffDCAz; // Bimodality coefficient of the DCAz distribution of tracks associated to a collision + std::map bimodalityCoeffDCAzBinned; // Bimodality coefficient of the DCAz distribution of tracks associated to a collision, binned + std::map bimodalityCoeffDCAzBinnedTrimmed1; // Bimodality coefficient of the DCAz distribution of tracks associated to a collision, binned and trimmed 1 + std::map bimodalityCoeffDCAzBinnedTrimmed2; // Bimodality coefficient of the DCAz distribution of tracks associated to a collision, binned and trimmed 2 + std::map bimodalityCoeffDCAzBinnedTrimmed3; // Bimodality coefficient of the DCAz distribution of tracks associated to a collision, binned and trimmed 3 + std::map meanDCAz; + std::map meanDCAzBinnedTrimmed1; + std::map meanDCAzBinnedTrimmed2; + std::map meanDCAzBinnedTrimmed3; + std::map rmsDCAz; + std::map rmsDCAzBinnedTrimmed1; + std::map rmsDCAzBinnedTrimmed2; + std::map rmsDCAzBinnedTrimmed3; + std::map skewnessDCAz; + std::map kurtosisDCAz; + std::map fraction100umDCAz; // fraction of tracks with |DCAz|>100um + std::map fraction200umDCAz; // fraction of tracks with |DCAz|>200um + std::map fraction500umDCAz; // fraction of tracks with |DCAz|>500um + std::map fraction1mmDCAz; // fraction of tracks with |DCAz|>1mm + std::map fraction2mmDCAz; // fraction of tracks with |DCAz|>2mm + std::map fraction5mmDCAz; // fraction of tracks with |DCAz|>5mm + std::map fraction10mmDCAz; // fraction of tracks with |DCAz|>10mm + std::map nPeaksDCAz; // number of peaks in the DCAz distribution of tracks associated to a collision + std::map nPeaksDCAzTrimmed1; // number of peaks in the binned DCAz distribution (trimmed 1) + std::map nPeaksDCAzTrimmed2; // number of peaks in the binned DCAz distribution (trimmed 2) + std::map nPeaksDCAzTrimmed3; // number of peaks in the binned DCAz distribution (trimmed 3) + } fCollMergingTag; + void init(o2::framework::InitContext& context) { // CCDB configuration @@ -851,6 +882,72 @@ struct TableMaker { return mu; } + template + void computeCollMergingTag(TEvents const& collisions, TTracks const& tracks, Preslice& preslice) + { + // This function uses the standard track-collision association to compute quantities related to collision merging + // clear the maps for this time frame + fCollMergingTag.bimodalityCoeffDCAz.clear(); + fCollMergingTag.bimodalityCoeffDCAzBinned.clear(); + fCollMergingTag.bimodalityCoeffDCAzBinnedTrimmed1.clear(); + fCollMergingTag.bimodalityCoeffDCAzBinnedTrimmed2.clear(); + fCollMergingTag.bimodalityCoeffDCAzBinnedTrimmed3.clear(); + fCollMergingTag.meanDCAz.clear(); + fCollMergingTag.meanDCAzBinnedTrimmed1.clear(); + fCollMergingTag.meanDCAzBinnedTrimmed2.clear(); + fCollMergingTag.meanDCAzBinnedTrimmed3.clear(); + fCollMergingTag.rmsDCAz.clear(); + fCollMergingTag.rmsDCAzBinnedTrimmed1.clear(); + fCollMergingTag.rmsDCAzBinnedTrimmed2.clear(); + fCollMergingTag.rmsDCAzBinnedTrimmed3.clear(); + fCollMergingTag.skewnessDCAz.clear(); + fCollMergingTag.kurtosisDCAz.clear(); + fCollMergingTag.fraction100umDCAz.clear(); + fCollMergingTag.fraction200umDCAz.clear(); + fCollMergingTag.fraction500umDCAz.clear(); + fCollMergingTag.fraction1mmDCAz.clear(); + fCollMergingTag.fraction2mmDCAz.clear(); + fCollMergingTag.fraction5mmDCAz.clear(); + fCollMergingTag.fraction10mmDCAz.clear(); + fCollMergingTag.nPeaksDCAz.clear(); + fCollMergingTag.nPeaksDCAzTrimmed1.clear(); + fCollMergingTag.nPeaksDCAzTrimmed2.clear(); + fCollMergingTag.nPeaksDCAzTrimmed3.clear(); + + for (const auto& collision : collisions) { + // make a slice for this collision and compute the DCAz based event quantities + auto thisCollTracks = tracks.sliceBy(preslice, collision.globalIndex()); + VarManager::FillEventTracks(thisCollTracks); // fill the VarManager arrays with the information of the tracks associated to this collision, needed for the cuts and histograms + // add the computed variables to the maps with the collision index as key + fCollMergingTag.bimodalityCoeffDCAz[collision.globalIndex()] = VarManager::fgValues[VarManager::kDCAzBimodalityCoefficient]; + fCollMergingTag.bimodalityCoeffDCAzBinned[collision.globalIndex()] = VarManager::fgValues[VarManager::kDCAzBimodalityCoefficientBinned]; + fCollMergingTag.bimodalityCoeffDCAzBinnedTrimmed1[collision.globalIndex()] = VarManager::fgValues[VarManager::kDCAzBimodalityCoefficientBinnedTrimmed1]; + fCollMergingTag.bimodalityCoeffDCAzBinnedTrimmed2[collision.globalIndex()] = VarManager::fgValues[VarManager::kDCAzBimodalityCoefficientBinnedTrimmed2]; + fCollMergingTag.bimodalityCoeffDCAzBinnedTrimmed3[collision.globalIndex()] = VarManager::fgValues[VarManager::kDCAzBimodalityCoefficientBinnedTrimmed3]; + fCollMergingTag.meanDCAz[collision.globalIndex()] = VarManager::fgValues[VarManager::kDCAzMean]; + fCollMergingTag.meanDCAzBinnedTrimmed1[collision.globalIndex()] = VarManager::fgValues[VarManager::kDCAzMeanBinnedTrimmed1]; + fCollMergingTag.meanDCAzBinnedTrimmed2[collision.globalIndex()] = VarManager::fgValues[VarManager::kDCAzMeanBinnedTrimmed2]; + fCollMergingTag.meanDCAzBinnedTrimmed3[collision.globalIndex()] = VarManager::fgValues[VarManager::kDCAzMeanBinnedTrimmed3]; + fCollMergingTag.rmsDCAz[collision.globalIndex()] = VarManager::fgValues[VarManager::kDCAzRMS]; + fCollMergingTag.rmsDCAzBinnedTrimmed1[collision.globalIndex()] = VarManager::fgValues[VarManager::kDCAzRMSBinnedTrimmed1]; + fCollMergingTag.rmsDCAzBinnedTrimmed2[collision.globalIndex()] = VarManager::fgValues[VarManager::kDCAzRMSBinnedTrimmed2]; + fCollMergingTag.rmsDCAzBinnedTrimmed3[collision.globalIndex()] = VarManager::fgValues[VarManager::kDCAzRMSBinnedTrimmed3]; + fCollMergingTag.skewnessDCAz[collision.globalIndex()] = VarManager::fgValues[VarManager::kDCAzSkewness]; + fCollMergingTag.kurtosisDCAz[collision.globalIndex()] = VarManager::fgValues[VarManager::kDCAzKurtosis]; + fCollMergingTag.fraction100umDCAz[collision.globalIndex()] = VarManager::fgValues[VarManager::kDCAzFracAbove100um]; + fCollMergingTag.fraction200umDCAz[collision.globalIndex()] = VarManager::fgValues[VarManager::kDCAzFracAbove200um]; + fCollMergingTag.fraction500umDCAz[collision.globalIndex()] = VarManager::fgValues[VarManager::kDCAzFracAbove500um]; + fCollMergingTag.fraction1mmDCAz[collision.globalIndex()] = VarManager::fgValues[VarManager::kDCAzFracAbove1mm]; + fCollMergingTag.fraction2mmDCAz[collision.globalIndex()] = VarManager::fgValues[VarManager::kDCAzFracAbove2mm]; + fCollMergingTag.fraction5mmDCAz[collision.globalIndex()] = VarManager::fgValues[VarManager::kDCAzFracAbove5mm]; + fCollMergingTag.fraction10mmDCAz[collision.globalIndex()] = VarManager::fgValues[VarManager::kDCAzFracAbove10mm]; + fCollMergingTag.nPeaksDCAz[collision.globalIndex()] = VarManager::fgValues[VarManager::kDCAzNPeaks]; + fCollMergingTag.nPeaksDCAzTrimmed1[collision.globalIndex()] = VarManager::fgValues[VarManager::kDCAzNPeaksTrimmed1]; + fCollMergingTag.nPeaksDCAzTrimmed2[collision.globalIndex()] = VarManager::fgValues[VarManager::kDCAzNPeaksTrimmed2]; + fCollMergingTag.nPeaksDCAzTrimmed3[collision.globalIndex()] = VarManager::fgValues[VarManager::kDCAzNPeaksTrimmed3]; + } + } + template void skimCollisions(TEvents const& collisions, TBCs const& bcs, TZdcs const& /*zdcs*/, @@ -993,6 +1090,7 @@ struct TableMaker { VarManager::fgValues[VarManager::kNTPCmedianTimeShortA] = fOccup.oMedianTimeShortA[collision.globalIndex()]; VarManager::fgValues[VarManager::kNTPCmedianTimeShortC] = fOccup.oMedianTimeShortC[collision.globalIndex()]; } + if (fDoDetailedQA) { fHistMan->FillHistClass("Event_BeforeCuts", VarManager::fgValues); } @@ -1118,7 +1216,19 @@ struct TableMaker { fOccup.oMeanTimeShortA[collision.globalIndex()], fOccup.oMeanTimeShortC[collision.globalIndex()], fOccup.oMedianTimeShortA[collision.globalIndex()], fOccup.oMedianTimeShortC[collision.globalIndex()]); } - + mergingTable(fCollMergingTag.bimodalityCoeffDCAz[collision.globalIndex()], fCollMergingTag.bimodalityCoeffDCAzBinned[collision.globalIndex()], + fCollMergingTag.bimodalityCoeffDCAzBinnedTrimmed1[collision.globalIndex()], fCollMergingTag.bimodalityCoeffDCAzBinnedTrimmed2[collision.globalIndex()], fCollMergingTag.bimodalityCoeffDCAzBinnedTrimmed3[collision.globalIndex()], + fCollMergingTag.meanDCAz[collision.globalIndex()], fCollMergingTag.meanDCAzBinnedTrimmed1[collision.globalIndex()], fCollMergingTag.meanDCAzBinnedTrimmed2[collision.globalIndex()], fCollMergingTag.meanDCAzBinnedTrimmed3[collision.globalIndex()], + fCollMergingTag.rmsDCAz[collision.globalIndex()], fCollMergingTag.rmsDCAzBinnedTrimmed1[collision.globalIndex()], fCollMergingTag.rmsDCAzBinnedTrimmed2[collision.globalIndex()], fCollMergingTag.rmsDCAzBinnedTrimmed3[collision.globalIndex()], + fCollMergingTag.skewnessDCAz[collision.globalIndex()], fCollMergingTag.kurtosisDCAz[collision.globalIndex()], + fCollMergingTag.fraction100umDCAz[collision.globalIndex()], fCollMergingTag.fraction200umDCAz[collision.globalIndex()], + fCollMergingTag.fraction500umDCAz[collision.globalIndex()], fCollMergingTag.fraction1mmDCAz[collision.globalIndex()], + fCollMergingTag.fraction2mmDCAz[collision.globalIndex()], fCollMergingTag.fraction5mmDCAz[collision.globalIndex()], + fCollMergingTag.fraction10mmDCAz[collision.globalIndex()], + fCollMergingTag.nPeaksDCAz[collision.globalIndex()], fCollMergingTag.nPeaksDCAzTrimmed1[collision.globalIndex()], + fCollMergingTag.nPeaksDCAzTrimmed2[collision.globalIndex()], fCollMergingTag.nPeaksDCAzTrimmed3[collision.globalIndex()]); + + // fCollIndexMap[collision.globalIndex()] = event.lastIndex(); } } @@ -1783,6 +1893,8 @@ struct TableMaker { MyBarrelTracksWithCov const& tracksBarrel, TrackAssoc const& trackAssocs) { + computeOccupancyEstimators(collisions, tracksPosWithCov, tracksNegWithCov, presliceWithCov, bcs); + computeCollMergingTag(collisions, tracksBarrel, presliceWithCov); fullSkimming(collisions, bcs, nullptr, tracksBarrel, nullptr, nullptr, trackAssocs, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr); } @@ -1792,6 +1904,7 @@ struct TableMaker { TrackAssoc const& trackAssocs) { computeOccupancyEstimators(collisions, tracksPosWithCovNoTOF, tracksNegWithCovNoTOF, presliceWithCovNoTOF, bcs); + computeCollMergingTag(collisions, tracksBarrel, presliceWithCovNoTOF); fullSkimming(collisions, bcs, nullptr, tracksBarrel, nullptr, nullptr, trackAssocs, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr); } @@ -1801,6 +1914,7 @@ struct TableMaker { TrackAssoc const& trackAssocs, aod::FT0s& ft0s, aod::FV0As& fv0as, aod::FDDs& fdds) { computeOccupancyEstimators(collisions, tracksPosWithCov, tracksNegWithCov, presliceWithCov, bcs); + computeCollMergingTag(collisions, tracksBarrel, presliceWithCov); fullSkimming(collisions, bcs, zdcs, tracksBarrel, nullptr, nullptr, trackAssocs, nullptr, nullptr, nullptr, ft0s, fv0as, fdds); } @@ -1810,6 +1924,7 @@ struct TableMaker { TrackAssoc const& trackAssocs) { computeOccupancyEstimators(collisions, tracksPos, tracksNeg, preslice, bcs); + computeCollMergingTag(collisions, tracksBarrel, preslice); fullSkimming(collisions, bcs, nullptr, tracksBarrel, nullptr, nullptr, trackAssocs, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr); } @@ -1819,6 +1934,7 @@ struct TableMaker { TrackAssoc const& trackAssocs) { computeOccupancyEstimators(collisions, tracksPosNoTOF, tracksNegNoTOF, presliceNoTOF, bcs); + computeCollMergingTag(collisions, tracksBarrel, presliceNoTOF); fullSkimming(collisions, bcs, nullptr, tracksBarrel, nullptr, nullptr, trackAssocs, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr); }