From a635c8b7e2adba91ea6246a6afb4f41ac4c4f8c0 Mon Sep 17 00:00:00 2001 From: Luca Barioglio Date: Wed, 20 May 2026 15:21:21 +0200 Subject: [PATCH 1/3] Add quadratic response option --- PWGLF/TableProducer/QC/flowQC.cxx | 86 ++++++++++++++----------------- 1 file changed, 40 insertions(+), 46 deletions(-) diff --git a/PWGLF/TableProducer/QC/flowQC.cxx b/PWGLF/TableProducer/QC/flowQC.cxx index 2d590e0df69..5b72d66d136 100644 --- a/PWGLF/TableProducer/QC/flowQC.cxx +++ b/PWGLF/TableProducer/QC/flowQC.cxx @@ -89,8 +89,6 @@ enum methods { }; static const std::vector suffixes = {"EP", "Qvec"}; -std::shared_ptr hQxQy[kNmethods][kNqVecDetectors]; -std::shared_ptr hNormQxQy[kNmethods][kNqVecDetectors]; std::shared_ptr hPsi[kNmethods][kNqVecDetectors]; std::shared_ptr hDeltaPsi[kNmethods][kNqVecDetectors][kNqVecDetectors]; std::shared_ptr hScalarProduct[kNmethods][kNqVecDetectors][kNqVecDetectors]; @@ -122,6 +120,7 @@ struct flowQC { float mBz = 0.f; Configurable cfgHarmonic{"cfgHarmonic", 2.f, "Harmonics for flow analysis"}; + Configurable cfgQuadraticResponse{"cfgQuadraticResponse", false, "Use quadratic response for Q-vector quantities"}; // Flow analysis using CollWithEPandQvec = soa::Joinadd(Form("hQxQy_%s_%s", qVecDetectorNames[iQvecDet].c_str(), suffixes[iMethod].c_str()), "", HistType::kTH3F, {centAxis, QxAxis, QyAxis}); - hNormQxQy[iMethod][iQvecDet] = registry->add(Form("hNormQxQy_%s_%s", qVecDetectorNames[iQvecDet].c_str(), suffixes[iMethod].c_str()), "", HistType::kTH3F, {centAxis, NormQxAxis, NormQyAxis}); hPsi[iMethod][iQvecDet] = registry->add(Form("hPsi_%s_%s", qVecDetectorNames[iQvecDet].c_str(), suffixes[iMethod].c_str()), "", HistType::kTH2F, {centAxis, psiAxis}); for (int jQvecDet = iQvecDet + 1; jQvecDet < qVecDetectors::kNqVecDetectors; jQvecDet++) { @@ -215,12 +208,12 @@ struct flowQC { hDeltaPsi[iMethod][iQvecDet][jQvecDet] = registry->add(Form("hDeltaPsi_%s_%s_%s", qVecDetectorNames[iQvecDet].c_str(), qVecDetectorNames[jQvecDet].c_str(), suffixes[iMethod].c_str()), "", HistType::kTH2F, {centAxis, {cfgDeltaPhiBins, Form("#psi_{%s} - #psi_{%s}", qVecDetectorNames[iQvecDet].c_str(), qVecDetectorNames[jQvecDet].c_str())}}); // Scalar-product histograms - auto spLabel = Form("#vec{Q}_{%.0f}^{%s} #upoint #vec{Q}_{%.0f}^{%s}", cfgHarmonic.value, qVecDetectorNames[iQvecDet].c_str(), cfgHarmonic.value, qVecDetectorNames[jQvecDet].c_str()); + auto spLabel = Form("#vec{%s}_{%.0f}^{%s} #upoint #vec{%s}_{%.0f}^{%s}", qLabel, cfgHarmonic.value, qVecDetectorNames[iQvecDet].c_str(), qLabel, cfgHarmonic.value, qVecDetectorNames[jQvecDet].c_str()); hScalarProduct[iMethod][iQvecDet][jQvecDet] = registry->add(Form("hScalarProduct_%s_%s_%s", qVecDetectorNames[iQvecDet].c_str(), qVecDetectorNames[jQvecDet].c_str(), suffixes[iMethod].c_str()), "", HistType::kTH2F, {centAxis, {cfgQvecBins, spLabel}}); // Normalised scalar-product histograms - auto normSpLabel = Form("#frac{#vec{Q}_{%.0f}^{%s} #upoint #vec{Q}_{%.0f}^{%s}}{||#vec{Q}_{%.0f}^{%s}|| ||#vec{Q}_{%.0f}^{%s}||}", cfgHarmonic.value, qVecDetectorNames[iQvecDet].c_str(), cfgHarmonic.value, qVecDetectorNames[jQvecDet].c_str(), cfgHarmonic.value, qVecDetectorNames[iQvecDet].c_str(), cfgHarmonic.value, qVecDetectorNames[jQvecDet].c_str()); + auto normSpLabel = Form("#frac{#vec{%s}_{%.0f}^{%s} #upoint #vec{%s}_{%.0f}^{%s}}{||#vec{%s}_{%.0f}^{%s}|| ||#vec{%s}_{%.0f}^{%s}||}", qLabel, cfgHarmonic.value, qVecDetectorNames[iQvecDet].c_str(), qLabel, cfgHarmonic.value, qVecDetectorNames[jQvecDet].c_str(), qLabel, cfgHarmonic.value, qVecDetectorNames[iQvecDet].c_str(), qLabel, cfgHarmonic.value, qVecDetectorNames[jQvecDet].c_str()); hNormalisedScalarProduct[iMethod][iQvecDet][jQvecDet] = registry->add(Form("hNormalisedScalarProduct_%s_%s_%s", qVecDetectorNames[iQvecDet].c_str(), qVecDetectorNames[jQvecDet].c_str(), suffixes[iMethod].c_str()), "", HistType::kTH2F, {centAxis, {cfgQvecBins, normSpLabel}}); } @@ -270,74 +263,75 @@ struct flowQC { float centrality = getCentrality(collision); + const bool quadraticResponse = cfgQuadraticResponse; + auto maybeSquare = [quadraticResponse](float value) { + return quadraticResponse ? value * value : value; + }; + // EP method - float QmodFT0A_EP = collision.qFT0A(); + float QmodFT0A_EP = maybeSquare(collision.qFT0A()); float psiFT0A_EP = collision.psiFT0A(); - float QxFT0A_EP = QmodFT0A_EP * std::cos(cfgHarmonic.value * psiFT0A_EP); - float QyFT0A_EP = QmodFT0A_EP * std::sin(cfgHarmonic.value * psiFT0A_EP); - float QmodFT0C_EP = collision.qFT0C(); + float QmodFT0C_EP = maybeSquare(collision.qFT0C()); float psiFT0C_EP = collision.psiFT0C(); - float QxFT0C_EP = QmodFT0C_EP * std::cos(cfgHarmonic.value * psiFT0C_EP); - float QyFT0C_EP = QmodFT0C_EP * std::sin(cfgHarmonic.value * psiFT0C_EP); - float QmodTPCl_EP = collision.qTPCL(); + float QmodTPCl_EP = maybeSquare(collision.qTPCL()); float psiTPCl_EP = collision.psiTPCL(); - float QxTPCl_EP = QmodTPCl_EP * std::cos(cfgHarmonic.value * psiTPCl_EP); - float QyTPCl_EP = QmodTPCl_EP * std::sin(cfgHarmonic.value * psiTPCl_EP); - float QmodTPCr_EP = collision.qTPCR(); + float QmodTPCr_EP = maybeSquare(collision.qTPCR()); float psiTPCr_EP = collision.psiTPCR(); - float QxTPCr_EP = QmodTPCr_EP * std::cos(cfgHarmonic.value * psiTPCr_EP); - float QyTPCr_EP = QmodTPCr_EP * std::sin(cfgHarmonic.value * psiTPCr_EP); - float QmodTPC_EP = collision.qTPC(); + float QmodTPC_EP = maybeSquare(collision.qTPC()); float psiTPC_EP = collision.psiTPC(); - float QxTPC_EP = QmodTPC_EP * std::cos(cfgHarmonic.value * psiTPC_EP); - float QyTPC_EP = QmodTPC_EP * std::sin(cfgHarmonic.value * psiTPC_EP); // Qvec method - float QxFT0A_Qvec = collision.qvecFT0AReVec()[cfgHarmonic.value - 2]; - float QyFT0A_Qvec = collision.qvecFT0AImVec()[cfgHarmonic.value - 2]; + float QxFT0A_Qvec_raw = collision.qvecFT0AReVec()[cfgHarmonic.value - 2]; + float QyFT0A_Qvec_raw = collision.qvecFT0AImVec()[cfgHarmonic.value - 2]; + float psiFT0A_Qvec = computeEventPlane(QyFT0A_Qvec_raw, QxFT0A_Qvec_raw); + float QxFT0A_Qvec = maybeSquare(QxFT0A_Qvec_raw); + float QyFT0A_Qvec = maybeSquare(QyFT0A_Qvec_raw); float QmodFT0A_Qvec = std::hypot(QxFT0A_Qvec, QyFT0A_Qvec); - float psiFT0A_Qvec = computeEventPlane(QyFT0A_Qvec, QxFT0A_Qvec); - float QxFT0C_Qvec = collision.qvecFT0CReVec()[cfgHarmonic.value - 2]; - float QyFT0C_Qvec = collision.qvecFT0CImVec()[cfgHarmonic.value - 2]; + float QxFT0C_Qvec_raw = collision.qvecFT0CReVec()[cfgHarmonic.value - 2]; + float QyFT0C_Qvec_raw = collision.qvecFT0CImVec()[cfgHarmonic.value - 2]; + float psiFT0C_Qvec = computeEventPlane(QyFT0C_Qvec_raw, QxFT0C_Qvec_raw); + float QxFT0C_Qvec = maybeSquare(QxFT0C_Qvec_raw); + float QyFT0C_Qvec = maybeSquare(QyFT0C_Qvec_raw); float QmodFT0C_Qvec = std::hypot(QxFT0C_Qvec, QyFT0C_Qvec); - float psiFT0C_Qvec = computeEventPlane(QyFT0C_Qvec, QxFT0C_Qvec); - float QxTPCl_Qvec = collision.qvecTPCnegReVec()[cfgHarmonic.value - 2]; - float QyTPCl_Qvec = collision.qvecTPCnegImVec()[cfgHarmonic.value - 2]; + float QxTPCl_Qvec_raw = collision.qvecTPCnegReVec()[cfgHarmonic.value - 2]; + float QyTPCl_Qvec_raw = collision.qvecTPCnegImVec()[cfgHarmonic.value - 2]; + float psiTPCl_Qvec = computeEventPlane(QyTPCl_Qvec_raw, QxTPCl_Qvec_raw); + float QxTPCl_Qvec = maybeSquare(QxTPCl_Qvec_raw); + float QyTPCl_Qvec = maybeSquare(QyTPCl_Qvec_raw); float QmodTPCl_Qvec = std::hypot(QxTPCl_Qvec, QyTPCl_Qvec); - float psiTPCl_Qvec = computeEventPlane(QyTPCl_Qvec, QxTPCl_Qvec); - float QxTPCr_Qvec = collision.qvecTPCposReVec()[cfgHarmonic.value - 2]; - float QyTPCr_Qvec = collision.qvecTPCposImVec()[cfgHarmonic.value - 2]; + float QxTPCr_Qvec_raw = collision.qvecTPCposReVec()[cfgHarmonic.value - 2]; + float QyTPCr_Qvec_raw = collision.qvecTPCposImVec()[cfgHarmonic.value - 2]; + float psiTPCr_Qvec = computeEventPlane(QyTPCr_Qvec_raw, QxTPCr_Qvec_raw); + float QxTPCr_Qvec = maybeSquare(QxTPCr_Qvec_raw); + float QyTPCr_Qvec = maybeSquare(QyTPCr_Qvec_raw); float QmodTPCr_Qvec = std::hypot(QxTPCr_Qvec, QyTPCr_Qvec); - float psiTPCr_Qvec = computeEventPlane(QyTPCr_Qvec, QxTPCr_Qvec); - float QxTPC_Qvec = collision.qvecTPCallReVec()[cfgHarmonic.value - 2]; - float QyTPC_Qvec = collision.qvecTPCallImVec()[cfgHarmonic.value - 2]; + float QxTPC_Qvec_raw = collision.qvecTPCallReVec()[cfgHarmonic.value - 2]; + float QyTPC_Qvec_raw = collision.qvecTPCallImVec()[cfgHarmonic.value - 2]; + float psiTPC_Qvec = computeEventPlane(QyTPC_Qvec_raw, QxTPC_Qvec_raw); + float QxTPC_Qvec = maybeSquare(QxTPC_Qvec_raw); + float QyTPC_Qvec = maybeSquare(QyTPC_Qvec_raw); float QmodTPC_Qvec = std::hypot(QxTPC_Qvec, QyTPC_Qvec); - float psiTPC_Qvec = computeEventPlane(QyTPC_Qvec, QxTPC_Qvec); - std::array vec_Qx[2] = {{QxFT0C_EP, QxFT0A_EP, QxTPCl_EP, QxTPCr_EP, QxTPC_EP}, {QxFT0C_Qvec, QxFT0A_Qvec, QxTPCl_Qvec, QxTPCr_Qvec, QxTPC_Qvec}}; - std::array vec_Qy[2] = {{QyFT0C_EP, QyFT0A_EP, QyTPCl_EP, QyTPCr_EP, QyTPC_EP}, {QyFT0C_Qvec, QyFT0A_Qvec, QyTPCl_Qvec, QyTPCr_Qvec, QyTPC_Qvec}}; std::array vec_Qmod[2] = {{QmodFT0C_EP, QmodFT0A_EP, QmodTPCl_EP, QmodTPCr_EP, QmodTPC_EP}, {QmodFT0C_Qvec, QmodFT0A_Qvec, QmodTPCl_Qvec, QmodTPCr_Qvec, QmodTPC_Qvec}}; std::array vec_Qpsi[2] = {{psiFT0C_EP, psiFT0A_EP, psiTPCl_EP, psiTPCr_EP, psiTPC_EP}, {psiFT0C_Qvec, psiFT0A_Qvec, psiTPCl_Qvec, psiTPCr_Qvec, psiTPC_Qvec}}; for (int iMethod = 0; iMethod < methods::kNmethods; iMethod++) { for (int iQvecDet = 0; iQvecDet < qVecDetectors::kNqVecDetectors; iQvecDet++) { - hQxQy[iMethod][iQvecDet]->Fill(centrality, vec_Qx[iMethod][iQvecDet], vec_Qy[iMethod][iQvecDet]); - hNormQxQy[iMethod][iQvecDet]->Fill(centrality, vec_Qx[iMethod][iQvecDet] / vec_Qmod[iMethod][iQvecDet], vec_Qy[iMethod][iQvecDet] / vec_Qmod[iMethod][iQvecDet]); hPsi[iMethod][iQvecDet]->Fill(centrality, vec_Qpsi[iMethod][iQvecDet]); for (int jQvecDet = iQvecDet + 1; jQvecDet < qVecDetectors::kNqVecDetectors; jQvecDet++) { // Q-vector azimuthal-angle differences hDeltaPsi[iMethod][iQvecDet][jQvecDet]->Fill(centrality, vec_Qpsi[iMethod][iQvecDet] - vec_Qpsi[iMethod][jQvecDet]); // Scalar-product histograms auto getSP = [&](int iDet1, int iDet2) { - return vec_Qx[iMethod][iDet1] * vec_Qx[iMethod][iDet2] + vec_Qy[iMethod][iDet1] * vec_Qy[iMethod][iDet2]; + return vec_Qmod[iMethod][iDet1] * vec_Qmod[iMethod][iDet2] * std::cos(cfgHarmonic.value * (vec_Qpsi[iMethod][iDet1] - vec_Qpsi[iMethod][iDet2])); }; hScalarProduct[iMethod][iQvecDet][jQvecDet]->Fill(centrality, getSP(iQvecDet, jQvecDet)); // Normalised scalar-product histograms From 48bc30fc430738376abb3572636bdea3e8224db8 Mon Sep 17 00:00:00 2001 From: Luca Barioglio Date: Wed, 20 May 2026 16:13:43 +0200 Subject: [PATCH 2/3] Fix harmonic --- PWGLF/TableProducer/QC/flowQC.cxx | 36 ++++++++++++++++--------------- 1 file changed, 19 insertions(+), 17 deletions(-) diff --git a/PWGLF/TableProducer/QC/flowQC.cxx b/PWGLF/TableProducer/QC/flowQC.cxx index 5b72d66d136..eff0c1db5c2 100644 --- a/PWGLF/TableProducer/QC/flowQC.cxx +++ b/PWGLF/TableProducer/QC/flowQC.cxx @@ -137,9 +137,9 @@ struct flowQC { return collision.sel8() && collision.posZ() > -cfgCutVertex && collision.posZ() < cfgCutVertex && collision.selection_bit(aod::evsel::kNoTimeFrameBorder) && collision.triggereventep() && collision.selection_bit(aod::evsel::kNoSameBunchPileup); } - float computeEventPlane(float y, float x) + float computeEventPlane(float y, float x, float harmonic) { - return 0.5 * TMath::ATan2(y, x); + return (1.f / harmonic) * TMath::ATan2(y, x); } void initCCDB(aod::BCsWithTimestamps::iterator const& bc) @@ -267,6 +267,8 @@ struct flowQC { auto maybeSquare = [quadraticResponse](float value) { return quadraticResponse ? value * value : value; }; + const float qvecHarmonic = quadraticResponse ? (cfgHarmonic.value / 2.f) : cfgHarmonic.value; + const int qvecHarmonicIndex = static_cast(qvecHarmonic) - 2; // EP method float QmodFT0A_EP = maybeSquare(collision.qFT0A()); @@ -285,37 +287,37 @@ struct flowQC { float psiTPC_EP = collision.psiTPC(); // Qvec method - float QxFT0A_Qvec_raw = collision.qvecFT0AReVec()[cfgHarmonic.value - 2]; - float QyFT0A_Qvec_raw = collision.qvecFT0AImVec()[cfgHarmonic.value - 2]; - float psiFT0A_Qvec = computeEventPlane(QyFT0A_Qvec_raw, QxFT0A_Qvec_raw); + float QxFT0A_Qvec_raw = collision.qvecFT0AReVec()[qvecHarmonicIndex]; + float QyFT0A_Qvec_raw = collision.qvecFT0AImVec()[qvecHarmonicIndex]; + float psiFT0A_Qvec = computeEventPlane(QyFT0A_Qvec_raw, QxFT0A_Qvec_raw, qvecHarmonic); float QxFT0A_Qvec = maybeSquare(QxFT0A_Qvec_raw); float QyFT0A_Qvec = maybeSquare(QyFT0A_Qvec_raw); float QmodFT0A_Qvec = std::hypot(QxFT0A_Qvec, QyFT0A_Qvec); - float QxFT0C_Qvec_raw = collision.qvecFT0CReVec()[cfgHarmonic.value - 2]; - float QyFT0C_Qvec_raw = collision.qvecFT0CImVec()[cfgHarmonic.value - 2]; - float psiFT0C_Qvec = computeEventPlane(QyFT0C_Qvec_raw, QxFT0C_Qvec_raw); + float QxFT0C_Qvec_raw = collision.qvecFT0CReVec()[qvecHarmonicIndex]; + float QyFT0C_Qvec_raw = collision.qvecFT0CImVec()[qvecHarmonicIndex]; + float psiFT0C_Qvec = computeEventPlane(QyFT0C_Qvec_raw, QxFT0C_Qvec_raw, qvecHarmonic); float QxFT0C_Qvec = maybeSquare(QxFT0C_Qvec_raw); float QyFT0C_Qvec = maybeSquare(QyFT0C_Qvec_raw); float QmodFT0C_Qvec = std::hypot(QxFT0C_Qvec, QyFT0C_Qvec); - float QxTPCl_Qvec_raw = collision.qvecTPCnegReVec()[cfgHarmonic.value - 2]; - float QyTPCl_Qvec_raw = collision.qvecTPCnegImVec()[cfgHarmonic.value - 2]; - float psiTPCl_Qvec = computeEventPlane(QyTPCl_Qvec_raw, QxTPCl_Qvec_raw); + float QxTPCl_Qvec_raw = collision.qvecTPCnegReVec()[qvecHarmonicIndex]; + float QyTPCl_Qvec_raw = collision.qvecTPCnegImVec()[qvecHarmonicIndex]; + float psiTPCl_Qvec = computeEventPlane(QyTPCl_Qvec_raw, QxTPCl_Qvec_raw, qvecHarmonic); float QxTPCl_Qvec = maybeSquare(QxTPCl_Qvec_raw); float QyTPCl_Qvec = maybeSquare(QyTPCl_Qvec_raw); float QmodTPCl_Qvec = std::hypot(QxTPCl_Qvec, QyTPCl_Qvec); - float QxTPCr_Qvec_raw = collision.qvecTPCposReVec()[cfgHarmonic.value - 2]; - float QyTPCr_Qvec_raw = collision.qvecTPCposImVec()[cfgHarmonic.value - 2]; - float psiTPCr_Qvec = computeEventPlane(QyTPCr_Qvec_raw, QxTPCr_Qvec_raw); + float QxTPCr_Qvec_raw = collision.qvecTPCposReVec()[qvecHarmonicIndex]; + float QyTPCr_Qvec_raw = collision.qvecTPCposImVec()[qvecHarmonicIndex]; + float psiTPCr_Qvec = computeEventPlane(QyTPCr_Qvec_raw, QxTPCr_Qvec_raw, qvecHarmonic); float QxTPCr_Qvec = maybeSquare(QxTPCr_Qvec_raw); float QyTPCr_Qvec = maybeSquare(QyTPCr_Qvec_raw); float QmodTPCr_Qvec = std::hypot(QxTPCr_Qvec, QyTPCr_Qvec); - float QxTPC_Qvec_raw = collision.qvecTPCallReVec()[cfgHarmonic.value - 2]; - float QyTPC_Qvec_raw = collision.qvecTPCallImVec()[cfgHarmonic.value - 2]; - float psiTPC_Qvec = computeEventPlane(QyTPC_Qvec_raw, QxTPC_Qvec_raw); + float QxTPC_Qvec_raw = collision.qvecTPCallReVec()[qvecHarmonicIndex]; + float QyTPC_Qvec_raw = collision.qvecTPCallImVec()[qvecHarmonicIndex]; + float psiTPC_Qvec = computeEventPlane(QyTPC_Qvec_raw, QxTPC_Qvec_raw, qvecHarmonic); float QxTPC_Qvec = maybeSquare(QxTPC_Qvec_raw); float QyTPC_Qvec = maybeSquare(QyTPC_Qvec_raw); float QmodTPC_Qvec = std::hypot(QxTPC_Qvec, QyTPC_Qvec); From ece3ae2a2c51a1a2617e88537c0901350bebb860 Mon Sep 17 00:00:00 2001 From: Luca Barioglio Date: Wed, 20 May 2026 16:41:39 +0200 Subject: [PATCH 3/3] Very minor change --- PWGLF/TableProducer/QC/flowQC.cxx | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/PWGLF/TableProducer/QC/flowQC.cxx b/PWGLF/TableProducer/QC/flowQC.cxx index eff0c1db5c2..ae946265d88 100644 --- a/PWGLF/TableProducer/QC/flowQC.cxx +++ b/PWGLF/TableProducer/QC/flowQC.cxx @@ -263,11 +263,10 @@ struct flowQC { float centrality = getCentrality(collision); - const bool quadraticResponse = cfgQuadraticResponse; - auto maybeSquare = [quadraticResponse](float value) { - return quadraticResponse ? value * value : value; + auto maybeSquare = [this](float value) { + return cfgQuadraticResponse ? value * value : value; }; - const float qvecHarmonic = quadraticResponse ? (cfgHarmonic.value / 2.f) : cfgHarmonic.value; + const float qvecHarmonic = cfgQuadraticResponse ? (cfgHarmonic.value / 2.f) : cfgHarmonic.value; const int qvecHarmonicIndex = static_cast(qvecHarmonic) - 2; // EP method