From 9d5cd63c89c9e03ea95827926077887582132e09 Mon Sep 17 00:00:00 2001 From: Shiqi Wang <3326561594@qq.com> Date: Tue, 18 Aug 2026 16:50:33 +0800 Subject: [PATCH 1/3] Add hist for 3rd order ESE calculation --- PWGCF/Flow/Tasks/flowEseTask.cxx | 489 +++++++++++++------------------ 1 file changed, 207 insertions(+), 282 deletions(-) diff --git a/PWGCF/Flow/Tasks/flowEseTask.cxx b/PWGCF/Flow/Tasks/flowEseTask.cxx index e36ad406927..33379961c96 100644 --- a/PWGCF/Flow/Tasks/flowEseTask.cxx +++ b/PWGCF/Flow/Tasks/flowEseTask.cxx @@ -46,18 +46,22 @@ #include // IWYU pragma: keep (do not replace with Math/Vector4Dfwd.h) #include #include +#include #include #include #include #include #include +#include #include #include #include +#include #include #include #include +#include #include using namespace o2; @@ -120,6 +124,7 @@ struct FlowEseTask { Configurable cfgnMods{"cfgnMods", 1, "The number of modulations of interest starting from 2"}; Configurable cfgNQvec{"cfgNQvec", 7, "The number of total Qvectors for looping over the task"}; + Configurable cfgEseHarmonic{"cfgEseHarmonic", 2, "Harmonic used only for ESE analysis (2 or 3)"}; Configurable cfgQvecDetName{"cfgQvecDetName", "FT0C", "The name of detector to be analyzed"}; Configurable cfgQvecRefAName{"cfgQvecRefAName", "TPCpos", "The name of detector for reference A"}; @@ -157,6 +162,7 @@ struct FlowEseTask { 15.370223f, 22.386745f, 28.322567f, 33.920535f, 39.547603f, 45.520061f, 52.251446f, 60.535625f, 72.660540f, 10.351031f, 15.094899f, 19.114684f, 22.915848f, 26.751668f, 30.837985f, 35.471803f, 41.208444f, 49.696178f}, "p10-p90 q2 cuts in each 10% centrality interval from 0 to 80%"}; + Configurable> cfgQ3PercentileCuts{"cfgQ3PercentileCuts", {}, "p10-p90 q3 cuts in each 10% centrality interval from 0 to 80%"}; ConfigurableAxis massAxis{"massAxis", {30, 1.1, 1.13}, "Invariant mass axis"}; ConfigurableAxis ptAxis{"ptAxis", {VARIABLE_WIDTH, 0.2, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0, 6.5, 8.0, 10.0, 100.0}, "Transverse momentum bins"}; @@ -165,22 +171,27 @@ struct FlowEseTask { ConfigurableAxis rapAxis{"rapAxis", {10, -0.5, 0.5}, "Rapidity axis"}; ConfigurableAxis qqAxis{"qqAxis", {100, -0.1, 0.1}, "qq axis"}; ConfigurableAxis multAxis{"multAxis", {300, 0, 2700}, "multiplicity"}; - ConfigurableAxis q2QaAxis{"q2QaAxis", {800, 0.0, 800.0}, "q2 axis for QA"}; + ConfigurableAxis q2QaAxis{"q2QaAxis", {8000, 0.0, 800.0}, "q2 axis for QA"}; static constexpr float MinAmplitudeThreshold = 1e-5f; static constexpr int ShiftLevel = 10; static constexpr int LambdaId = 3122; + static constexpr int SecondHarmonic = 2; + static constexpr int ThirdHarmonic = 3; static constexpr std::array CorrLevel = {2, 3, 4, 1}; static constexpr std::array CentBoundaries = {0.0f, 3.49f, 4.93f, 6.98f, 8.55f, 9.87f, 11.0f, 12.1f, 13.1f, 14.0f}; static constexpr std::array CentValues = {2.5f, 7.5f, 15.0f, 25.0f, 35.0f, 45.0f, 55.0f, 65.0f, 75.0f}; static constexpr float EtaAcceptance = 0.8f; static constexpr float CentUpperLimit = 80.0f; - static constexpr int NQ2CentBins = 8; - static constexpr int NQ2Groups = 10; - static constexpr int NQ2CutsPerCentBin = NQ2Groups - 1; - static constexpr float Q2CentMin = 0.0f; - static constexpr float Q2CentMax = 80.0f; - static constexpr float Q2CentBinWidth = 10.0f; + static constexpr int NEseCentBins = 8; + static constexpr int NEseGroups = 10; + static constexpr int NEseCutsPerCentBin = NEseGroups - 1; + static constexpr float EseCentMin = 0.0f; + static constexpr float EseCentMax = 80.0f; + static constexpr float EseCentBinWidth = 10.0f; + + TH2* histEventCountEseGroup = nullptr; + std::unordered_map histEse; EventPlaneHelper helperEP; @@ -234,60 +245,89 @@ struct FlowEseTask { return 0; } - int q2CentBin(float cent) const + int eseCentBin(float cent) const { - if (cent < Q2CentMin || cent >= Q2CentMax) { + if (cent < EseCentMin || cent >= EseCentMax) { return -1; } - return static_cast(cent / Q2CentBinWidth); + return static_cast(cent / EseCentBinWidth); } - const char* q2GroupSuffix(int group) const + std::string eseGroupSuffix(int group) const { - static constexpr std::array Q2GroupSuffixes = { - "q2p00_10", "q2p10_20", "q2p20_30", "q2p30_40", "q2p40_50", - "q2p50_60", "q2p60_70", "q2p70_80", "q2p80_90", "q2p90_100"}; - if (group < 0 || group >= NQ2Groups) { - return nullptr; + if (group < 0 || group >= NEseGroups) { + return {}; } - return Q2GroupSuffixes[group]; + return Form("q%dp%02d_%02d", cfgEseHarmonic.value, group * 10, (group + 1) * 10); + } + + const std::vector& esePercentileCuts() const + { + return cfgEseHarmonic.value == SecondHarmonic ? cfgQ2PercentileCuts.value : cfgQ3PercentileCuts.value; } - int q2Group(float cent, double q2) const + template + void addEseHistogram(const std::string& name, HistType type, const std::vector& axes) { - const int centBin = q2CentBin(cent); - if (centBin < 0 || !std::isfinite(q2)) { + histEse[name] = histos.add(name, "", type, axes).get(); + } + + template + T* getEseHistogram(const std::string& name) + { + const auto histogram = histEse.find(name); + if (histogram == histEse.end()) { + LOGF(fatal, "Could not find ESE histogram %s", name.c_str()); + } + return static_cast(histogram->second); + } + + int eseGroup(float cent, double qn) const + { + const int centBin = eseCentBin(cent); + if (centBin < 0 || !std::isfinite(qn)) { return -1; } - const int offset = centBin * NQ2CutsPerCentBin; - for (int iCut = 0; iCut < NQ2CutsPerCentBin; ++iCut) { - if (q2 < cfgQ2PercentileCuts->at(offset + iCut)) { + const auto& cuts = esePercentileCuts(); + const int offset = centBin * NEseCutsPerCentBin; + for (int iCut = 0; iCut < NEseCutsPerCentBin; ++iCut) { + if (qn < cuts.at(offset + iCut)) { return iCut; } } - return NQ2Groups - 1; + return NEseGroups - 1; } template - double getQ2(TCollision const& collision) + double getEseQ(TCollision const& collision) { + const int harmonicIndex = cfgEseHarmonic.value - SecondHarmonic; + if (collision.qvecFT0CReVec().size() <= static_cast(harmonicIndex) || collision.qvecFT0CImVec().size() <= static_cast(harmonicIndex)) { + LOGF(fatal, "FT0C Q-vector table does not contain harmonic %d", cfgEseHarmonic.value); + } + const double qx = collision.qvecFT0CReVec()[harmonicIndex]; + const double qy = collision.qvecFT0CImVec()[harmonicIndex]; if (cfgMultCor) { - return std::sqrt(collision.qvecFT0CReVec()[0] * collision.qvecFT0CReVec()[0] + collision.qvecFT0CImVec()[0] * collision.qvecFT0CImVec()[0]) * collision.sumAmplFT0C() / std::sqrt(collision.multFT0C()); + return std::sqrt(qx * qx + qy * qy) * collision.sumAmplFT0C() / std::sqrt(collision.multFT0C()); } - return std::sqrt(collision.qvecFT0CReVec()[0] * collision.qvecFT0CReVec()[0] + collision.qvecFT0CImVec()[0] * collision.qvecFT0CImVec()[0]) * std::sqrt(collision.sumAmplFT0C()); + return std::sqrt(qx * qx + qy * qy) * std::sqrt(collision.sumAmplFT0C()); } void init(o2::framework::InitContext&) { - if (cfgQ2PercentileCuts->size() != NQ2CentBins * NQ2CutsPerCentBin) { - LOGF(fatal, "cfgQ2PercentileCuts must contain %d values, got %d", NQ2CentBins * NQ2CutsPerCentBin, static_cast(cfgQ2PercentileCuts->size())); + if (cfgEseHarmonic.value != SecondHarmonic && cfgEseHarmonic.value != ThirdHarmonic) { + LOGF(fatal, "cfgEseHarmonic must be 2 or 3, got %d", cfgEseHarmonic.value); + } + const auto& eseCuts = esePercentileCuts(); + if (eseCuts.size() != NEseCentBins * NEseCutsPerCentBin) { + LOGF(fatal, "cfgQ%dPercentileCuts must contain %d values, got %d", cfgEseHarmonic.value, NEseCentBins * NEseCutsPerCentBin, static_cast(eseCuts.size())); } - for (int iCent = 0; iCent < NQ2CentBins; ++iCent) { - const int offset = iCent * NQ2CutsPerCentBin; - for (int iCut = 0; iCut < NQ2CutsPerCentBin; ++iCut) { - const float cut = cfgQ2PercentileCuts->at(offset + iCut); - if (!std::isfinite(cut) || (iCut > 0 && cut <= cfgQ2PercentileCuts->at(offset + iCut - 1))) { - LOGF(fatal, "Invalid q2 percentile cut at centrality bin %d, cut %d", iCent, iCut); + for (int iCent = 0; iCent < NEseCentBins; ++iCent) { + const int offset = iCent * NEseCutsPerCentBin; + for (int iCut = 0; iCut < NEseCutsPerCentBin; ++iCut) { + const float cut = eseCuts.at(offset + iCut); + if (!std::isfinite(cut) || (iCut > 0 && cut <= eseCuts.at(offset + iCut - 1))) { + LOGF(fatal, "Invalid q%d percentile cut at centrality bin %d, cut %d", cfgEseHarmonic.value, iCent, iCut); } } } @@ -302,19 +342,20 @@ struct FlowEseTask { AxisSpec shiftAxis = {10, 0, 10, "shift"}; AxisSpec basisAxis = {20, 0, 20, "basis"}; - AxisSpec q2GroupAxis = {10, 0.0, 10.0, "q_{2} percentile group"}; + AxisSpec eseGroupAxis = {10, 0.0, 10.0, Form("q_{%d} percentile group", cfgEseHarmonic.value)}; histos.add("histQvecCent", "", {HistType::kTH2F, {q2QaAxis, centQaAxis}}); - histos.add("histEventCountQ2Group", "", {HistType::kTH2F, {centQaAxis, q2GroupAxis}}); + histEventCountEseGroup = histos.add(Form("histEventCountQ%dGroup", cfgEseHarmonic.value), "", HistType::kTH2F, {centQaAxis, eseGroupAxis}).get(); histos.add(Form("histVertex"), "", {HistType::kTHnSparseF, {vertexAxis, vertexAxis, vertexAxis, centAxis}}); - for (int iGroup = 0; iGroup < NQ2Groups; ++iGroup) { - histos.add(Form("histV2_%s", q2GroupSuffix(iGroup)), "", {HistType::kTHnSparseF, {centAxis, ptAxis, cosAxis}}); - histos.add(Form("psi2/h_lambda_cos_%s", q2GroupSuffix(iGroup)), "", {HistType::kTHnSparseF, {massAxis, ptAxis, cosAxis, centAxis}}); - histos.add(Form("psi2/h_alambda_cos_%s", q2GroupSuffix(iGroup)), "", {HistType::kTHnSparseF, {massAxis, ptAxis, cosAxis, centAxis}}); - histos.add(Form("psi2/h_lambda_cos2_%s", q2GroupSuffix(iGroup)), "", {HistType::kTHnSparseF, {massAxis, ptAxis, cosAxis, centAxis}}); - histos.add(Form("psi2/h_alambda_cos2_%s", q2GroupSuffix(iGroup)), "", {HistType::kTHnSparseF, {massAxis, ptAxis, cosAxis, centAxis}}); - histos.add(Form("psi2/h_lambda_cossin_%s", q2GroupSuffix(iGroup)), "", {HistType::kTHnSparseF, {massAxis, ptAxis, cosAxis, centAxis}}); - histos.add(Form("psi2/h_alambda_cossin_%s", q2GroupSuffix(iGroup)), "", {HistType::kTHnSparseF, {massAxis, ptAxis, cosAxis, centAxis}}); + for (int iGroup = 0; iGroup < NEseGroups; ++iGroup) { + const auto suffix = eseGroupSuffix(iGroup); + addEseHistogram(Form("histV%d_%s", cfgEseHarmonic.value, suffix.c_str()), HistType::kTHnSparseF, {centAxis, ptAxis, cosAxis}); + addEseHistogram(Form("psi%d/h_lambda_cos_%s", cfgEseHarmonic.value, suffix.c_str()), HistType::kTHnSparseF, {massAxis, ptAxis, cosAxis, centAxis}); + addEseHistogram(Form("psi%d/h_alambda_cos_%s", cfgEseHarmonic.value, suffix.c_str()), HistType::kTHnSparseF, {massAxis, ptAxis, cosAxis, centAxis}); + addEseHistogram(Form("psi%d/h_lambda_cos2_%s", cfgEseHarmonic.value, suffix.c_str()), HistType::kTHnSparseF, {massAxis, ptAxis, cosAxis, centAxis}); + addEseHistogram(Form("psi%d/h_alambda_cos2_%s", cfgEseHarmonic.value, suffix.c_str()), HistType::kTHnSparseF, {massAxis, ptAxis, cosAxis, centAxis}); + addEseHistogram(Form("psi%d/h_lambda_cossin_%s", cfgEseHarmonic.value, suffix.c_str()), HistType::kTHnSparseF, {massAxis, ptAxis, cosAxis, centAxis}); + addEseHistogram(Form("psi%d/h_alambda_cossin_%s", cfgEseHarmonic.value, suffix.c_str()), HistType::kTHnSparseF, {massAxis, ptAxis, cosAxis, centAxis}); } histos.add("QA/CentDist", "", {HistType::kTH1F, {centQaAxis}}); histos.add("QA/PVzDist", "", {HistType::kTH1F, {pVzQaAxis}}); @@ -438,13 +479,6 @@ struct FlowEseTask { histos.add(Form("psi%d/QA/EPRes_Det_RefA", i), "", {HistType::kTH2F, {centQaAxis, cosAxis}}); histos.add(Form("psi%d/QA/EPRes_Det_RefB", i), "", {HistType::kTH2F, {centQaAxis, cosAxis}}); histos.add(Form("psi%d/QA/EPRes_RefA_RefB", i), "", {HistType::kTH2F, {centQaAxis, cosAxis}}); - if (i == CorrLevel[0]) { - for (int iGroup = 0; iGroup < NQ2Groups; ++iGroup) { - histos.add(Form("psi2/QA/EPRes_Det_RefA_%s", q2GroupSuffix(iGroup)), "", {HistType::kTH2F, {centQaAxis, cosAxis}}); - histos.add(Form("psi2/QA/EPRes_Det_RefB_%s", q2GroupSuffix(iGroup)), "", {HistType::kTH2F, {centQaAxis, cosAxis}}); - histos.add(Form("psi2/QA/EPRes_RefA_RefB_%s", q2GroupSuffix(iGroup)), "", {HistType::kTH2F, {centQaAxis, cosAxis}}); - } - } histos.add(Form("psi%d/QA/EP_FT0C_shifted", i), "", {HistType::kTH2F, {centQaAxis, epQaAxis}}); histos.add(Form("psi%d/QA/EP_FT0A_shifted", i), "", {HistType::kTH2F, {centQaAxis, epQaAxis}}); histos.add(Form("psi%d/QA/EP_FV0A_shifted", i), "", {HistType::kTH2F, {centQaAxis, epQaAxis}}); @@ -453,6 +487,12 @@ struct FlowEseTask { histos.add(Form("psi%d/QA/EPRes_FT0C_FV0A_shifted", i), "", {HistType::kTH2F, {centQaAxis, cosAxis}}); histos.add(Form("psi%d/QA/EPRes_FT0A_FV0A_shifted", i), "", {HistType::kTH2F, {centQaAxis, cosAxis}}); } + for (int iGroup = 0; iGroup < NEseGroups; ++iGroup) { + const auto suffix = eseGroupSuffix(iGroup); + addEseHistogram(Form("psi%d/QA/EPRes_Det_RefA_%s", cfgEseHarmonic.value, suffix.c_str()), HistType::kTH2F, {centQaAxis, cosAxis}); + addEseHistogram(Form("psi%d/QA/EPRes_Det_RefB_%s", cfgEseHarmonic.value, suffix.c_str()), HistType::kTH2F, {centQaAxis, cosAxis}); + addEseHistogram(Form("psi%d/QA/EPRes_RefA_RefB_%s", cfgEseHarmonic.value, suffix.c_str()), HistType::kTH2F, {centQaAxis, cosAxis}); + } } if (doprocessMcItsTpc) { @@ -469,6 +509,9 @@ struct FlowEseTask { for (auto i = 2; i < cfgnMods + 2; i++) { histos.add(Form("psi%d/ShiftFIT", i), "", kTProfile3D, {centQaAxis, basisAxis, shiftAxis}); } + if (cfgEseHarmonic.value >= cfgnMods.value + 2) { + histos.add(Form("psi%d/ShiftFIT", cfgEseHarmonic.value), "", kTProfile3D, {centQaAxis, basisAxis, shiftAxis}); + } } detId = getDetId(cfgQvecDetName); @@ -695,63 +738,6 @@ struct FlowEseTask { histos.fill(HIST("psi2/QA/EPRes_Det_RefA"), centrality, std::cos(std::atan2(collision.qvecIm()[qvecDetInd], collision.qvecRe()[qvecDetInd]) - std::atan2(collision.qvecIm()[qvecRefAInd], collision.qvecRe()[qvecRefAInd]))); histos.fill(HIST("psi2/QA/EPRes_Det_RefB"), centrality, std::cos(std::atan2(collision.qvecIm()[qvecDetInd], collision.qvecRe()[qvecDetInd]) - std::atan2(collision.qvecIm()[qvecRefBInd], collision.qvecRe()[qvecRefBInd]))); histos.fill(HIST("psi2/QA/EPRes_RefA_RefB"), centrality, std::cos(std::atan2(collision.qvecIm()[qvecRefAInd], collision.qvecRe()[qvecRefAInd]) - std::atan2(collision.qvecIm()[qvecRefBInd], collision.qvecRe()[qvecRefBInd]))); - const double epResDetRefA = std::cos(std::atan2(collision.qvecIm()[qvecDetInd], collision.qvecRe()[qvecDetInd]) - std::atan2(collision.qvecIm()[qvecRefAInd], collision.qvecRe()[qvecRefAInd])); - const double epResDetRefB = std::cos(std::atan2(collision.qvecIm()[qvecDetInd], collision.qvecRe()[qvecDetInd]) - std::atan2(collision.qvecIm()[qvecRefBInd], collision.qvecRe()[qvecRefBInd])); - const double epResRefARefB = std::cos(std::atan2(collision.qvecIm()[qvecRefAInd], collision.qvecRe()[qvecRefAInd]) - std::atan2(collision.qvecIm()[qvecRefBInd], collision.qvecRe()[qvecRefBInd])); - switch (q2Group(centrality, getQ2(collision))) { - case 0: - histos.fill(HIST("psi2/QA/EPRes_Det_RefA_q2p00_10"), centrality, epResDetRefA); - histos.fill(HIST("psi2/QA/EPRes_Det_RefB_q2p00_10"), centrality, epResDetRefB); - histos.fill(HIST("psi2/QA/EPRes_RefA_RefB_q2p00_10"), centrality, epResRefARefB); - break; - case 1: - histos.fill(HIST("psi2/QA/EPRes_Det_RefA_q2p10_20"), centrality, epResDetRefA); - histos.fill(HIST("psi2/QA/EPRes_Det_RefB_q2p10_20"), centrality, epResDetRefB); - histos.fill(HIST("psi2/QA/EPRes_RefA_RefB_q2p10_20"), centrality, epResRefARefB); - break; - case 2: - histos.fill(HIST("psi2/QA/EPRes_Det_RefA_q2p20_30"), centrality, epResDetRefA); - histos.fill(HIST("psi2/QA/EPRes_Det_RefB_q2p20_30"), centrality, epResDetRefB); - histos.fill(HIST("psi2/QA/EPRes_RefA_RefB_q2p20_30"), centrality, epResRefARefB); - break; - case 3: - histos.fill(HIST("psi2/QA/EPRes_Det_RefA_q2p30_40"), centrality, epResDetRefA); - histos.fill(HIST("psi2/QA/EPRes_Det_RefB_q2p30_40"), centrality, epResDetRefB); - histos.fill(HIST("psi2/QA/EPRes_RefA_RefB_q2p30_40"), centrality, epResRefARefB); - break; - case 4: - histos.fill(HIST("psi2/QA/EPRes_Det_RefA_q2p40_50"), centrality, epResDetRefA); - histos.fill(HIST("psi2/QA/EPRes_Det_RefB_q2p40_50"), centrality, epResDetRefB); - histos.fill(HIST("psi2/QA/EPRes_RefA_RefB_q2p40_50"), centrality, epResRefARefB); - break; - case 5: - histos.fill(HIST("psi2/QA/EPRes_Det_RefA_q2p50_60"), centrality, epResDetRefA); - histos.fill(HIST("psi2/QA/EPRes_Det_RefB_q2p50_60"), centrality, epResDetRefB); - histos.fill(HIST("psi2/QA/EPRes_RefA_RefB_q2p50_60"), centrality, epResRefARefB); - break; - case 6: - histos.fill(HIST("psi2/QA/EPRes_Det_RefA_q2p60_70"), centrality, epResDetRefA); - histos.fill(HIST("psi2/QA/EPRes_Det_RefB_q2p60_70"), centrality, epResDetRefB); - histos.fill(HIST("psi2/QA/EPRes_RefA_RefB_q2p60_70"), centrality, epResRefARefB); - break; - case 7: - histos.fill(HIST("psi2/QA/EPRes_Det_RefA_q2p70_80"), centrality, epResDetRefA); - histos.fill(HIST("psi2/QA/EPRes_Det_RefB_q2p70_80"), centrality, epResDetRefB); - histos.fill(HIST("psi2/QA/EPRes_RefA_RefB_q2p70_80"), centrality, epResRefARefB); - break; - case 8: - histos.fill(HIST("psi2/QA/EPRes_Det_RefA_q2p80_90"), centrality, epResDetRefA); - histos.fill(HIST("psi2/QA/EPRes_Det_RefB_q2p80_90"), centrality, epResDetRefB); - histos.fill(HIST("psi2/QA/EPRes_RefA_RefB_q2p80_90"), centrality, epResRefARefB); - break; - case 9: - histos.fill(HIST("psi2/QA/EPRes_Det_RefA_q2p90_100"), centrality, epResDetRefA); - histos.fill(HIST("psi2/QA/EPRes_Det_RefB_q2p90_100"), centrality, epResDetRefB); - histos.fill(HIST("psi2/QA/EPRes_RefA_RefB_q2p90_100"), centrality, epResRefARefB); - break; - default: - break; - } } else if (nmode == CorrLevel[1]) { histos.fill(HIST("psi3/QA/EP_Det"), centrality, std::atan2(collision.qvecIm()[qvecDetInd], collision.qvecRe()[qvecDetInd]) / static_cast(nmode)); histos.fill(HIST("psi3/QA/EP_RefA"), centrality, std::atan2(collision.qvecIm()[qvecRefAInd], collision.qvecRe()[qvecRefAInd]) / static_cast(nmode)); @@ -834,59 +820,76 @@ struct FlowEseTask { } } + template + void fillEseEPQA(TCollision const& collision, int eseGroupIndex) + { + if (eseGroupIndex < 0 || !cfgQAv0) { + return; + } + const int nmode = cfgEseHarmonic.value; + const int detIndex = detId * 4 + 3 + (nmode - 2) * cfgNQvec * 4; + const int refAIndex = refAId * 4 + 3 + (nmode - 2) * cfgNQvec * 4; + const int refBIndex = refBId * 4 + 3 + (nmode - 2) * cfgNQvec * 4; + if (collision.qvecAmp()[detId] < MinAmplitudeThreshold || collision.qvecAmp()[refAId] < MinAmplitudeThreshold || collision.qvecAmp()[refBId] < MinAmplitudeThreshold) { + return; + } + + const double detPhase = std::atan2(collision.qvecIm()[detIndex], collision.qvecRe()[detIndex]); + const double refAPhase = std::atan2(collision.qvecIm()[refAIndex], collision.qvecRe()[refAIndex]); + const double refBPhase = std::atan2(collision.qvecIm()[refBIndex], collision.qvecRe()[refBIndex]); + const auto suffix = eseGroupSuffix(eseGroupIndex); + auto* histEseEPResDetRefA = getEseHistogram(Form("psi%d/QA/EPRes_Det_RefA_%s", cfgEseHarmonic.value, suffix.c_str())); + auto* histEseEPResDetRefB = getEseHistogram(Form("psi%d/QA/EPRes_Det_RefB_%s", cfgEseHarmonic.value, suffix.c_str())); + auto* histEseEPResRefARefB = getEseHistogram(Form("psi%d/QA/EPRes_RefA_RefB_%s", cfgEseHarmonic.value, suffix.c_str())); + histEseEPResDetRefA->Fill(centrality, std::cos(detPhase - refAPhase)); + histEseEPResDetRefB->Fill(centrality, std::cos(detPhase - refBPhase)); + histEseEPResRefARefB->Fill(centrality, std::cos(refAPhase - refBPhase)); + } + template - void fillHistograms(TCollision const& collision, V0 const& V0s, TrackType const& track, int nmode) + void fillHistograms(TCollision const& collision, V0 const& V0s, TrackType const& track, int nmode, int eseGroupIndex = -1, bool fillRegular = true) { qvecDetInd = detId * 4 + 3 + (nmode - 2) * cfgNQvec * 4; qvecRefAInd = refAId * 4 + 3 + (nmode - 2) * cfgNQvec * 4; qvecRefBInd = refBId * 4 + 3 + (nmode - 2) * cfgNQvec * 4; - - const int q2PercentileGroup = nmode == CorrLevel[0] ? q2Group(centrality, getQ2(collision)) : -1; - - for (const auto& trk : track) { - if (!selectionTrack(trk)) { - continue; - } - if (nmode == CorrLevel[0]) { - const auto trackV2 = std::cos(static_cast(nmode) * (trk.phi() - helperEP.GetEventPlane(collision.qvecFT0CReVec()[0], collision.qvecFT0CImVec()[0], nmode))); - switch (q2PercentileGroup) { - case 0: - histos.fill(HIST("histV2_q2p00_10"), centrality, trk.pt(), trackV2); - break; - case 1: - histos.fill(HIST("histV2_q2p10_20"), centrality, trk.pt(), trackV2); - break; - case 2: - histos.fill(HIST("histV2_q2p20_30"), centrality, trk.pt(), trackV2); - break; - case 3: - histos.fill(HIST("histV2_q2p30_40"), centrality, trk.pt(), trackV2); - break; - case 4: - histos.fill(HIST("histV2_q2p40_50"), centrality, trk.pt(), trackV2); - break; - case 5: - histos.fill(HIST("histV2_q2p50_60"), centrality, trk.pt(), trackV2); - break; - case 6: - histos.fill(HIST("histV2_q2p60_70"), centrality, trk.pt(), trackV2); - break; - case 7: - histos.fill(HIST("histV2_q2p70_80"), centrality, trk.pt(), trackV2); - break; - case 8: - histos.fill(HIST("histV2_q2p80_90"), centrality, trk.pt(), trackV2); - break; - case 9: - histos.fill(HIST("histV2_q2p90_100"), centrality, trk.pt(), trackV2); - break; - default: - break; + const bool fillEse = eseGroupIndex >= 0 && nmode == cfgEseHarmonic.value; + + if (fillEse) { + const auto suffix = eseGroupSuffix(eseGroupIndex); + auto* histEseVn = getEseHistogram(Form("histV%d_%s", cfgEseHarmonic.value, suffix.c_str())); + const int harmonicIndex = nmode - 2; + const double esePlane = helperEP.GetEventPlane(collision.qvecFT0CReVec()[harmonicIndex], collision.qvecFT0CImVec()[harmonicIndex], nmode); + for (const auto& trk : track) { + if (!selectionTrack(trk)) { + continue; } + const double values[3] = {centrality, trk.pt(), std::cos(static_cast(nmode) * (trk.phi() - esePlane))}; + histEseVn->Fill(values); } } - histos.fill(HIST("histVertex"), collision.posX(), collision.posY(), collision.posZ(), collision.centFT0C()); + if (!fillRegular && !cfgFullCheck) { + return; + } + if (fillRegular) { + histos.fill(HIST("histVertex"), collision.posX(), collision.posY(), collision.posZ(), collision.centFT0C()); + } + + THnSparse* histEseLambdaCos = nullptr; + THnSparse* histEseAntiLambdaCos = nullptr; + THnSparse* histEseLambdaCos2 = nullptr; + THnSparse* histEseAntiLambdaCos2 = nullptr; + THnSparse* histEseLambdaCosSin = nullptr; + THnSparse* histEseAntiLambdaCosSin = nullptr; + if (fillEse && cfgFullCheck) { + const auto suffix = eseGroupSuffix(eseGroupIndex); + histEseLambdaCos = getEseHistogram(Form("psi%d/h_lambda_cos_%s", cfgEseHarmonic.value, suffix.c_str())); + histEseAntiLambdaCos = getEseHistogram(Form("psi%d/h_alambda_cos_%s", cfgEseHarmonic.value, suffix.c_str())); + histEseLambdaCos2 = getEseHistogram(Form("psi%d/h_lambda_cos2_%s", cfgEseHarmonic.value, suffix.c_str())); + histEseAntiLambdaCos2 = getEseHistogram(Form("psi%d/h_alambda_cos2_%s", cfgEseHarmonic.value, suffix.c_str())); + histEseLambdaCosSin = getEseHistogram(Form("psi%d/h_lambda_cossin_%s", cfgEseHarmonic.value, suffix.c_str())); + histEseAntiLambdaCosSin = getEseHistogram(Form("psi%d/h_alambda_cossin_%s", cfgEseHarmonic.value, suffix.c_str())); + } for (const auto& v0 : V0s) { auto postrack = v0.template posTrack_as(); @@ -898,7 +901,7 @@ struct FlowEseTask { double nTPCSigmaNegPr = negtrack.tpcNSigmaPr(); double nTPCSigmaPosPi = postrack.tpcNSigmaPi(); - if (cfgQAv0 && nmode == CorrLevel[0]) { + if (fillRegular && cfgQAv0 && nmode == CorrLevel[0]) { histos.fill(HIST("QA/nsigma_tpc_pt_ppr"), postrack.pt(), nTPCSigmaPosPr); histos.fill(HIST("QA/nsigma_tpc_pt_ppi"), postrack.pt(), nTPCSigmaPosPi); @@ -970,29 +973,50 @@ struct FlowEseTask { continue; } - if (lambdaTag) { - histos.fill(HIST("QA/ptspec_l"), v0.mLambda(), v0.pt(), centrality); - if (cfgEffCor) { - histos.fill(HIST("QA/ptspecCor_l"), v0.mLambda(), v0.pt(), centrality, - 1.0 / effMap->GetBinContent(effMap->GetXaxis()->FindBin(v0.pt()), effMap->GetYaxis()->FindBin(centrality))); + if (fillRegular) { + if (lambdaTag) { + histos.fill(HIST("QA/ptspec_l"), v0.mLambda(), v0.pt(), centrality); + if (cfgEffCor) { + histos.fill(HIST("QA/ptspecCor_l"), v0.mLambda(), v0.pt(), centrality, + 1.0 / effMap->GetBinContent(effMap->GetXaxis()->FindBin(v0.pt()), effMap->GetYaxis()->FindBin(centrality))); + } } - } - if (aLambdaTag) { - histos.fill(HIST("QA/ptspec_al"), v0.mAntiLambda(), v0.pt(), centrality); - if (cfgEffCor) { - histos.fill(HIST("QA/ptspecCor_al"), v0.mAntiLambda(), v0.pt(), centrality, - 1.0 / effMap->GetBinContent(effMap->GetXaxis()->FindBin(v0.pt()), effMap->GetYaxis()->FindBin(centrality))); + if (aLambdaTag) { + histos.fill(HIST("QA/ptspec_al"), v0.mAntiLambda(), v0.pt(), centrality); + if (cfgEffCor) { + histos.fill(HIST("QA/ptspecCor_al"), v0.mAntiLambda(), v0.pt(), centrality, + 1.0 / effMap->GetBinContent(effMap->GetXaxis()->FindBin(v0.pt()), effMap->GetYaxis()->FindBin(centrality))); + } } } double weight = 1.0; weight *= cfgEffCor ? 1.0 / effMap->GetBinContent(effMap->GetXaxis()->FindBin(v0.pt()), effMap->GetYaxis()->FindBin(centrality)) : 1.; weight *= cfgAccCor ? 1.0 / accMap->GetBinContent(accMap->GetXaxis()->FindBin(v0.pt()), accMap->GetYaxis()->FindBin(v0.yLambda())) : 1.; + if (fillEse && cfgFullCheck) { + const double mass = lambdaTag ? v0.mLambda() : v0.mAntiLambda(); + const double cosValues[4] = {mass, v0.pt(), angle * weight, centrality}; + const double cos2Values[4] = {mass, v0.pt(), angle * angle, centrality}; + const double cosSinValues[4] = {mass, v0.pt(), angle * std::sin(relphi) * weight, centrality}; + if (lambdaTag) { + histEseLambdaCos->Fill(cosValues); + histEseLambdaCos2->Fill(cos2Values); + histEseLambdaCosSin->Fill(cosSinValues); + } else { + histEseAntiLambdaCos->Fill(cosValues); + histEseAntiLambdaCos2->Fill(cos2Values); + histEseAntiLambdaCosSin->Fill(cosSinValues); + } + } + if (!fillRegular) { + continue; + } + double qvecMag = 1.0; if (cfgUSESP) qvecMag *= std::sqrt(std::pow(collision.qvecIm()[3 + (nmode - 2) * 28], 2) + std::pow(collision.qvecRe()[3 + (nmode - 2) * 28], 2)); - if (nmode == CorrLevel[0] && cfgFullCheck) { //////////// + if (fillRegular && nmode == CorrLevel[0] && cfgFullCheck) { //////////// if (lambdaTag) { histos.fill(HIST("psi2/h_lambda_cos"), v0.mLambda(), v0.pt(), angle * weight, centrality, relphi); histos.fill(HIST("psi2/h_lambda_cos2"), v0.mLambda(), v0.pt(), angle * angle, centrality, relphi); @@ -1001,60 +1025,6 @@ struct FlowEseTask { histos.fill(HIST("psi2/h_lambda_vncos"), v0.mLambda(), v0.pt(), qvecMag * std::cos(relphi) * weight, centrality); histos.fill(HIST("psi2/h_lambda_vnsin"), v0.mLambda(), v0.pt(), std::sin(relphi), centrality); - switch (q2PercentileGroup) { - case 0: - histos.fill(HIST("psi2/h_lambda_cos_q2p00_10"), v0.mLambda(), v0.pt(), angle * weight, centrality); - histos.fill(HIST("psi2/h_lambda_cos2_q2p00_10"), v0.mLambda(), v0.pt(), angle * angle, centrality); - histos.fill(HIST("psi2/h_lambda_cossin_q2p00_10"), v0.mLambda(), v0.pt(), angle * std::sin(relphi) * weight, centrality); - break; - case 1: - histos.fill(HIST("psi2/h_lambda_cos_q2p10_20"), v0.mLambda(), v0.pt(), angle * weight, centrality); - histos.fill(HIST("psi2/h_lambda_cos2_q2p10_20"), v0.mLambda(), v0.pt(), angle * angle, centrality); - histos.fill(HIST("psi2/h_lambda_cossin_q2p10_20"), v0.mLambda(), v0.pt(), angle * std::sin(relphi) * weight, centrality); - break; - case 2: - histos.fill(HIST("psi2/h_lambda_cos_q2p20_30"), v0.mLambda(), v0.pt(), angle * weight, centrality); - histos.fill(HIST("psi2/h_lambda_cos2_q2p20_30"), v0.mLambda(), v0.pt(), angle * angle, centrality); - histos.fill(HIST("psi2/h_lambda_cossin_q2p20_30"), v0.mLambda(), v0.pt(), angle * std::sin(relphi) * weight, centrality); - break; - case 3: - histos.fill(HIST("psi2/h_lambda_cos_q2p30_40"), v0.mLambda(), v0.pt(), angle * weight, centrality); - histos.fill(HIST("psi2/h_lambda_cos2_q2p30_40"), v0.mLambda(), v0.pt(), angle * angle, centrality); - histos.fill(HIST("psi2/h_lambda_cossin_q2p30_40"), v0.mLambda(), v0.pt(), angle * std::sin(relphi) * weight, centrality); - break; - case 4: - histos.fill(HIST("psi2/h_lambda_cos_q2p40_50"), v0.mLambda(), v0.pt(), angle * weight, centrality); - histos.fill(HIST("psi2/h_lambda_cos2_q2p40_50"), v0.mLambda(), v0.pt(), angle * angle, centrality); - histos.fill(HIST("psi2/h_lambda_cossin_q2p40_50"), v0.mLambda(), v0.pt(), angle * std::sin(relphi) * weight, centrality); - break; - case 5: - histos.fill(HIST("psi2/h_lambda_cos_q2p50_60"), v0.mLambda(), v0.pt(), angle * weight, centrality); - histos.fill(HIST("psi2/h_lambda_cos2_q2p50_60"), v0.mLambda(), v0.pt(), angle * angle, centrality); - histos.fill(HIST("psi2/h_lambda_cossin_q2p50_60"), v0.mLambda(), v0.pt(), angle * std::sin(relphi) * weight, centrality); - break; - case 6: - histos.fill(HIST("psi2/h_lambda_cos_q2p60_70"), v0.mLambda(), v0.pt(), angle * weight, centrality); - histos.fill(HIST("psi2/h_lambda_cos2_q2p60_70"), v0.mLambda(), v0.pt(), angle * angle, centrality); - histos.fill(HIST("psi2/h_lambda_cossin_q2p60_70"), v0.mLambda(), v0.pt(), angle * std::sin(relphi) * weight, centrality); - break; - case 7: - histos.fill(HIST("psi2/h_lambda_cos_q2p70_80"), v0.mLambda(), v0.pt(), angle * weight, centrality); - histos.fill(HIST("psi2/h_lambda_cos2_q2p70_80"), v0.mLambda(), v0.pt(), angle * angle, centrality); - histos.fill(HIST("psi2/h_lambda_cossin_q2p70_80"), v0.mLambda(), v0.pt(), angle * std::sin(relphi) * weight, centrality); - break; - case 8: - histos.fill(HIST("psi2/h_lambda_cos_q2p80_90"), v0.mLambda(), v0.pt(), angle * weight, centrality); - histos.fill(HIST("psi2/h_lambda_cos2_q2p80_90"), v0.mLambda(), v0.pt(), angle * angle, centrality); - histos.fill(HIST("psi2/h_lambda_cossin_q2p80_90"), v0.mLambda(), v0.pt(), angle * std::sin(relphi) * weight, centrality); - break; - case 9: - histos.fill(HIST("psi2/h_lambda_cos_q2p90_100"), v0.mLambda(), v0.pt(), angle * weight, centrality); - histos.fill(HIST("psi2/h_lambda_cos2_q2p90_100"), v0.mLambda(), v0.pt(), angle * angle, centrality); - histos.fill(HIST("psi2/h_lambda_cossin_q2p90_100"), v0.mLambda(), v0.pt(), angle * std::sin(relphi) * weight, centrality); - break; - default: - break; - } if (cfgRapidityDep) { histos.fill(HIST("psi2/h_lambda_cos2_rap"), v0.mLambda(), v0.pt(), angle * angle, centrality, v0.yLambda(), weight); } @@ -1105,60 +1075,6 @@ struct FlowEseTask { histos.fill(HIST("psi2/h_alambda_vncos"), v0.mAntiLambda(), v0.pt(), qvecMag * std::cos(relphi) * weight, centrality); histos.fill(HIST("psi2/h_alambda_vnsin"), v0.mAntiLambda(), v0.pt(), std::sin(relphi), centrality); - switch (q2PercentileGroup) { - case 0: - histos.fill(HIST("psi2/h_alambda_cos_q2p00_10"), v0.mAntiLambda(), v0.pt(), angle * weight, centrality); - histos.fill(HIST("psi2/h_alambda_cos2_q2p00_10"), v0.mAntiLambda(), v0.pt(), angle * angle, centrality); - histos.fill(HIST("psi2/h_alambda_cossin_q2p00_10"), v0.mAntiLambda(), v0.pt(), angle * std::sin(relphi) * weight, centrality); - break; - case 1: - histos.fill(HIST("psi2/h_alambda_cos_q2p10_20"), v0.mAntiLambda(), v0.pt(), angle * weight, centrality); - histos.fill(HIST("psi2/h_alambda_cos2_q2p10_20"), v0.mAntiLambda(), v0.pt(), angle * angle, centrality); - histos.fill(HIST("psi2/h_alambda_cossin_q2p10_20"), v0.mAntiLambda(), v0.pt(), angle * std::sin(relphi) * weight, centrality); - break; - case 2: - histos.fill(HIST("psi2/h_alambda_cos_q2p20_30"), v0.mAntiLambda(), v0.pt(), angle * weight, centrality); - histos.fill(HIST("psi2/h_alambda_cos2_q2p20_30"), v0.mAntiLambda(), v0.pt(), angle * angle, centrality); - histos.fill(HIST("psi2/h_alambda_cossin_q2p20_30"), v0.mAntiLambda(), v0.pt(), angle * std::sin(relphi) * weight, centrality); - break; - case 3: - histos.fill(HIST("psi2/h_alambda_cos_q2p30_40"), v0.mAntiLambda(), v0.pt(), angle * weight, centrality); - histos.fill(HIST("psi2/h_alambda_cos2_q2p30_40"), v0.mAntiLambda(), v0.pt(), angle * angle, centrality); - histos.fill(HIST("psi2/h_alambda_cossin_q2p30_40"), v0.mAntiLambda(), v0.pt(), angle * std::sin(relphi) * weight, centrality); - break; - case 4: - histos.fill(HIST("psi2/h_alambda_cos_q2p40_50"), v0.mAntiLambda(), v0.pt(), angle * weight, centrality); - histos.fill(HIST("psi2/h_alambda_cos2_q2p40_50"), v0.mAntiLambda(), v0.pt(), angle * angle, centrality); - histos.fill(HIST("psi2/h_alambda_cossin_q2p40_50"), v0.mAntiLambda(), v0.pt(), angle * std::sin(relphi) * weight, centrality); - break; - case 5: - histos.fill(HIST("psi2/h_alambda_cos_q2p50_60"), v0.mAntiLambda(), v0.pt(), angle * weight, centrality); - histos.fill(HIST("psi2/h_alambda_cos2_q2p50_60"), v0.mAntiLambda(), v0.pt(), angle * angle, centrality); - histos.fill(HIST("psi2/h_alambda_cossin_q2p50_60"), v0.mAntiLambda(), v0.pt(), angle * std::sin(relphi) * weight, centrality); - break; - case 6: - histos.fill(HIST("psi2/h_alambda_cos_q2p60_70"), v0.mAntiLambda(), v0.pt(), angle * weight, centrality); - histos.fill(HIST("psi2/h_alambda_cos2_q2p60_70"), v0.mAntiLambda(), v0.pt(), angle * angle, centrality); - histos.fill(HIST("psi2/h_alambda_cossin_q2p60_70"), v0.mAntiLambda(), v0.pt(), angle * std::sin(relphi) * weight, centrality); - break; - case 7: - histos.fill(HIST("psi2/h_alambda_cos_q2p70_80"), v0.mAntiLambda(), v0.pt(), angle * weight, centrality); - histos.fill(HIST("psi2/h_alambda_cos2_q2p70_80"), v0.mAntiLambda(), v0.pt(), angle * angle, centrality); - histos.fill(HIST("psi2/h_alambda_cossin_q2p70_80"), v0.mAntiLambda(), v0.pt(), angle * std::sin(relphi) * weight, centrality); - break; - case 8: - histos.fill(HIST("psi2/h_alambda_cos_q2p80_90"), v0.mAntiLambda(), v0.pt(), angle * weight, centrality); - histos.fill(HIST("psi2/h_alambda_cos2_q2p80_90"), v0.mAntiLambda(), v0.pt(), angle * angle, centrality); - histos.fill(HIST("psi2/h_alambda_cossin_q2p80_90"), v0.mAntiLambda(), v0.pt(), angle * std::sin(relphi) * weight, centrality); - break; - case 9: - histos.fill(HIST("psi2/h_alambda_cos_q2p90_100"), v0.mAntiLambda(), v0.pt(), angle * weight, centrality); - histos.fill(HIST("psi2/h_alambda_cos2_q2p90_100"), v0.mAntiLambda(), v0.pt(), angle * angle, centrality); - histos.fill(HIST("psi2/h_alambda_cossin_q2p90_100"), v0.mAntiLambda(), v0.pt(), angle * std::sin(relphi) * weight, centrality); - break; - default: - break; - } if (cfgRapidityDep) { histos.fill(HIST("psi2/h_alambda_cos2_rap"), v0.mAntiLambda(), v0.pt(), angle * angle, centrality, v0.yLambda(), weight); } @@ -1201,7 +1117,7 @@ struct FlowEseTask { histos.fill(HIST("psi2/QA/sinPhi_al"), v0.mAntiLambda(), v0.pt(), std::sin(v0.phi() * 2.0), centrality); } } - } else if (nmode == CorrLevel[1]) { + } else if (fillRegular && nmode == CorrLevel[1]) { if (lambdaTag) { histos.fill(HIST("psi3/h_lambda_cos"), v0.mLambda(), v0.pt(), angle * weight, centrality, relphi); histos.fill(HIST("psi3/h_lambda_cos2"), v0.mLambda(), v0.pt(), angle * angle, centrality, relphi); @@ -1232,7 +1148,7 @@ struct FlowEseTask { histos.fill(HIST("psi3/h_alambda_coscos"), v0.mAntiLambda(), v0.pt(), angle * std::cos(relphi), centrality, weight); } } - } else if (nmode == CorrLevel[2]) { + } else if (fillRegular && nmode == CorrLevel[2]) { if (lambdaTag) { histos.fill(HIST("psi4/h_lambda_cos"), v0.mLambda(), v0.pt(), angle * weight, centrality, relphi); histos.fill(HIST("psi4/h_lambda_cos2"), v0.mLambda(), v0.pt(), angle * angle, centrality, relphi); @@ -1281,10 +1197,11 @@ struct FlowEseTask { } histos.fill(HIST("QA/CentDist"), centrality, 1.0); histos.fill(HIST("QA/PVzDist"), collision.posZ(), 1.0); - histos.fill(HIST("histQvecCent"), getQ2(collision), centrality); - const int q2PercentileGroup = q2Group(centrality, getQ2(collision)); - if (q2PercentileGroup >= 0) { - histos.fill(HIST("histEventCountQ2Group"), centrality, q2PercentileGroup + 0.5); + const double eseQ = getEseQ(collision); + histos.fill(HIST("histQvecCent"), eseQ, centrality); + const int eseGroupIndex = eseGroup(centrality, eseQ); + if (eseGroupIndex >= 0) { + histEventCountEseGroup->Fill(centrality, eseGroupIndex + 0.5); } if (cfgShiftCorr) { @@ -1292,7 +1209,8 @@ struct FlowEseTask { currentRunNumber = bc.runNumber(); if (currentRunNumber != lastRunNumber) { shiftprofile.clear(); - for (int i = 2; i < cfgnMods + 2; i++) { + const int maxShiftHarmonic = std::max(cfgnMods.value + 1, cfgEseHarmonic.value); + for (int i = 2; i <= maxShiftHarmonic; i++) { fullCCDBShiftCorrPath = cfgShiftPath; fullCCDBShiftCorrPath += "/v"; fullCCDBShiftCorrPath += std::to_string(i); @@ -1309,6 +1227,13 @@ struct FlowEseTask { if (cfgAccCor) { accMap = ccdb->getForTimeStamp(cfgAccCorPath.value, bc.timestamp()); } + fillEseEPQA(collision, eseGroupIndex); + if (cfgShiftCorrDef && cfgEseHarmonic.value >= cfgnMods.value + 2) { + fillShiftCorrection(collision, cfgEseHarmonic.value); + } + if (eseGroupIndex >= 0 && cfgEseHarmonic.value >= cfgnMods.value + 2) { + fillHistograms(collision, V0s, tracks, cfgEseHarmonic.value, eseGroupIndex, false); + } for (int i = 2; i < cfgnMods + 2; i++) { if (cfgShiftCorrDef) { fillShiftCorrection(collision, i); @@ -1316,7 +1241,7 @@ struct FlowEseTask { if (cfgQAv0) { fillEPQA(collision, i); } - fillHistograms(collision, V0s, tracks, i); + fillHistograms(collision, V0s, tracks, i, i == cfgEseHarmonic.value ? eseGroupIndex : -1); } // FIXME: need to fill different histograms for different harmonic } PROCESS_SWITCH(FlowEseTask, processData, "Process Event for data", true); From d8c5b05debdf764bd5d409442856b48ca601a8ce Mon Sep 17 00:00:00 2001 From: Shiqi Wang <3326561594@qq.com> Date: Tue, 18 Aug 2026 19:28:24 +0800 Subject: [PATCH 2/3] Add hist for 3rd order ESE calculation --- PWGCF/Flow/Tasks/flowEseTask.cxx | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/PWGCF/Flow/Tasks/flowEseTask.cxx b/PWGCF/Flow/Tasks/flowEseTask.cxx index 33379961c96..77bab7594d7 100644 --- a/PWGCF/Flow/Tasks/flowEseTask.cxx +++ b/PWGCF/Flow/Tasks/flowEseTask.cxx @@ -278,6 +278,7 @@ struct FlowEseTask { const auto histogram = histEse.find(name); if (histogram == histEse.end()) { LOGF(fatal, "Could not find ESE histogram %s", name.c_str()); + return nullptr; } return static_cast(histogram->second); } @@ -1016,7 +1017,7 @@ struct FlowEseTask { if (cfgUSESP) qvecMag *= std::sqrt(std::pow(collision.qvecIm()[3 + (nmode - 2) * 28], 2) + std::pow(collision.qvecRe()[3 + (nmode - 2) * 28], 2)); - if (fillRegular && nmode == CorrLevel[0] && cfgFullCheck) { //////////// + if (nmode == CorrLevel[0] && cfgFullCheck) { //////////// if (lambdaTag) { histos.fill(HIST("psi2/h_lambda_cos"), v0.mLambda(), v0.pt(), angle * weight, centrality, relphi); histos.fill(HIST("psi2/h_lambda_cos2"), v0.mLambda(), v0.pt(), angle * angle, centrality, relphi); @@ -1117,7 +1118,7 @@ struct FlowEseTask { histos.fill(HIST("psi2/QA/sinPhi_al"), v0.mAntiLambda(), v0.pt(), std::sin(v0.phi() * 2.0), centrality); } } - } else if (fillRegular && nmode == CorrLevel[1]) { + } else if (nmode == CorrLevel[1]) { if (lambdaTag) { histos.fill(HIST("psi3/h_lambda_cos"), v0.mLambda(), v0.pt(), angle * weight, centrality, relphi); histos.fill(HIST("psi3/h_lambda_cos2"), v0.mLambda(), v0.pt(), angle * angle, centrality, relphi); @@ -1148,7 +1149,7 @@ struct FlowEseTask { histos.fill(HIST("psi3/h_alambda_coscos"), v0.mAntiLambda(), v0.pt(), angle * std::cos(relphi), centrality, weight); } } - } else if (fillRegular && nmode == CorrLevel[2]) { + } else if (nmode == CorrLevel[2]) { if (lambdaTag) { histos.fill(HIST("psi4/h_lambda_cos"), v0.mLambda(), v0.pt(), angle * weight, centrality, relphi); histos.fill(HIST("psi4/h_lambda_cos2"), v0.mLambda(), v0.pt(), angle * angle, centrality, relphi); From dbc63b8a8fc4b72abd590f7cbcd3a0ba63a5fda9 Mon Sep 17 00:00:00 2001 From: Shiqi Wang <3326561594@qq.com> Date: Wed, 19 Aug 2026 09:24:04 +0800 Subject: [PATCH 3/3] Add 3rd order ESE calculation and fix CI warnings --- PWGCF/Flow/Tasks/flowEseTask.cxx | 128 ++++++++++++++++++++----------- 1 file changed, 82 insertions(+), 46 deletions(-) diff --git a/PWGCF/Flow/Tasks/flowEseTask.cxx b/PWGCF/Flow/Tasks/flowEseTask.cxx index 77bab7594d7..20553572aa7 100644 --- a/PWGCF/Flow/Tasks/flowEseTask.cxx +++ b/PWGCF/Flow/Tasks/flowEseTask.cxx @@ -46,6 +46,8 @@ #include // IWYU pragma: keep (do not replace with Math/Vector4Dfwd.h) #include #include +#include +#include #include #include #include @@ -303,7 +305,7 @@ struct FlowEseTask { double getEseQ(TCollision const& collision) { const int harmonicIndex = cfgEseHarmonic.value - SecondHarmonic; - if (collision.qvecFT0CReVec().size() <= static_cast(harmonicIndex) || collision.qvecFT0CImVec().size() <= static_cast(harmonicIndex)) { + if (collision.qvecFT0CReVec().size() <= static_cast(harmonicIndex) || collision.qvecFT0CImVec().size() <= static_cast(harmonicIndex)) { LOGF(fatal, "FT0C Q-vector table does not contain harmonic %d", cfgEseHarmonic.value); } const double qx = collision.qvecFT0CReVec()[harmonicIndex]; @@ -585,31 +587,42 @@ struct FlowEseTask { template bool selectionV0(TCollision const& collision, V0 const& candidate, int lambdaTag) { - if (candidate.v0radius() < cfgv0radiusMin) + if (candidate.v0radius() < cfgv0radiusMin) { return false; + } if (lambdaTag) { - if (std::abs(candidate.dcapostopv()) < cfgDCAPrToPVMin) + if (std::abs(candidate.dcapostopv()) < cfgDCAPrToPVMin) { return false; - if (std::abs(candidate.dcanegtopv()) < cfgDCAPiToPVMin) + } + if (std::abs(candidate.dcanegtopv()) < cfgDCAPiToPVMin) { return false; + } } else { - if (std::abs(candidate.dcapostopv()) < cfgDCAPiToPVMin) + if (std::abs(candidate.dcapostopv()) < cfgDCAPiToPVMin) { return false; - if (std::abs(candidate.dcanegtopv()) < cfgDCAPrToPVMin) + } + if (std::abs(candidate.dcanegtopv()) < cfgDCAPrToPVMin) { return false; + } } - if (candidate.v0cosPA() < cfgv0CosPA) + if (candidate.v0cosPA() < cfgv0CosPA) { return false; - if (std::abs(candidate.dcaV0daughters()) > cfgDCAV0Dau) + } + if (std::abs(candidate.dcaV0daughters()) > cfgDCAV0Dau) { return false; - if (candidate.pt() < cfgV0PtMin) + } + if (candidate.pt() < cfgV0PtMin) { return false; - if (candidate.yLambda() < cfgV0EtaMin) + } + if (candidate.yLambda() < cfgV0EtaMin) { return false; - if (candidate.yLambda() > cfgV0EtaMax) + } + if (candidate.yLambda() > cfgV0EtaMax) { return false; - if (candidate.distovertotmom(collision.posX(), collision.posY(), collision.posZ()) * massLambda > cfgV0LifeTime) + } + if (candidate.distovertotmom(collision.posX(), collision.posY(), collision.posZ()) * massLambda > cfgV0LifeTime) { return false; + } return true; } @@ -617,20 +630,27 @@ struct FlowEseTask { template bool isSelectedV0Daughter(T const& track, int pid) // pid 0: proton, pid 1: pion { - if (track.tpcNClsFound() < cfgDaughTPCnclsMin) + if (track.tpcNClsFound() < cfgDaughTPCnclsMin) { return false; - if (pid == 0 && std::abs(track.tpcNSigmaPr()) > cfgDaughPIDCutsTPCPr) + } + if (pid == 0 && std::abs(track.tpcNSigmaPr()) > cfgDaughPIDCutsTPCPr) { return false; - if (pid == 1 && std::abs(track.tpcNSigmaPi()) > cfgDaughPIDCutsTPCPi) + } + if (pid == 1 && std::abs(track.tpcNSigmaPi()) > cfgDaughPIDCutsTPCPi) { return false; - if (track.eta() > cfgDaughEtaMax) + } + if (track.eta() > cfgDaughEtaMax) { return false; - if (track.eta() < cfgDaughEtaMin) + } + if (track.eta() < cfgDaughEtaMin) { return false; - if (pid == 0 && track.pt() < cfgDaughPrPt) + } + if (pid == 0 && track.pt() < cfgDaughPrPt) { return false; - if (pid == 1 && track.pt() < cfgDaughPiPt) + } + if (pid == 1 && track.pt() < cfgDaughPiPt) { return false; + } return true; } @@ -652,24 +672,33 @@ struct FlowEseTask { template bool selectionTrack(TrackType const& track) { - if (track.pt() < cfgMinPt) + if (track.pt() < cfgMinPt) { return false; - if (std::abs(track.eta()) > cfgMaxEta) + } + if (std::abs(track.eta()) > cfgMaxEta) { return false; - if (!track.passedITSNCls()) + } + if (!track.passedITSNCls()) { return false; - if (!track.passedITSChi2NDF()) + } + if (!track.passedITSChi2NDF()) { return false; - if (!track.passedITSHits()) + } + if (!track.passedITSHits()) { return false; - if (!track.passedTPCCrossedRowsOverNCls()) + } + if (!track.passedTPCCrossedRowsOverNCls()) { return false; - if (!track.passedTPCChi2NDF()) + } + if (!track.passedTPCChi2NDF()) { return false; - if (!track.passedDCAxy()) + } + if (!track.passedDCAxy()) { return false; - if (!track.passedDCAz()) + } + if (!track.passedDCAz()) { return false; + } return true; } @@ -720,8 +749,9 @@ struct FlowEseTask { qvecRefAInd = refAId * 4 + 3 + (nmode - 2) * cfgNQvec * 4; qvecRefBInd = refBId * 4 + 3 + (nmode - 2) * cfgNQvec * 4; - if (collision.qvecAmp()[detId] < MinAmplitudeThreshold || collision.qvecAmp()[refAId] < MinAmplitudeThreshold || collision.qvecAmp()[refBId] < MinAmplitudeThreshold) + if (collision.qvecAmp()[detId] < MinAmplitudeThreshold || collision.qvecAmp()[refAId] < MinAmplitudeThreshold || collision.qvecAmp()[refBId] < MinAmplitudeThreshold) { return; + } if (nmode == CorrLevel[0]) { histos.fill(HIST("psi2/QA/EP_Det"), centrality, std::atan2(collision.qvecIm()[qvecDetInd], collision.qvecRe()[qvecDetInd]) / static_cast(nmode)); @@ -864,8 +894,8 @@ struct FlowEseTask { if (!selectionTrack(trk)) { continue; } - const double values[3] = {centrality, trk.pt(), std::cos(static_cast(nmode) * (trk.phi() - esePlane))}; - histEseVn->Fill(values); + const std::array values = {centrality, trk.pt(), std::cos(static_cast(nmode) * (trk.phi() - esePlane))}; + histEseVn->Fill(values.data()); } } @@ -920,11 +950,13 @@ struct FlowEseTask { aLambdaTag = 1; } - if (lambdaTag == aLambdaTag) + if (lambdaTag == aLambdaTag) { continue; + } - if (!selectionV0(collision, v0, lambdaTag)) + if (!selectionV0(collision, v0, lambdaTag)) { continue; + } if (lambdaTag) { protonVec = ROOT::Math::PxPyPzMVector(v0.pxpos(), v0.pypos(), v0.pzpos(), massPr); @@ -996,17 +1028,17 @@ struct FlowEseTask { if (fillEse && cfgFullCheck) { const double mass = lambdaTag ? v0.mLambda() : v0.mAntiLambda(); - const double cosValues[4] = {mass, v0.pt(), angle * weight, centrality}; - const double cos2Values[4] = {mass, v0.pt(), angle * angle, centrality}; - const double cosSinValues[4] = {mass, v0.pt(), angle * std::sin(relphi) * weight, centrality}; + const std::array cosValues = {mass, v0.pt(), angle * weight, centrality}; + const std::array cos2Values = {mass, v0.pt(), angle * angle, centrality}; + const std::array cosSinValues = {mass, v0.pt(), angle * std::sin(relphi) * weight, centrality}; if (lambdaTag) { - histEseLambdaCos->Fill(cosValues); - histEseLambdaCos2->Fill(cos2Values); - histEseLambdaCosSin->Fill(cosSinValues); + histEseLambdaCos->Fill(cosValues.data()); + histEseLambdaCos2->Fill(cos2Values.data()); + histEseLambdaCosSin->Fill(cosSinValues.data()); } else { - histEseAntiLambdaCos->Fill(cosValues); - histEseAntiLambdaCos2->Fill(cos2Values); - histEseAntiLambdaCosSin->Fill(cosSinValues); + histEseAntiLambdaCos->Fill(cosValues.data()); + histEseAntiLambdaCos2->Fill(cos2Values.data()); + histEseAntiLambdaCosSin->Fill(cosSinValues.data()); } } if (!fillRegular) { @@ -1014,8 +1046,9 @@ struct FlowEseTask { } double qvecMag = 1.0; - if (cfgUSESP) + if (cfgUSESP) { qvecMag *= std::sqrt(std::pow(collision.qvecIm()[3 + (nmode - 2) * 28], 2) + std::pow(collision.qvecRe()[3 + (nmode - 2) * 28], 2)); + } if (nmode == CorrLevel[0] && cfgFullCheck) { //////////// if (lambdaTag) { @@ -1291,12 +1324,15 @@ struct FlowEseTask { float deltaPhi = mcParticle.phi() - mcCollision.eventPlaneAngle(); // focus on bulk: e, mu, pi, k, p int pdgCode = std::abs(mcParticle.pdgCode()); - if (pdgCode != LambdaId) + if (pdgCode != LambdaId) { continue; - if (!mcParticle.isPhysicalPrimary()) + } + if (!mcParticle.isPhysicalPrimary()) { continue; - if (std::abs(mcParticle.eta()) > EtaAcceptance) // main acceptance + } + if (std::abs(mcParticle.eta()) > EtaAcceptance) { // main acceptance continue; + } histos.fill(HIST("hSparseMCGenWeight"), centclass, RecoDecay::constrainAngle(deltaPhi, 0, 2), std::pow(std::cos(2.0 * RecoDecay::constrainAngle(deltaPhi, 0, 2)), 2.0), mcParticle.pt(), mcParticle.eta()); nCh++; bool validGlobal = false;