Skip to content

Commit f817854

Browse files
lbariogllbariogl
committed
Add quadratic response option
1 parent 43dfb1c commit f817854

1 file changed

Lines changed: 29 additions & 21 deletions

File tree

‎PWGLF/TableProducer/QC/flowQC.cxx‎

Lines changed: 29 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -122,6 +122,7 @@ struct flowQC {
122122
float mBz = 0.f;
123123

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

126127
// Flow analysis
127128
using CollWithEPandQvec = soa::Join<aod::Collisions,
@@ -183,11 +184,13 @@ struct flowQC {
183184

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

186-
const AxisSpec QxAxis{cfgQvecBins, Form("Q_{%.0f,x}", cfgHarmonic.value)};
187-
const AxisSpec QyAxis{cfgQvecBins, Form("Q_{%.0f,y}", cfgHarmonic.value)};
187+
const char* qLabel = cfgQuadraticResponse ? "Q^{2}" : "Q";
188188

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)};
189+
const AxisSpec QxAxis{cfgQvecBins, Form("%s_{%.0f,x}", qLabel, cfgHarmonic.value)};
190+
const AxisSpec QyAxis{cfgQvecBins, Form("%s_{%.0f,y}", qLabel, cfgHarmonic.value)};
191+
192+
const AxisSpec NormQxAxis{cfgQvecBins, Form("#frac{%s_{%.0f,x}}{||#vec{%s}_{%.0f}||}", qLabel, cfgHarmonic.value, qLabel, cfgHarmonic.value)};
193+
const AxisSpec NormQyAxis{cfgQvecBins, Form("#frac{%s_{%.0f,y}}{||#vec{%s}_{%.0f}||}", qLabel, cfgHarmonic.value, qLabel, cfgHarmonic.value)};
191194

192195
const AxisSpec psiAxis{cfgPhiBins, Form("#psi_{%.0f}", cfgHarmonic.value)};
193196
const AxisSpec psiCompAxis{cfgPhiBins, Form("#psi_{%.0f}^{EP} - #psi_{%.0f}^{Qvec}", cfgHarmonic.value, cfgHarmonic.value)};
@@ -215,12 +218,12 @@ struct flowQC {
215218
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())}});
216219

217220
// 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());
221+
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());
219222

220223
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}});
221224

222225
// 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());
226+
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());
224227

