diff --git a/PWGCF/Flow/Tasks/flowSP.cxx b/PWGCF/Flow/Tasks/flowSP.cxx index 4ea5937776d..c7b1a28b004 100644 --- a/PWGCF/Flow/Tasks/flowSP.cxx +++ b/PWGCF/Flow/Tasks/flowSP.cxx @@ -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"); @@ -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"); @@ -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> cfgEvSelsMultPv{"cfgEvSelsMultPv", std::vector{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> cfgEvSelsMult{"cfgEvSelsMult", std::vector{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> cfgEvSelsMultPv{"cfgEvSelsMultPv", std::vector{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> cfgEvSelsMult{"cfgEvSelsMult", std::vector{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; @@ -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; @@ -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 fWeights{GFWWeights("weights")}; @@ -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 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}"}; @@ -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("incl/QA/after/hPt_Eta", "", kTH3D, {axisPt, axisEta, axisCent}); @@ -511,6 +523,10 @@ struct FlowSP { registry.add("incl/vnC", "", kTProfile3D, {axisPt, axisEtaVn, axisCentrality}); registry.add("incl/vnA", "", kTProfile3D, {axisPt, axisEtaVn, axisCentrality}); } + if (cfgFillMeanPT) { + registry.add("incl/meanPT/meanRelPtA", "", kTProfile2D, {axisEtaVn, axisCentrality}); + registry.add("incl/meanPT/meanRelPtC", "", kTProfile2D, {axisEtaVn, axisCentrality}); + } if (cfgFillPID) { registry.add("incl/pion/vnC", "", kTProfile3D, {axisPt, axisEtaVn, axisCentrality}); registry.add("incl/pion/vnA", "", kTProfile3D, {axisPt, axisEtaVn, axisCentrality}); @@ -710,7 +726,7 @@ struct FlowSP { int nWeights = 3; - if (cfguseNUA1D) { + if (cfgUseNUA1D) { if (cfgCCDB_NUA.value.empty() == false) { TList* listCorrections = ccdb->getForTimeStamp(cfgCCDB_NUA, timestamp); cfg.mAcceptance.push_back(reinterpret_cast(listCorrections->FindObject("weights"))); @@ -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(cfgCCDB_NUA, timestamp); if (!hNUA2D) { @@ -754,9 +770,9 @@ struct FlowSP { // Get Efficiency correction if (cfgCCDB_NUE2D.value.empty() == false) { TList* listCorrections = ccdb->getForTimeStamp(cfgCCDB_NUE2D, timestamp); - cfg.mEfficiency2D.push_back(reinterpret_cast(listCorrections->FindObject("Efficiency2D"))); - cfg.mEfficiency2D.push_back(reinterpret_cast(listCorrections->FindObject("Efficiency2D_pos"))); - cfg.mEfficiency2D.push_back(reinterpret_cast(listCorrections->FindObject("Efficiency2D_neg"))); + cfg.mEfficiency2D.push_back(reinterpret_cast(listCorrections->FindObject("Efficiency"))); + cfg.mEfficiency2D.push_back(reinterpret_cast(listCorrections->FindObject("Efficiency_pos"))); + cfg.mEfficiency2D.push_back(reinterpret_cast(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()); @@ -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)); @@ -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 { @@ -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 @@ -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); @@ -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) @@ -1278,6 +1302,17 @@ struct FlowSP { fillEventQA(collision, tracks); + float trackPt = 0; + float trackPx = 0; + int trackSize = 0; + + 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); + + double sumPt = 0; + int ptCounter = 0; + for (const auto& track : tracks) { ParticleType trackPID = (cfgFillPID || cfgFillPIDQA) ? getTrackPID(track) : kUnidentified; @@ -1333,13 +1368,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(track); if (cfgFillPIDQA) { switch (trackPID) { @@ -1365,9 +1412,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(track); @@ -1435,6 +1484,47 @@ 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 nEntries = hRelEtaPt->GetEntries(); + + 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);