Skip to content

Commit 4bf9224

Browse files
authored
[PWGLF] Add quadratic response option (#16323)
1 parent 57ab50c commit 4bf9224

1 file changed

Lines changed: 43 additions & 48 deletions

File tree

PWGLF/TableProducer/QC/flowQC.cxx

Lines changed: 43 additions & 48 deletions
Original file line numberDiff line numberDiff line change
@@ -89,8 +89,6 @@ enum methods {
8989
};
9090
static const std::vector<std::string> suffixes = {"EP", "Qvec"};
9191

92-
std::shared_ptr<TH3> hQxQy[kNmethods][kNqVecDetectors];
93-
std::shared_ptr<TH3> hNormQxQy[kNmethods][kNqVecDetectors];
9492
std::shared_ptr<TH2> hPsi[kNmethods][kNqVecDetectors];
9593
std::shared_ptr<TH2> hDeltaPsi[kNmethods][kNqVecDetectors][kNqVecDetectors];
9694
std::shared_ptr<TH2> hScalarProduct[kNmethods][kNqVecDetectors][kNqVecDetectors];
@@ -122,6 +120,7 @@ struct flowQC {
122120
float mBz = 0.f;
123121

124122
Configurable<float> cfgHarmonic{"cfgHarmonic", 2.f, "Harmonics for flow analysis"};
123+
Configurable<bool> cfgQuadraticResponse{"cfgQuadraticResponse", false, "Use quadratic response for Q-vector quantities"};
125124

126125
// Flow analysis
127126
using CollWithEPandQvec = soa::Join<aod::Collisions,
@@ -138,9 +137,9 @@ struct flowQC {
138137
return collision.sel8() && collision.posZ() > -cfgCutVertex && collision.posZ() < cfgCutVertex && collision.selection_bit(aod::evsel::kNoTimeFrameBorder) && collision.triggereventep() && collision.selection_bit(aod::evsel::kNoSameBunchPileup);
139138
}
140139

141-
float computeEventPlane(float y, float x)
140+
float computeEventPlane(float y, float x, float harmonic)
142141
{
143-
return 0.5 * TMath::ATan2(y, x);
142+
return (1.f / harmonic) * TMath::ATan2(y, x);
144143
}
145144

146145
void initCCDB(aod::BCsWithTimestamps::iterator const& bc)
@@ -183,11 +182,7 @@ struct flowQC {
183182

184183
const AxisSpec centAxis{cfgCentralityBins, fmt::format("{} percentile", (std::string)centDetectorNames[cfgCentralityEstimator])};
185184

186-
const AxisSpec QxAxis{cfgQvecBins, Form("Q_{%.0f,x}", cfgHarmonic.value)};
187-
const AxisSpec QyAxis{cfgQvecBins, Form("Q_{%.0f,y}", cfgHarmonic.value)};
188-
189-
const AxisSpec NormQxAxis{cfgQvecBins, Form("#frac{Q_{%.0f,x}}{||#vec{Q_{%.0f}}||}", cfgHarmonic.value, cfgHarmonic.value)};
190-
const AxisSpec NormQyAxis{cfgQvecBins, Form("#frac{Q_{%.0f,y}}{||#vec{Q_{%.0f}}||}", cfgHarmonic.value, cfgHarmonic.value)};
185+
const char* qLabel = cfgQuadraticResponse ? "Q^{2}" : "Q";
191186

192187
const AxisSpec psiAxis{cfgPhiBins, Form("#psi_{%.0f}", cfgHarmonic.value)};
193188
const AxisSpec psiCompAxis{cfgPhiBins, Form("#psi_{%.0f}^{EP} - #psi_{%.0f}^{Qvec}", cfgHarmonic.value, cfgHarmonic.value)};
@@ -206,21 +201,19 @@ struct flowQC {
206201
HistogramRegistry* registry = (iMethod == methods::kEP) ? &flow_ep : &flow_qvec;
207202

208203
for (int iQvecDet = 0; iQvecDet < qVecDetectors::kNqVecDetectors; iQvecDet++) {
209-
hQxQy[iMethod][iQvecDet] = registry->add<TH3>(Form("hQxQy_%s_%s", qVecDetectorNames[iQvecDet].c_str(), suffixes[iMethod].c_str()), "", HistType::kTH3F, {centAxis, QxAxis, QyAxis});
210-
hNormQxQy[iMethod][iQvecDet] = registry->add<TH3>(Form("hNormQxQy_%s_%s", qVecDetectorNames[iQvecDet].c_str(), suffixes[iMethod].c_str()), "", HistType::kTH3F, {centAxis, NormQxAxis, NormQyAxis});
211204
hPsi[iMethod][iQvecDet] = registry->add<TH2>(Form("hPsi_%s_%s", qVecDetectorNames[iQvecDet].c_str(), suffixes[iMethod].c_str()), "", HistType::kTH2F, {centAxis, psiAxis});
212205
for (int jQvecDet = iQvecDet + 1; jQvecDet < qVecDetectors::kNqVecDetectors; jQvecDet++) {
213206

214207
// Q-vector azimuthal-angle differences
215208
hDeltaPsi[iMethod][iQvecDet][jQvecDet] = registry->add<TH2>(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())}});
216209

217210
// Scalar-product histograms
218-
auto spLabel = Form("#vec{Q}_{%.0f}^{%s} #upoint #vec{Q}_{%.0f}^{%s}", cfgHarmonic.value, qVecDetectorNames[iQvecDet].c_str(), cfgHarmonic.value, qVecDetectorNames[jQvecDet].c_str());
211+
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());
219212

220213
hScalarProduct[iMethod][iQvecDet][jQvecDet] = registry->add<TH2>(Form("hScalarProduct_%s_%s_%s", qVecDetectorNames[iQvecDet].c_str(), qVecDetectorNames[jQvecDet].c_str(), suffixes[iMethod].c_str()), "", HistType::kTH2F, {centAxis, {cfgQvecBins, spLabel}});
221214

222215
// Normalised scalar-product histograms
223-
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());
216+
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());
224217