225228
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}});
226229
}
@@ -270,55 +273,60 @@ struct flowQC {
270273

271274
float centrality = getCentrality(collision);
272275

276+
const bool quadraticResponse = cfgQuadraticResponse;
277+
auto maybeSquare = [quadraticResponse](float value) {
278+
return quadraticResponse ? value * value : value;
279+
};
280+
273281
// EP method
274-
float QmodFT0A_EP = collision.qFT0A();
282+
float QmodFT0A_EP = maybeSquare(collision.qFT0A());
275283
float psiFT0A_EP = collision.psiFT0A();
276284
float QxFT0A_EP = QmodFT0A_EP * std::cos(cfgHarmonic.value * psiFT0A_EP);
277285
float QyFT0A_EP = QmodFT0A_EP * std::sin(cfgHarmonic.value * psiFT0A_EP);
278286

279-
float QmodFT0C_EP = collision.qFT0C();
287+
float QmodFT0C_EP = maybeSquare(collision.qFT0C());
280288
float psiFT0C_EP = collision.psiFT0C();
281289
float QxFT0C_EP = QmodFT0C_EP * std::cos(cfgHarmonic.value * psiFT0C_EP);
282290
float QyFT0C_EP = QmodFT0C_EP * std::sin(cfgHarmonic.value * psiFT0C_EP);
283291

284-
float QmodTPCl_EP = collision.qTPCL();
292+
float QmodTPCl_EP = maybeSquare(collision.qTPCL());
285293
float psiTPCl_EP = collision.psiTPCL();
286294
float QxTPCl_EP = QmodTPCl_EP * std::cos(cfgHarmonic.value * psiTPCl_EP);
287295
float QyTPCl_EP = QmodTPCl_EP * std::sin(cfgHarmonic.value * psiTPCl_EP);
288296

289-
float QmodTPCr_EP = collision.qTPCR();
297+
float QmodTPCr_EP = maybeSquare(collision.qTPCR());
290298
float psiTPCr_EP = collision.psiTPCR();
291299
float QxTPCr_EP = QmodTPCr_EP * std::cos(cfgHarmonic.value * psiTPCr_EP);
292300
float QyTPCr_EP = QmodTPCr_EP * std::sin(cfgHarmonic.value * psiTPCr_EP);
293301

294-
float QmodTPC_EP = collision.qTPC();
302+
float QmodTPC_EP = maybeSquare(collision.qTPC());
295303
float psiTPC_EP = collision.psiTPC();
296304
float QxTPC_EP = QmodTPC_EP * std::cos(cfgHarmonic.value * psiTPC_EP);
297305
float QyTPC_EP = QmodTPC_EP * std::sin(cfgHarmonic.value * psiTPC_EP);
298306

299307
// Qvec method
300-
float QxFT0A_Qvec = collision.qvecFT0AReVec()[cfgHarmonic.value - 2];
301-
float QyFT0A_Qvec = collision.qvecFT0AImVec()[cfgHarmonic.value - 2];
308+
float QxFT0A_Qvec = maybeSquare(collision.qvecFT0AReVec()[cfgHarmonic.value - 2]);
309+
float QyFT0A_Qvec = maybeSquare(collision.qvecFT0AImVec()[cfgHarmonic.value - 2]);
302310
float QmodFT0A_Qvec = std::hypot(QxFT0A_Qvec, QyFT0A_Qvec);
303311
float psiFT0A_Qvec = computeEventPlane(QyFT0A_Qvec, QxFT0A_Qvec);
304312

305-
float QxFT0C_Qvec = collision.qvecFT0CReVec()[cfgHarmonic.value - 2];
306-
float QyFT0C_Qvec = collision.qvecFT0CImVec()[cfgHarmonic.value - 2];
313+
float QxFT0C_Qvec = maybeSquare(collision.qvecFT0CReVec()[cfgHarmonic.value - 2]);
314+
float QyFT0C_Qvec = maybeSquare(collision.qvecFT0CImVec()[cfgHarmonic.value - 2]);
307315
float QmodFT0C_Qvec = std::hypot(QxFT0C_Qvec, QyFT0C_Qvec);
308316
float psiFT0C_Qvec = computeEventPlane(QyFT0C_Qvec, QxFT0C_Qvec);
309317

310-
float QxTPCl_Qvec = collision.qvecTPCnegReVec()[cfgHarmonic.value - 2];
311-
float QyTPCl_Qvec = collision.qvecTPCnegImVec()[cfgHarmonic.value - 2];
318+
float QxTPCl_Qvec = maybeSquare(collision.qvecTPCnegReVec()[cfgHarmonic.value - 2]);
319+
float QyTPCl_Qvec = maybeSquare(collision.qvecTPCnegImVec()[cfgHarmonic.value - 2]);
312320
float QmodTPCl_Qvec = std::hypot(QxTPCl_Qvec, QyTPCl_Qvec);
313321
float psiTPCl_Qvec = computeEventPlane(QyTPCl_Qvec, QxTPCl_Qvec);
314322

315-
float QxTPCr_Qvec = collision.qvecTPCposReVec()[cfgHarmonic.value - 2];
316-
float QyTPCr_Qvec = collision.qvecTPCposImVec()[cfgHarmonic.value - 2];
323+
float QxTPCr_Qvec = maybeSquare(collision.qvecTPCposReVec()[cfgHarmonic.value - 2]);
324+
float QyTPCr_Qvec = maybeSquare(collision.qvecTPCposImVec()[cfgHarmonic.value - 2]);
317325
float QmodTPCr_Qvec = std::hypot(QxTPCr_Qvec, QyTPCr_Qvec);
318326
float psiTPCr_Qvec = computeEventPlane(QyTPCr_Qvec, QxTPCr_Qvec);
319327

320-
float QxTPC_Qvec = collision.qvecTPCallReVec()[cfgHarmonic.value - 2];
321-
float QyTPC_Qvec = collision.qvecTPCallImVec()[cfgHarmonic.value - 2];
328+
float QxTPC_Qvec = maybeSquare(collision.qvecTPCallReVec()[cfgHarmonic.value - 2]);
329+
float QyTPC_Qvec = maybeSquare(collision.qvecTPCallImVec()[cfgHarmonic.value - 2]);
322330
float QmodTPC_Qvec = std::hypot(QxTPC_Qvec, QyTPC_Qvec);
323331
float psiTPC_Qvec = computeEventPlane(QyTPC_Qvec, QxTPC_Qvec);
324332

0 commit comments

Comments
 (0)