diff --git a/PWGLF/TableProducer/QC/flowQC.cxx b/PWGLF/TableProducer/QC/flowQC.cxx index 2d590e0df69..ae946265d88 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::Join -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) @@ -183,11 +182,7 @@ struct flowQC { const AxisSpec centAxis{cfgCentralityBins, fmt::format("{} percentile", (std::string)centDetectorNames[cfgCentralityEstimator])}; - const AxisSpec QxAxis{cfgQvecBins, Form("Q_{%.0f,x}", cfgHarmonic.value)}; - const AxisSpec QyAxis{cfgQvecBins, Form("Q_{%.0f,y}", cfgHarmonic.value)}; - - const AxisSpec NormQxAxis{cfgQvecBins, Form("#frac{Q_{%.0f,x}}{||#vec{Q_{%.0f}}||}", cfgHarmonic.value, cfgHarmonic.value)}; - const AxisSpec NormQyAxis{cfgQvecBins, Form("#frac{Q_{%.0f,y}}{||#vec{Q_{%.0f}}||}", cfgHarmonic.value, cfgHarmonic.value)}; + const char* qLabel = cfgQuadraticResponse ? "Q^{2}" : "Q"; const AxisSpec psiAxis{cfgPhiBins, Form("#psi_{%.0f}", cfgHarmonic.value)}; const AxisSpec psiCompAxis{cfgPhiBins, Form("#psi_{%.0f}^{EP} - #psi_{%.0f}^{Qvec}", cfgHarmonic.value, cfgHarmonic.value)}; @@ -206,8 +201,6 @@ struct flowQC { HistogramRegistry* registry = (iMethod == methods::kEP) ? &flow_ep : &flow_qvec; for (int iQvecDet = 0; iQvecDet < qVecDetectors::kNqVecDetectors; iQvecDet++) { - hQxQy[iMethod][iQvecDet] = registry->add(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,76 @@ struct flowQC { float centrality = getCentrality(collision); + auto maybeSquare = [this](float value) { + return cfgQuadraticResponse ? value * value : value; + }; + const float qvecHarmonic = cfgQuadraticResponse ? (cfgHarmonic.value / 2.f) : cfgHarmonic.value; + const int qvecHarmonicIndex = static_cast(qvecHarmonic) - 2; + // 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()[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 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()[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 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()[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 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()[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 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()[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); - 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