Skip to content
Merged
Changes from all commits
Commits
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
146 changes: 114 additions & 32 deletions PWGCF/Flow/Tasks/flowSP.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -89,6 +89,7 @@ struct FlowSP {
O2_DEFINE_CONFIGURABLE(cfgFillPIDQA, bool, false, "Fill histograms for PID QA");
O2_DEFINE_CONFIGURABLE(cfgFillEventPlaneQA, bool, false, "Fill histograms for Event Plane QA");
O2_DEFINE_CONFIGURABLE(cfgFillQABefore, bool, false, "Fill QA histograms before cuts, only for processData");
O2_DEFINE_CONFIGURABLE(cfgFillMeanPT, bool, false, "Fill histograms for mean PX/PT");
// Flags to make and fill histograms
O2_DEFINE_CONFIGURABLE(cfgFillGeneralV1Histos, bool, true, "Fill histograms for vn analysis");
O2_DEFINE_CONFIGURABLE(cfgFillMixedHarmonics, bool, true, "Flag to make and fill histos for mixed harmonics");
Expand Down Expand Up @@ -119,9 +120,10 @@ struct FlowSP {
O2_DEFINE_CONFIGURABLE(cfgFillWeights, bool, true, "Fill NUA weights");
O2_DEFINE_CONFIGURABLE(cfgFillWeightsPOS, bool, true, "Fill NUA weights only for positive charges");
O2_DEFINE_CONFIGURABLE(cfgFillWeightsNEG, bool, true, "Fill NUA weights only for negative charges");
O2_DEFINE_CONFIGURABLE(cfguseNUA1D, bool, false, "Use 1D NUA weights (only phi)");
O2_DEFINE_CONFIGURABLE(cfguseNUA2D, bool, true, "Use 2D NUA weights (phi and eta)");
O2_DEFINE_CONFIGURABLE(cfguseNUE2D, bool, true, "Use 2D NUE weights (pt and eta)");
O2_DEFINE_CONFIGURABLE(cfgUseNUA1D, bool, false, "Use 1D NUA weights (only phi)");
O2_DEFINE_CONFIGURABLE(cfgUseNUA2D, bool, true, "Use 2D NUA weights (phi and eta)");
O2_DEFINE_CONFIGURABLE(cfgUseNUE2D, bool, true, "Use 2D NUE weights");
O2_DEFINE_CONFIGURABLE(cfgUseNUE2Deta, bool, true, "Use 2D NUE weights TRUE: (pt and eta) FALSE: (pt and centrality)");
// Additional track Selections
O2_DEFINE_CONFIGURABLE(cfgTrackSelsUseAdditionalTrackCut, bool, false, "Bool to enable Additional Track Cut");
O2_DEFINE_CONFIGURABLE(cfgTrackSelsDoDCApt, bool, false, "Apply Pt dependent DCAz cut");
Expand All @@ -144,8 +146,8 @@ struct FlowSP {
ConfigurableAxis axisNch = {"axisNch", {400, 0, 4000}, "Global N_{ch}"};
ConfigurableAxis axisMultpv = {"axisMultpv", {400, 0, 4000}, "N_{ch} (PV)"};
// Configurables containing vector
Configurable<std::vector<double>> cfgEvSelsMultPv{"cfgEvSelsMultPv", std::vector<double>{2228.05, -75.5988, 0.976695, -0.00585275, 1.40738e-05, 3795.65, -136.988, 2.12393, -0.017028, 5.78679e-05}, "Multiplicity cuts (PV) first 5 parameters cutLOW last 5 cutHIGH (Default is +-2sigma pass5) "};
Configurable<std::vector<double>> cfgEvSelsMult{"cfgEvSelsMult", std::vector<double>{1308.86, -41.9314, 0.488423, -0.00248178, 4.71554e-06, 2973.55, -103.092, 1.47673, -0.0106685, 3.29348e-05}, "Multiplicity cuts (Global) first 5 parameters cutLOW last 5 cutHIGH (Default is +-2sigma pass5) "};
Configurable<std::vector<double>> cfgEvSelsMultPv{"cfgEvSelsMultPv", std::vector<double>{2223.49, -75.1444, 0.963572, -0.00570399, 1.34877e-05, 3790.99, -137.064, 2.13044, -0.017122, 5.82834e-05}, "Multiplicity cuts (PV) first 5 parameters cutLOW last 5 cutHIGH (Default is +-2sigma pass5) "};
Configurable<std::vector<double>> cfgEvSelsMult{"cfgEvSelsMult", std::vector<double>{1301.56, -41.4615, 0.478224, -0.00239449, 4.46966e-06, 2967.6, -102.927, 1.47488, -0.0106534, 3.28622e-05}, "Multiplicity cuts (Global) first 5 parameters cutLOW last 5 cutHIGH (Default is +-2sigma pass5) "};

Filter collisionFilter = nabs(aod::collision::posZ) < cfgEvSelsVtxZ;
Filter trackFilter = nabs(aod::track::eta) < cfgTrackSelsEta && aod::track::pt > cfgTrackSelsPtmin&& aod::track::pt < cfgTrackSelsPtmax && ((requireGlobalTrackInFilter()) || (aod::track::isGlobalTrackSDD == (uint8_t) true)) && nabs(aod::track::dcaXY) < cfgTrackSelsDCAxy&& nabs(aod::track::dcaZ) < cfgTrackSelsDCAz;
Expand Down Expand Up @@ -187,10 +189,12 @@ struct FlowSP {
TProfile* hcorrQQy = nullptr;
TProfile* hEvPlaneRes = nullptr;
TH1D* hCentrality = nullptr;
TProfile2D* hRelPt = nullptr;

bool clQQ = false;
bool clEvPlaneRes = false;
bool clCentrality = false;
bool clRelPt = false;

} cfg;

Expand Down Expand Up @@ -218,6 +222,12 @@ struct FlowSP {
double vy = 0;
double vz = 0;
int charge = 0;
double relPt = 1.;
double psiA = 0;
double psiC = 0;
double psiFull = 0;
double trackPxA = 0;
double trackPxC = 0;
} spm;

OutputObj<GFWWeights> fWeights{GFWWeights("weights")};
Expand Down Expand Up @@ -328,6 +338,7 @@ struct FlowSP {
AxisSpec axisNsigma = {100, -10, 10, "Nsigma for TPC and TOF"};
AxisSpec axisdEdx = {300, 0, 300, "dEdx for PID"};
AxisSpec axisBeta = {150, 0, 1.5, "Beta for PID"};
AxisSpec axisCharge = {3, 0, 3, "Charge: 0 = inclusive, 1 = positive, 2 = negative"};

std::vector<double> ptbinning = {0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2, 2.2, 2.4, 2.6, 2.8, 3, 3.5, 4, 5, 6, 8, 10};
AxisSpec axisPt = {ptbinning, "#it{p}_{T} GeV/#it{c}"};
Expand Down Expand Up @@ -428,6 +439,7 @@ struct FlowSP {
histos.add("incl/QA/after/hSharedClusters_pt", "", {HistType::kTH2D, {axisPt, axisShCl}});
histos.add("incl/QA/after/hCrossedRows_pt", "", {HistType::kTH2D, {axisPt, axisCl}});
histos.add("incl/QA/after/hCrossedRows_vs_SharedClusters", "", {HistType::kTH2D, {axisCl, axisShCl}});
histos.add("incl/QA/after/hMeanPtEta", "", {HistType::kTProfile2D, {axisEta, axisCent}});

if (cfgTrackSelDoTrackQAvsCent) {
histos.add<TH3>("incl/QA/after/hPt_Eta", "", kTH3D, {axisPt, axisEta, axisCent});
Expand Down Expand Up @@ -511,6 +523,10 @@ struct FlowSP {
registry.add<TProfile3D>("incl/vnC", "", kTProfile3D, {axisPt, axisEtaVn, axisCentrality});
registry.add<TProfile3D>("incl/vnA", "", kTProfile3D, {axisPt, axisEtaVn, axisCentrality});
}
if (cfgFillMeanPT) {
registry.add<TProfile2D>("incl/meanPT/meanRelPtA", "", kTProfile2D, {axisEtaVn, axisCentrality});
registry.add<TProfile2D>("incl/meanPT/meanRelPtC", "", kTProfile2D, {axisEtaVn, axisCentrality});
}
if (cfgFillPID) {
registry.add<TProfile3D>("incl/pion/vnC", "", kTProfile3D, {axisPt, axisEtaVn, axisCentrality});
registry.add<TProfile3D>("incl/pion/vnA", "", kTProfile3D, {axisPt, axisEtaVn, axisCentrality});
Expand Down Expand Up @@ -710,7 +726,7 @@ struct FlowSP {

int nWeights = 3;

if (cfguseNUA1D) {
if (cfgUseNUA1D) {
if (cfgCCDB_NUA.value.empty() == false) {
TList* listCorrections = ccdb->getForTimeStamp<TList>(cfgCCDB_NUA, timestamp);
cfg.mAcceptance.push_back(reinterpret_cast<GFWWeights*>(listCorrections->FindObject("weights")));
Expand All @@ -724,7 +740,7 @@ struct FlowSP {
} else {
LOGF(info, "cfgCCDB_NUA empty! No corrections loaded");
}
} else if (cfguseNUA2D) {
} else if (cfgUseNUA2D) {
if (cfgCCDB_NUA.value.empty() == false) {
TH3D* hNUA2D = ccdb->getForTimeStamp<TH3D>(cfgCCDB_NUA, timestamp);
if (!hNUA2D) {
Expand Down Expand Up @@ -754,9 +770,9 @@ struct FlowSP {
// Get Efficiency correction
if (cfgCCDB_NUE2D.value.empty() == false) {
TList* listCorrections = ccdb->getForTimeStamp<TList>(cfgCCDB_NUE2D, timestamp);
cfg.mEfficiency2D.push_back(reinterpret_cast<TH2D*>(listCorrections->FindObject("Efficiency2D")));
cfg.mEfficiency2D.push_back(reinterpret_cast<TH2D*>(listCorrections->FindObject("Efficiency2D_pos")));
cfg.mEfficiency2D.push_back(reinterpret_cast<TH2D*>(listCorrections->FindObject("Efficiency2D_neg")));
cfg.mEfficiency2D.push_back(reinterpret_cast<TH2D*>(listCorrections->FindObject("Efficiency")));
cfg.mEfficiency2D.push_back(reinterpret_cast<TH2D*>(listCorrections->FindObject("Efficiency_pos")));
cfg.mEfficiency2D.push_back(reinterpret_cast<TH2D*>(listCorrections->FindObject("Efficiency_neg")));
int sizeEff = cfg.mEfficiency2D.size();
if (sizeEff < nWeights)
LOGF(fatal, "Could not load efficiency histogram for trigger particles from %s", cfgCCDB_NUE.value.c_str());
Expand All @@ -769,14 +785,21 @@ struct FlowSP {
}

// From Generic Framework
bool setCurrentParticleWeights(int pID, int spec, const float& phi, const float& eta, const float& pt, const float& vtxz)
bool setCurrentParticleWeights(int pID, int spec, const float& phi, const float& eta, const float& pt, const float& vtxz, const float& centrality)
{
float eff = 1.;
int sizeEff = cfg.mEfficiency.size();
if (sizeEff > pID) {
if (cfguseNUE2D) {
int binx = cfg.mEfficiency2D[pID]->GetXaxis()->FindBin(eta);
int biny = cfg.mEfficiency2D[pID]->GetYaxis()->FindBin(pt);
if (cfgUseNUE2D) {
int binx;
int biny;
if (cfgUseNUE2Deta) {
biny = cfg.mEfficiency2D[pID]->GetYaxis()->FindBin(pt);
binx = cfg.mEfficiency2D[pID]->GetXaxis()->FindBin(eta);
} else {
binx = cfg.mEfficiency2D[pID]->GetXaxis()->FindBin(pt);
biny = cfg.mEfficiency2D[pID]->GetYaxis()->FindBin(centrality);
}
eff = cfg.mEfficiency2D[pID]->GetBinContent(binx, biny);
} else {
eff = cfg.mEfficiency[pID]->GetBinContent(cfg.mEfficiency[pID]->FindBin(pt));
Expand All @@ -789,14 +812,14 @@ struct FlowSP {

spm.weff[pID][spec] = 1. / eff;

if (cfguseNUA1D) {
if (cfgUseNUA1D) {
int sizeAcc = cfg.mAcceptance.size();
if (sizeAcc > pID) {
spm.wacc[pID][spec] = cfg.mAcceptance[pID]->getNUA(phi, eta, vtxz);
} else {
spm.wacc[pID][spec] = 1;
}
} else if (cfguseNUA2D) {
} else if (cfgUseNUA2D) {
if (cfg.mAcceptance2D.size() > 0) {
spm.wacc[pID][spec] = getNUA2D(cfg.mAcceptance2D[0], eta, phi, vtxz);
} else {
Expand Down Expand Up @@ -1068,6 +1091,7 @@ struct FlowSP {
histos.fill(HIST(Charge[ct]) + HIST(Species[par]) + HIST("QA/") + HIST(Time[ft]) + HIST("hSharedClusters_pt"), track.pt(), track.tpcFractionSharedCls(), spm.wacc[ct][par] * spm.weff[ct][par]);
histos.fill(HIST(Charge[ct]) + HIST(Species[par]) + HIST("QA/") + HIST(Time[ft]) + HIST("hCrossedRows_pt"), track.pt(), track.tpcNClsFound(), spm.wacc[ct][par] * spm.weff[ct][par]);
histos.fill(HIST(Charge[ct]) + HIST(Species[par]) + HIST("QA/") + HIST(Time[ft]) + HIST("hCrossedRows_vs_SharedClusters"), track.tpcNClsFound(), track.tpcFractionSharedCls(), spm.wacc[ct][par] * spm.weff[ct][par]);
histos.fill(HIST(Charge[ct]) + HIST(Species[par]) + HIST("QA/") + HIST(Time[ft]) + HIST("hMeanPtEta"), track.eta(), spm.centrality, track.pt(), spm.wacc[ct][par] * spm.weff[ct][par]);
}

template <FillType ft, ChargeType ct, typename TrackObject>
Expand Down Expand Up @@ -1204,11 +1228,11 @@ struct FlowSP {
spm.vz = collision.posZ();
float vtxz = collision.posZ();

double psiA = 1.0 * std::atan2(spm.qyA, spm.qxA);
double psiC = 1.0 * std::atan2(spm.qyC, spm.qxC);
spm.psiA = 1.0 * std::atan2(spm.qyA, spm.qxA);
spm.psiC = 1.0 * std::atan2(spm.qyC, spm.qxC);

// https://twiki.cern.ch/twiki/pub/ALICE/DirectedFlowAnalysisNote/vn_ZDC_ALICE_INT_NOTE_version02.pdf
double psiFull = 1.0 * std::atan2(spm.qyA + spm.qyC, spm.qxA + spm.qxC);
spm.psiFull = 1.0 * std::atan2(spm.qyA + spm.qyC, spm.qxA + spm.qxC);

// always fill these histograms!
registry.fill(HIST("QQCorrelations/qAqCXY"), spm.centrality, spm.qxA * spm.qxC + spm.qyA * spm.qyC);
Expand All @@ -1222,14 +1246,14 @@ struct FlowSP {
histos.fill(HIST("QA/hCentFull"), spm.centrality, 1);
}
if (cfgFillEventPlaneQA) {
histos.fill(HIST("QA/hSPplaneA"), psiA, 1);
histos.fill(HIST("QA/hSPplaneC"), psiC, 1);
histos.fill(HIST("QA/hSPplaneFull"), psiFull, 1);
histos.fill(HIST("QA/hCosPhiACosPhiC"), spm.centrality, std::cos(psiA) * std::cos(psiC));
histos.fill(HIST("QA/hSinPhiASinPhiC"), spm.centrality, std::sin(psiA) * std::sin(psiC));
histos.fill(HIST("QA/hSinPhiACosPhiC"), spm.centrality, std::sin(psiA) * std::cos(psiC));
histos.fill(HIST("QA/hCosPhiASinsPhiC"), spm.centrality, std::cos(psiA) * std::sin(psiC));
histos.fill(HIST("QA/hFullEvPlaneRes"), spm.centrality, -1 * std::cos(psiA - psiC));
histos.fill(HIST("QA/hSPplaneA"), spm.psiA, 1);
histos.fill(HIST("QA/hSPplaneC"), spm.psiC, 1);
histos.fill(HIST("QA/hSPplaneFull"), spm.psiFull, 1);
histos.fill(HIST("QA/hCosPhiACosPhiC"), spm.centrality, std::cos(spm.psiA) * std::cos(spm.psiC));
histos.fill(HIST("QA/hSinPhiASinPhiC"), spm.centrality, std::sin(spm.psiA) * std::sin(spm.psiC));
histos.fill(HIST("QA/hSinPhiACosPhiC"), spm.centrality, std::sin(spm.psiA) * std::cos(spm.psiC));
histos.fill(HIST("QA/hCosPhiASinsPhiC"), spm.centrality, std::cos(spm.psiA) * std::sin(spm.psiC));
histos.fill(HIST("QA/hFullEvPlaneRes"), spm.centrality, -1 * std::cos(spm.psiA - spm.psiC));
}

if (spm.centrality > cfgCentMax || spm.centrality < cfgCentMin)
Expand Down Expand Up @@ -1278,6 +1302,10 @@ struct FlowSP {

fillEventQA<kAfter>(collision, tracks);

TProfile2D* hRelEtaPt = new TProfile2D("hRelEtaPt", "hRelEtaPt", 8, -.8, .8, 3, 0, 3);
TProfile2D* ptV1A = new TProfile2D("ptV1A", "ptV1A", 8, -.8, .8, 3, 0, 3);
TProfile2D* ptV1C = new TProfile2D("ptV1C", "ptV1C", 8, -.8, .8, 3, 0, 3);

for (const auto& track : tracks) {

ParticleType trackPID = (cfgFillPID || cfgFillPIDQA) ? getTrackPID(track) : kUnidentified;
Expand Down Expand Up @@ -1333,13 +1361,25 @@ struct FlowSP {
}

// Set weff and wacc for inclusive, negative and positive hadrons
if (!setCurrentParticleWeights(kInclusive, kUnidentified, phi, track.eta(), track.pt(), vtxz))
if (!setCurrentParticleWeights(kInclusive, kUnidentified, phi, track.eta(), track.pt(), vtxz, spm.centrality))
continue;
if (!setCurrentParticleWeights(spm.charge, kUnidentified, phi, track.eta(), track.pt(), vtxz))
if (!setCurrentParticleWeights(spm.charge, kUnidentified, phi, track.eta(), track.pt(), vtxz, spm.centrality))
continue;

histos.fill(HIST("hTrackCount"), trackSel_ParticleWeights);

double weight = spm.wacc[0][0] * spm.weff[0][0] * spm.centWeight;
double weight_charged = spm.wacc[spm.charge][0] * spm.weff[spm.charge][0] * spm.centWeight;

hRelEtaPt->Fill(track.eta(), kInclusive, track.pt(), weight);
hRelEtaPt->Fill(track.eta(), spm.charge, track.pt(), weight_charged);

ptV1A->Fill(track.eta(), kInclusive, track.pt() * ((spm.uy * spm.qyA + spm.ux * spm.qxA) / std::sqrt(std::fabs(spm.corrQQ))), weight);
ptV1A->Fill(track.eta(), spm.charge, track.pt() * ((spm.uy * spm.qyA + spm.ux * spm.qxA) / std::sqrt(std::fabs(spm.corrQQ))), weight_charged);

ptV1C->Fill(track.eta(), kInclusive, track.pt() * ((spm.uy * spm.qyC + spm.ux * spm.qxC) / std::sqrt(std::fabs(spm.corrQQ))), weight);
ptV1C->Fill(track.eta(), spm.charge, track.pt() * ((spm.uy * spm.qyC + spm.ux * spm.qxC) / std::sqrt(std::fabs(spm.corrQQ))), weight_charged);

fillAllQA<kAfter, kUnidentified>(track);
if (cfgFillPIDQA) {
switch (trackPID) {
Expand All @@ -1365,9 +1405,11 @@ struct FlowSP {
spm.uyMH = std::sin(cfgHarmMixed * phi);
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

spm.vnA = std::cos(cfgHarm * (phi - psiA)) / evPlaneRes;
spm.vnC = std::cos(cfgHarm * (phi - psiC)) / evPlaneRes;
spm.vnFull = std::cos(cfgHarm * (phi - psiFull)) / evPlaneRes;
spm.vnA = std::cos(cfgHarm * (phi - spm.psiA)) / evPlaneRes;
spm.vnC = std::cos(cfgHarm * (phi - spm.psiC)) / evPlaneRes;
spm.vnFull = std::cos(cfgHarm * (phi - spm.psiFull)) / evPlaneRes;

// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

fillHistograms<kInclusive, kUnidentified>(track);

Expand Down Expand Up @@ -1435,6 +1477,46 @@ struct FlowSP {
} // end of fillPID

} // end of track loop

// Now we want to fill the final relPt histogram
// Loop over all eta and fill bins
if (cfgFillMeanPT) {
for (int i = 0; i < hRelEtaPt->GetNbinsX(); i++) {
double eta = hRelEtaPt->GetXaxis()->GetBinCenter(i);

int bin = hRelEtaPt->FindBin(eta, kInclusive);

double drelPt = hRelEtaPt->GetBinContent(bin);
double dptV1A = ptV1A->GetBinContent(bin);
double dptV1C = ptV1C->GetBinContent(bin);
if (drelPt)
registry.fill(HIST("incl/meanPT/meanRelPtA"), eta, spm.centrality, dptV1A / drelPt, 1);
if (drelPt)
registry.fill(HIST("incl/meanPT/meanRelPtC"), eta, spm.centrality, dptV1C / drelPt, 1);

bin = hRelEtaPt->FindBin(eta, kPositive);
double drelPt_pos = hRelEtaPt->GetBinContent(bin);
double dptV1A_pos = ptV1A->GetBinContent(bin);
double dptV1C_pos = ptV1C->GetBinContent(bin);
if (drelPt_pos)
registry.fill(HIST("pos/meanPT/meanRelPtA"), eta, spm.centrality, dptV1A_pos / drelPt_pos, 1);
if (drelPt_pos)
registry.fill(HIST("pos/meanPT/meanRelPtC"), eta, spm.centrality, dptV1C_pos / drelPt_pos, 1);

bin = hRelEtaPt->FindBin(eta, kNegative);
double drelPt_neg = hRelEtaPt->GetBinContent(bin);
double dptV1A_neg = ptV1A->GetBinContent(bin);
double dptV1C_neg = ptV1C->GetBinContent(bin);
if (drelPt_neg)
registry.fill(HIST("neg/meanPT/meanRelPtA"), eta, spm.centrality, dptV1A_neg / drelPt_neg, 1);
if (drelPt_neg)
registry.fill(HIST("neg/meanPT/meanRelPtC"), eta, spm.centrality, dptV1C_neg / drelPt_neg, 1);
}
}

delete hRelEtaPt;
delete ptV1A;
delete ptV1C;
}

PROCESS_SWITCH(FlowSP, processData, "Process analysis for non-derived data", true);
Expand Down
Loading