Skip to content

[PWGJE] add psi4, fix psi3 bug using Q2 #16776

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
nzardosh merged 2 commits into from
Jun 24, 2026
+83 −59
Merged

[PWGJE] add psi4, fix psi3 bug using Q2 #16776

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
993db1a
add psi4, fix psi3 bug using Q2
joachimckh Jun 23, 2026
7e72f7a
zz
joachimckh Jun 23, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
142 changes: 83 additions & 59 deletions PWGJE/Tasks/jetSpectraEseTask.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -452,20 +452,32 @@ struct JetSpectraEseTask {
if (doprocessESEEPData) {
LOGF(info, "JetSpectraEseTask::init() - Event Plane Process");
registry.add("eventQA/hPsi2FT0C", ";Centrality; #Psi_{2}", {HistType::kTH2F, {{centAxis}, {150, -2.5, 2.5}}});
registry.addClone("eventQA/hPsi2FT0C", "hPsi2FT0A");
registry.addClone("eventQA/hPsi2FT0C", "hPsi2FV0A");
registry.addClone("eventQA/hPsi2FT0C", "hPsi2TPCpos");
registry.addClone("eventQA/hPsi2FT0C", "hPsi2TPCneg");
registry.addClone("eventQA/hPsi2FT0C", "eventQA/hPsi2FT0A");
registry.addClone("eventQA/hPsi2FT0C", "eventQA/hPsi2FV0A");
registry.addClone("eventQA/hPsi2FT0C", "eventQA/hPsi2TPCpos");
registry.addClone("eventQA/hPsi2FT0C", "eventQA/hPsi2TPCneg");

registry.addClone("eventQA/hPsi2FT0C", "eventQA/hPsi3FT0C");
registry.addClone("eventQA/hPsi2FT0C", "eventQA/hPsi3FT0A");

registry.addClone("eventQA/hPsi2FT0C", "eventQA/hPsi4FT0C");
registry.addClone("eventQA/hPsi2FT0C", "eventQA/hPsi4FT0A");

registry.add("eventQA/hCosPsi2AmC", ";Centrality;cos(2(#Psi_{2}^{A}-#Psi_{2}^{B}));#it{q}_{2}", {HistType::kTH3F, {{centAxis}, {cosAxis}, {eseAxis}}});
registry.addClone("eventQA/hCosPsi2AmC", "hCosPsi2AmB");
registry.addClone("eventQA/hCosPsi2AmC", "hCosPsi2BmC");
registry.addClone("eventQA/hCosPsi2AmC", "eventQA/hCosPsi2AmB");
registry.addClone("eventQA/hCosPsi2AmC", "eventQA/hCosPsi2BmC");

registry.addClone("eventQA/hCosPsi2AmC", "eventQA/hCosPsi4AmC");
registry.addClone("eventQA/hCosPsi4AmC", "eventQA/hCosPsi4AmB");
registry.addClone("eventQA/hCosPsi4AmC", "eventQA/hCosPsi4BmC");

registry.add("eventQA/hQvecUncorV2", ";Centrality;Q_x;Q_y", {HistType::kTH3F, {{centAxis}, {qvecAxis}, {qvecAxis}}});
registry.addClone("eventQA/hQvecUncorV2", "hQvecRectrV2");
registry.addClone("eventQA/hQvecUncorV2", "hQvecTwistV2");
registry.addClone("eventQA/hQvecUncorV2", "hQvecFinalV2");
registry.addClone("eventQA/hPsi2FT0C", "hEPUncorV2");
registry.addClone("eventQA/hPsi2FT0C", "hEPRectrV2");
registry.addClone("eventQA/hPsi2FT0C", "hEPTwistV2");
registry.addClone("eventQA/hQvecUncorV2", "eventQA/hQvecRectrV2");
registry.addClone("eventQA/hQvecUncorV2", "eventQA/hQvecTwistV2");
registry.addClone("eventQA/hQvecUncorV2", "eventQA/hQvecFinalV2");
registry.addClone("eventQA/hPsi2FT0C", "eventQA/hEPUncorV2");
registry.addClone("eventQA/hPsi2FT0C", "eventQA/hEPRectrV2");
registry.addClone("eventQA/hPsi2FT0C", "eventQA/hEPTwistV2");
}
if (doprocessESEBackground) {
LOGF(info, "JetSpectraEseTask::init() - Background Process");
Expand Down Expand Up @@ -1193,67 +1205,79 @@ struct JetSpectraEseTask {

static constexpr float InvalidValue = 999.;

// template <bool FillAllPsi, bool FillHist, typename EPCol>
template <EventPlaneFiller P, typename EPCol>
EventPlane procEP(EPCol const& vec)
{
constexpr std::array<float, 2> AmpCut{LowFT0Cut, 0.0};
auto computeEP = [&AmpCut](const std::vector<float>& vec, auto det, float n) { return vec[2] > AmpCut[det] ? (1.0 / n) * std::atan2(vec[1], vec[0]) : InvalidValue; };
std::map<std::string, float> epMap;
std::map<std::string, float> ep3Map;
auto vec1{qVecNoESE<DetID::FT0A, P.hist>(vec)};
epMap["FT0A"] = computeEP(vec1, 0, 2.0);
ep3Map["FT0A"] = computeEP(vec1, 0, 3.0);
auto vec2{qVecNoESE<DetID::FT0C, false>(vec)};
epMap["FT0C"] = computeEP(vec2, 0, 2.0);
ep3Map["FT0C"] = computeEP(vec2, 0, 3.0);
epMap["FV0A"] = computeEP(qVecNoESE<DetID::FV0A, false>(vec), 0, 2.0);
epMap["TPCpos"] = computeEP(qVecNoESE<DetID::TPCpos, false>(vec), 1, 2.0);
epMap["TPCneg"] = computeEP(qVecNoESE<DetID::TPCneg, false>(vec), 1, 2.0);
auto computeEP = [](const std::vector<float>& qVec, float minAmp, float harmonic) {
return qVec[2] > minAmp ? std::atan2(qVec[1], qVec[0]) / harmonic : InvalidValue;
};

std::map<std::string, float> epMap{
{"FT0A", computeEP(qVecNoESE<DetID::FT0A, P.hist>(vec, 2), LowFT0Cut, 2.0f)},
{"FT0C", computeEP(qVecNoESE<DetID::FT0C, false>(vec, 2), LowFT0Cut, 2.0f)},
{"FV0A", computeEP(qVecNoESE<DetID::FV0A, false>(vec), LowFT0Cut, 2.0f)},
{"TPCpos", computeEP(qVecNoESE<DetID::TPCpos, false>(vec), 0.0f, 2.0f)},
{"TPCneg", computeEP(qVecNoESE<DetID::TPCneg, false>(vec), 0.0f, 2.0f)}};
std::map<std::string, float> ep3Map{
{"FT0A", computeEP(qVecNoESE<DetID::FT0A, false>(vec, 3), LowFT0Cut, 3.0f)},
{"FT0C", computeEP(qVecNoESE<DetID::FT0C, false>(vec, 3), LowFT0Cut, 3.0f)}};
std::map<std::string, float> ep4Map{
{"FT0A", computeEP(qVecNoESE<DetID::FT0A, false>(vec, 4), LowFT0Cut, 4.0f)},
{"FT0C", computeEP(qVecNoESE<DetID::FT0C, false>(vec, 4), LowFT0Cut, 4.0f)},
{"FV0A", computeEP(qVecNoESE<DetID::FV0A, false>(vec, 4), LowFT0Cut, 4.0f)},
{"TPCpos", computeEP(qVecNoESE<DetID::TPCpos, false>(vec, 4), 0.0f, 4.0f)},
{"TPCneg", computeEP(qVecNoESE<DetID::TPCneg, false>(vec, 4), 0.0f, 4.0f)}};

if constexpr (P.psi) {
if constexpr (P.hist)
fillEP(/*std::make_index_sequence<5>{},*/ vec, epMap);
auto cosPsi = [](float psiX, float psiY) { return (static_cast<double>(psiX) == InvalidValue || static_cast<double>(psiY) == InvalidValue) ? InvalidValue : std::cos(2.0 * (psiX - psiY)); };
std::array<float, 3> epCorrContainer{};
epCorrContainer[0] = cosPsi(epMap.at(cfgEPRefA), epMap.at(cfgEPRefC));
epCorrContainer[1] = cosPsi(epMap.at(cfgEPRefA), epMap.at(cfgEPRefB));
epCorrContainer[2] = cosPsi(epMap.at(cfgEPRefB), epMap.at(cfgEPRefC));
if constexpr (P.hist)
fillEPCos(/*std::make_index_sequence<3>{},*/ vec, epCorrContainer);
}
EventPlane localPlane;
localPlane.psi2 = epMap.at(cfgEPRefA);
localPlane.psi3 = ep3Map.at(cfgEPRefA);
return localPlane;
// return epMap.at(cfgEPRefA);
if constexpr (P.hist) {
fillEP(vec, epMap, ep3Map, ep4Map);
}

auto cosPsi = [](float psiX, float psiY) {
return psiX == InvalidValue || psiY == InvalidValue ? InvalidValue : std::cos(2.0f * (psiX - psiY));
};
const std::array<float, 3> epCorrContainer{
cosPsi(epMap.at(cfgEPRefA), epMap.at(cfgEPRefC)),
cosPsi(epMap.at(cfgEPRefA), epMap.at(cfgEPRefB)),
cosPsi(epMap.at(cfgEPRefB), epMap.at(cfgEPRefC))};
const std::array<float, 3> epCorrContainer4{
cosPsi(ep4Map.at(cfgEPRefA), ep4Map.at(cfgEPRefC)),
cosPsi(ep4Map.at(cfgEPRefA), ep4Map.at(cfgEPRefB)),
cosPsi(ep4Map.at(cfgEPRefB), ep4Map.at(cfgEPRefC))};

if constexpr (P.hist) {
fillEPCos(vec, epCorrContainer, epCorrContainer4);
}
}
return {epMap.at(cfgEPRefA), ep3Map.at(cfgEPRefA)};
}
template </*std::size_t... Idx,*/ typename collision>
void fillEPCos(/*const std::index_sequence<Idx...>&,*/ const collision& col, const std::array<float, 3>& Corr)
template <typename collision>
void fillEPCos(const collision& col, const std::array<float, 3>& Corr, const std::array<float, 3>& Corr4)
{
// static constexpr std::string CosList[] = {"hCosPsi2AmC", "hCosPsi2AmB", "hCosPsi2BmC"};
// (registry.fill(HIST(CosList[Idx]), col.centrality(), Corr[Idx], col.qPERCFT0C()[0]), ...);
registry.fill(HIST("eventQA/hCosPsi2AmC"), col.centFT0M(), Corr[0], col.qPERCFT0C()[0]);
registry.fill(HIST("eventQA/hCosPsi2AmB"), col.centFT0M(), Corr[1], col.qPERCFT0C()[0]);
registry.fill(HIST("eventQA/hCosPsi2BmC"), col.centFT0M(), Corr[2], col.qPERCFT0C()[0]);

registry.fill(HIST("eventQA/hCosPsi4AmC"), col.centFT0M(), Corr4[0], col.qPERCFT0C()[0]);
registry.fill(HIST("eventQA/hCosPsi4AmB"), col.centFT0M(), Corr4[1], col.qPERCFT0C()[0]);
registry.fill(HIST("eventQA/hCosPsi4BmC"), col.centFT0M(), Corr4[2], col.qPERCFT0C()[0]);
}

template </*std::size_t... Idx,*/ typename collision>
void fillEP(/*const std::index_sequence<Idx...>&,*/ const collision& col, const std::map<std::string, float>& epMap)
template <typename collision>
void fillEP(const collision& col, const std::map<std::string, float>& epMap, const std::map<std::string, float>& ep3Map, const std::map<std::string, float>& ep4Map)
{
// static constexpr std::string_view EpList[] = {"hPsi2FT0A", "hPsi2FV0A", "hPsi2FT0C", "hPsi2TPCpos", "hPsi2TPCneg"};
// (registry.fill(HIST(EpList[Idx]), col.centrality(), epMap.at(std::string(RemovePrefix(EpList[Idx])))), ...);
registry.fill(HIST("eventQA/hPsi2FT0A"), col.centFT0M(), epMap.at("FT0A"));
registry.fill(HIST("eventQA/hPsi2FV0A"), col.centFT0M(), epMap.at("FV0A"));
registry.fill(HIST("eventQA/hPsi2FT0C"), col.centFT0M(), epMap.at("FT0C"));
registry.fill(HIST("eventQA/hPsi2TPCpos"), col.centFT0M(), epMap.at("TPCpos"));
registry.fill(HIST("eventQA/hPsi2TPCneg"), col.centFT0M(), epMap.at("TPCneg"));
}
constexpr std::string_view RemovePrefix(std::string_view str)
{
constexpr std::string_view Prefix{"hPsi2"};
return str.substr(Prefix.size());
}

registry.fill(HIST("eventQA/hPsi3FT0A"), col.centFT0M(), ep3Map.at("FT0A"));
registry.fill(HIST("eventQA/hPsi3FT0C"), col.centFT0M(), ep3Map.at("FT0C"));

registry.fill(HIST("eventQA/hPsi4FT0A"), col.centFT0M(), ep4Map.at("FT0A"));
registry.fill(HIST("eventQA/hPsi4FT0C"), col.centFT0M(), ep4Map.at("FT0C"));
}
constexpr int detIDN(const DetID id)
{
switch (id) {
Expand All @@ -1275,14 +1299,14 @@ struct JetSpectraEseTask {
return -1;
}

const int secondHarmonic{2};
template <DetID id, bool fill, typename Col>
std::vector<float> qVecNoESE(Col collision)
std::vector<float> qVecNoESE(Col collision, int nmode = 2)
{
// const int nmode{2};
int detId{detIDN(id)};
int detInd{detId * 4 /*+ cfgnTotalSystem * 4 * (nmode - 2)*/};
int detInd{detId * 4 + cfgnTotalSystem * 4 * (nmode - 2)};
if constexpr (fill) {
if (collision.qvecAmp()[detInd] > LowFT0Cut) {
if (collision.qvecAmp()[detInd] > LowFT0Cut && nmode == secondHarmonic) {
registry.fill(HIST("eventQA/hQvecUncorV2"), collision.centFT0M(), collision.qvecRe()[detInd], collision.qvecIm()[detInd]);
registry.fill(HIST("eventQA/hQvecRectrV2"), collision.centFT0M(), collision.qvecRe()[detInd + 1], collision.qvecIm()[detInd + 1]);
registry.fill(HIST("eventQA/hQvecTwistV2"), collision.centFT0M(), collision.qvecRe()[detInd + 2], collision.qvecIm()[detInd + 2]);
Expand Down
Loading