225218
hNormalisedScalarProduct[iMethod][iQvecDet][jQvecDet] = registry->add<TH2>(Form("hNormalisedScalarProduct_%s_%s_%s", qVecDetectorNames[iQvecDet].c_str(), qVecDetectorNames[jQvecDet].c_str(), suffixes[iMethod].c_str()), "", HistType::kTH2F, {centAxis, {cfgQvecBins, normSpLabel}});
226219
}
@@ -270,74 +263,76 @@ struct flowQC {
270263

271264
float centrality = getCentrality(collision);
272265

266+
auto maybeSquare = [this](float value) {
267+
return cfgQuadraticResponse ? value * value : value;
268+
};
269+
const float qvecHarmonic = cfgQuadraticResponse ? (cfgHarmonic.value / 2.f) : cfgHarmonic.value;
270+
const int qvecHarmonicIndex = static_cast<int>(qvecHarmonic) - 2;
271+
273272
// EP method
274-
float QmodFT0A_EP = collision.qFT0A();
273+
float QmodFT0A_EP = maybeSquare(collision.qFT0A());
275274
float psiFT0A_EP = collision.psiFT0A();
276-
float QxFT0A_EP = QmodFT0A_EP * std::cos(cfgHarmonic.value * psiFT0A_EP);
277-
float QyFT0A_EP = QmodFT0A_EP * std::sin(cfgHarmonic.value * psiFT0A_EP);
278275

279-
float QmodFT0C_EP = collision.qFT0C();
276+
float QmodFT0C_EP = maybeSquare(collision.qFT0C());
280277
float psiFT0C_EP = collision.psiFT0C();
281-
float QxFT0C_EP = QmodFT0C_EP * std::cos(cfgHarmonic.value * psiFT0C_EP);
282-
float QyFT0C_EP = QmodFT0C_EP * std::sin(cfgHarmonic.value * psiFT0C_EP);
283278

284-
float QmodTPCl_EP = collision.qTPCL();
279+
float QmodTPCl_EP = maybeSquare(collision.qTPCL());
285280
float psiTPCl_EP = collision.psiTPCL();
286-
float QxTPCl_EP = QmodTPCl_EP * std::cos(cfgHarmonic.value * psiTPCl_EP);
287-
float QyTPCl_EP = QmodTPCl_EP * std::sin(cfgHarmonic.value * psiTPCl_EP);
288281

289-
float QmodTPCr_EP = collision.qTPCR();
282+
float QmodTPCr_EP = maybeSquare(collision.qTPCR());
290283
float psiTPCr_EP = collision.psiTPCR();
291-
float QxTPCr_EP = QmodTPCr_EP * std::cos(cfgHarmonic.value * psiTPCr_EP);
292-
float QyTPCr_EP = QmodTPCr_EP * std::sin(cfgHarmonic.value * psiTPCr_EP);
293284

294-
float QmodTPC_EP = collision.qTPC();
285+
float QmodTPC_EP = maybeSquare(collision.qTPC());
295286
float psiTPC_EP = collision.psiTPC();
296-
float QxTPC_EP = QmodTPC_EP * std::cos(cfgHarmonic.value * psiTPC_EP);
297-
float QyTPC_EP = QmodTPC_EP * std::sin(cfgHarmonic.value * psiTPC_EP);
298287

299288
// Qvec method
300-
float QxFT0A_Qvec = collision.qvecFT0AReVec()[cfgHarmonic.value - 2];
301-
float QyFT0A_Qvec = collision.qvecFT0AImVec()[cfgHarmonic.value - 2];
289+
float QxFT0A_Qvec_raw = collision.qvecFT0AReVec()[qvecHarmonicIndex];
290+
float QyFT0A_Qvec_raw = collision.qvecFT0AImVec()[qvecHarmonicIndex];
291+
float psiFT0A_Qvec = computeEventPlane(QyFT0A_Qvec_raw, QxFT0A_Qvec_raw, qvecHarmonic);
292+
float QxFT0A_Qvec = maybeSquare(QxFT0A_Qvec_raw);
293+
float QyFT0A_Qvec = maybeSquare(QyFT0A_Qvec_raw);
302294
float QmodFT0A_Qvec = std::hypot(QxFT0A_Qvec, QyFT0A_Qvec);
303-
float psiFT0A_Qvec = computeEventPlane(QyFT0A_Qvec, QxFT0A_Qvec);
304295

305-
float QxFT0C_Qvec = collision.qvecFT0CReVec()[cfgHarmonic.value - 2];
306-
float QyFT0C_Qvec = collision.qvecFT0CImVec()[cfgHarmonic.value - 2];
296+
float QxFT0C_Qvec_raw = collision.qvecFT0CReVec()[qvecHarmonicIndex];
297+
float QyFT0C_Qvec_raw = collision.qvecFT0CImVec()[qvecHarmonicIndex];
298+
float psiFT0C_Qvec = computeEventPlane(QyFT0C_Qvec_raw, QxFT0C_Qvec_raw, qvecHarmonic);
299+
float QxFT0C_Qvec = maybeSquare(QxFT0C_Qvec_raw);
300+
float QyFT0C_Qvec = maybeSquare(QyFT0C_Qvec_raw);
307301
float QmodFT0C_Qvec = std::hypot(QxFT0C_Qvec, QyFT0C_Qvec);
308-
float psiFT0C_Qvec = computeEventPlane(QyFT0C_Qvec, QxFT0C_Qvec);
309302

310-
float QxTPCl_Qvec = collision.qvecTPCnegReVec()[cfgHarmonic.value - 2];
311-
float QyTPCl_Qvec = collision.qvecTPCnegImVec()[cfgHarmonic.value - 2];
303+
float QxTPCl_Qvec_raw = collision.qvecTPCnegReVec()[qvecHarmonicIndex];
304+
float QyTPCl_Qvec_raw = collision.qvecTPCnegImVec()[qvecHarmonicIndex];
305+
float psiTPCl_Qvec = computeEventPlane(QyTPCl_Qvec_raw, QxTPCl_Qvec_raw, qvecHarmonic);
306+
float QxTPCl_Qvec = maybeSquare(QxTPCl_Qvec_raw);
307+
float QyTPCl_Qvec = maybeSquare(QyTPCl_Qvec_raw);
312308
float QmodTPCl_Qvec = std::hypot(QxTPCl_Qvec, QyTPCl_Qvec);
313-
float psiTPCl_Qvec = computeEventPlane(QyTPCl_Qvec, QxTPCl_Qvec);
314309

315-
float QxTPCr_Qvec = collision.qvecTPCposReVec()[cfgHarmonic.value - 2];
316-
float QyTPCr_Qvec = collision.qvecTPCposImVec()[cfgHarmonic.value - 2];
310+
float QxTPCr_Qvec_raw = collision.qvecTPCposReVec()[qvecHarmonicIndex];
311+
float QyTPCr_Qvec_raw = collision.qvecTPCposImVec()[qvecHarmonicIndex];
312+
float psiTPCr_Qvec = computeEventPlane(QyTPCr_Qvec_raw, QxTPCr_Qvec_raw, qvecHarmonic);
313+
float QxTPCr_Qvec = maybeSquare(QxTPCr_Qvec_raw);
314+
float QyTPCr_Qvec = maybeSquare(QyTPCr_Qvec_raw);
317315
float QmodTPCr_Qvec = std::hypot(QxTPCr_Qvec, QyTPCr_Qvec);
318-
float psiTPCr_Qvec = computeEventPlane(QyTPCr_Qvec, QxTPCr_Qvec);
319316

320-
float QxTPC_Qvec = collision.qvecTPCallReVec()[cfgHarmonic.value - 2];
321-
float QyTPC_Qvec = collision.qvecTPCallImVec()[cfgHarmonic.value - 2];
317+
float QxTPC_Qvec_raw = collision.qvecTPCallReVec()[qvecHarmonicIndex];
318+
float QyTPC_Qvec_raw = collision.qvecTPCallImVec()[qvecHarmonicIndex];
319+
float psiTPC_Qvec = computeEventPlane(QyTPC_Qvec_raw, QxTPC_Qvec_raw, qvecHarmonic);
320+
float QxTPC_Qvec = maybeSquare(QxTPC_Qvec_raw);
321+
float QyTPC_Qvec = maybeSquare(QyTPC_Qvec_raw);
322322
float QmodTPC_Qvec = std::hypot(QxTPC_Qvec, QyTPC_Qvec);
323-
float psiTPC_Qvec = computeEventPlane(QyTPC_Qvec, QxTPC_Qvec);
324323

325-
std::array<float, qVecDetectors::kNqVecDetectors> vec_Qx[2] = {{QxFT0C_EP, QxFT0A_EP, QxTPCl_EP, QxTPCr_EP, QxTPC_EP}, {QxFT0C_Qvec, QxFT0A_Qvec, QxTPCl_Qvec, QxTPCr_Qvec, QxTPC_Qvec}};
326-
std::array<float, qVecDetectors::kNqVecDetectors> vec_Qy[2] = {{QyFT0C_EP, QyFT0A_EP, QyTPCl_EP, QyTPCr_EP, QyTPC_EP}, {QyFT0C_Qvec, QyFT0A_Qvec, QyTPCl_Qvec, QyTPCr_Qvec, QyTPC_Qvec}};
327324
std::array<float, qVecDetectors::kNqVecDetectors> vec_Qmod[2] = {{QmodFT0C_EP, QmodFT0A_EP, QmodTPCl_EP, QmodTPCr_EP, QmodTPC_EP}, {QmodFT0C_Qvec, QmodFT0A_Qvec, QmodTPCl_Qvec, QmodTPCr_Qvec, QmodTPC_Qvec}};
328325
std::array<float, qVecDetectors::kNqVecDetectors> vec_Qpsi[2] = {{psiFT0C_EP, psiFT0A_EP, psiTPCl_EP, psiTPCr_EP, psiTPC_EP}, {psiFT0C_Qvec, psiFT0A_Qvec, psiTPCl_Qvec, psiTPCr_Qvec, psiTPC_Qvec}};
329326

330327
for (int iMethod = 0; iMethod < methods::kNmethods; iMethod++) {
331328
for (int iQvecDet = 0; iQvecDet < qVecDetectors::kNqVecDetectors; iQvecDet++) {
332-
hQxQy[iMethod][iQvecDet]->Fill(centrality, vec_Qx[iMethod][iQvecDet], vec_Qy[iMethod][iQvecDet]);
333-
hNormQxQy[iMethod][iQvecDet]->Fill(centrality, vec_Qx[iMethod][iQvecDet] / vec_Qmod[iMethod][iQvecDet], vec_Qy[iMethod][iQvecDet] / vec_Qmod[iMethod][iQvecDet]);
334329
hPsi[iMethod][iQvecDet]->Fill(centrality, vec_Qpsi[iMethod][iQvecDet]);
335330
for (int jQvecDet = iQvecDet + 1; jQvecDet < qVecDetectors::kNqVecDetectors; jQvecDet++) {
336331
// Q-vector azimuthal-angle differences
337332
hDeltaPsi[iMethod][iQvecDet][jQvecDet]->Fill(centrality, vec_Qpsi[iMethod][iQvecDet] - vec_Qpsi[iMethod][jQvecDet]);
338333
// Scalar-product histograms
339334
auto getSP = [&](int iDet1, int iDet2) {
340-
return vec_Qx[iMethod][iDet1] * vec_Qx[iMethod][iDet2] + vec_Qy[iMethod][iDet1] * vec_Qy[iMethod][iDet2];
335+
return vec_Qmod[iMethod][iDet1] * vec_Qmod[iMethod][iDet2] * std::cos(cfgHarmonic.value * (vec_Qpsi[iMethod][iDet1] - vec_Qpsi[iMethod][iDet2]));
341336
};
342337
hScalarProduct[iMethod][iQvecDet][jQvecDet]->Fill(centrality, getSP(iQvecDet, jQvecDet));
343338
// Normalised scalar-product histograms

0 commit comments

Comments
 (0)