From 7b593bd7b2c1f4a978535f292f3d94b355653adf Mon Sep 17 00:00:00 2001 From: losiash Date: Fri, 8 May 2026 19:47:08 +0200 Subject: [PATCH] PWGCF: update pp ESE flow fluctuation task --- PWGCF/Flow/Tasks/flowFlucGfwPp.cxx | 168 ++++++++++++++++++++++++----- 1 file changed, 141 insertions(+), 27 deletions(-) diff --git a/PWGCF/Flow/Tasks/flowFlucGfwPp.cxx b/PWGCF/Flow/Tasks/flowFlucGfwPp.cxx index 1405b0a60b0..14de3f9093e 100644 --- a/PWGCF/Flow/Tasks/flowFlucGfwPp.cxx +++ b/PWGCF/Flow/Tasks/flowFlucGfwPp.cxx @@ -11,7 +11,7 @@ /// \file flowFlucGfwPp.cxx /// \brief GFW task for Event Shape Engineering studies in pp collisions -/// \author Emil Gorm Nielsen, NBI, emil.gorm.nielsen@cern.ch, Wenya Wu, TUM, wenya.wu@cern.ch +/// \author Wenya Wu, TUM, wenya.wu@cern.ch, David Krylenkov, TUM #include "PWGCF/GenericFramework/Core/FlowContainer.h" #include "PWGCF/GenericFramework/Core/GFW.h" @@ -23,7 +23,6 @@ #include "Common/DataModel/Centrality.h" #include "Common/DataModel/EventSelection.h" #include "Common/DataModel/Multiplicity.h" -#include "Common/DataModel/Qvectors.h" #include "Common/DataModel/TrackSelectionTables.h" #include @@ -106,6 +105,7 @@ std::vector firstRunsOfFill; struct FlowFlucGfwPp { static constexpr int kInvalidQnBin = -999; static constexpr float kInvalidQnSeparator = -999.f; + static constexpr float kTPCSubeventEtaGapHalfWidth = 0.1f; static constexpr int kRequireBothEtaSides = 1; static constexpr int kRequireFullFourParticleTracks = 2; @@ -120,6 +120,7 @@ struct FlowFlucGfwPp { O2_DEFINE_CONFIGURABLE(cfgIsMC, bool, false, "Is MC event") O2_DEFINE_CONFIGURABLE(cfgCentEstimator, int, 0, "0:FT0C; 1:FT0CVariant1; 2:FT0M; 3:FV0A, 4:NTPV, 5:NGlobals, 6:MFT") O2_DEFINE_CONFIGURABLE(cfgUseNch, bool, false, "Do correlations as function of Nch") + O2_DEFINE_CONFIGURABLE(cfgQvecQA, bool, false, "Enable filling QA for q-Vec of TPC") O2_DEFINE_CONFIGURABLE(cfgFillWeights, bool, false, "Fill NUA weights") O2_DEFINE_CONFIGURABLE(cfgRunByRun, bool, false, "Fill histograms on a run-by-run basis") O2_DEFINE_CONFIGURABLE(cfgFillFlowRunByRun, bool, false, "Fill flow profile run-by-run (only for v22)") @@ -165,6 +166,10 @@ struct FlowFlucGfwPp { O2_DEFINE_CONFIGURABLE(cfgCentMax, float, 100, "Maximum of centrality or multiplicity"); O2_DEFINE_CONFIGURABLE(cfgEvtSelCent, bool, true, "Choose event selector as centrality(true) or multicplity(false)"); O2_DEFINE_CONFIGURABLE(cfgUseNegativeEtaHalfForq2, bool, true, "If true, use -eta half for q2 selection; otherwise use +eta half"); + O2_DEFINE_CONFIGURABLE(cfgMinPtOnTPC, float, 0.2, "minimum transverse momentum selection for TPC tracks participating in Q-vector reconstruction"); + O2_DEFINE_CONFIGURABLE(cfgMaxPtOnTPC, float, 5., "maximum transverse momentum selection for TPC tracks participating in Q-vector reconstruction"); + O2_DEFINE_CONFIGURABLE(cfgEtaMax, float, 0.8, "Maximum pseudorapidiy for charged track"); + O2_DEFINE_CONFIGURABLE(cfgEtaMin, float, -0.8, "Minimum pseudorapidiy for charged track"); Configurable> qnBinSeparator{"qnBinSeparator", std::vector{-999.f, -999.f, -999.f}, "Qn bin separator"}; struct : ConfigurableGroup { @@ -325,11 +330,26 @@ struct FlowFlucGfwPp { o2::framework::expressions::Filter mcCollFilter = nabs(aod::mccollision::posZ) < cfgVtxZ; o2::framework::expressions::Filter mcParticlesFilter = (aod::mcparticle::eta > o2::analysis::gfwflowflucpp::etalow && aod::mcparticle::eta < o2::analysis::gfwflowflucpp::etaup && aod::mcparticle::pt > o2::analysis::gfwflowflucpp::ptlow && aod::mcparticle::pt < o2::analysis::gfwflowflucpp::ptup); - using GFWTracks = soa::Filtered>; + using GFWTracks = soa::Filtered>; void init(InitContext const&) { LOGF(info, "FlowFlucGfwPp::init()"); + if (static_cast(cfgMinPtOnTPC) < static_cast(cfgPtmin) || + static_cast(cfgMaxPtOnTPC) > static_cast(cfgPtmax) || + static_cast(cfgEtaMax) > static_cast(cfgEta) || + static_cast(cfgEtaMin) < -static_cast(cfgEta)) { + LOGF(warning, + "The Q-vector loop sees only tracks that passed the main trackFilter. " + "[pt %.3g, %.3g], [eta %.3g, %.3g] and input cuts [pt %.3g, %.3g], |eta|<%.3g", + static_cast(cfgMinPtOnTPC), + static_cast(cfgMaxPtOnTPC), + static_cast(cfgEtaMin), + static_cast(cfgEtaMax), + static_cast(cfgPtmin), + static_cast(cfgPtmax), + static_cast(cfgEta)); + } o2::analysis::gfwflowflucpp::regions.SetNames(cfgRegions->GetNames()); o2::analysis::gfwflowflucpp::regions.SetEtaMin(cfgRegions->GetEtaMin()); o2::analysis::gfwflowflucpp::regions.SetEtaMax(cfgRegions->GetEtaMax()); @@ -427,6 +447,19 @@ struct FlowFlucGfwPp { int ptbins = o2::analysis::gfwflowflucpp::ptbinning.size() - 1; fSecondAxis = (cfgTimeDependent) ? new TAxis(timeAxis.binEdges.size() - 1, &(timeAxis.binEdges[0])) : new TAxis(ptbins, &o2::analysis::gfwflowflucpp::ptbinning[0]); + if (cfgQvecQA && (doprocessData || doprocessq2)) { + // qVectorsTable-equivalent TPC-track QA for the in-task raw TPC Q-vector loop. + AxisSpec qVecAxisPt = {40, 0.0, 4.0}; + AxisSpec qVecAxisEta = {32, -0.8, 0.8}; + AxisSpec qVecAxisPhi = {32, 0, constants::math::TwoPI}; + AxisSpec qVecAxisCent = {20, 0, 100}; + AxisSpec qVecAxisMulti = {20, 0, 1000}; + AxisSpec qVecAxis = {20, -10, 10}; + registry.add("qvecQA/ChTracks", ";pT;#eta;#phi", {HistType::kTHnSparseF, {qVecAxisPt, qVecAxisEta, qVecAxisPhi}}); + registry.add("qvecQA/CountEvt", ";Centrality;TrkSize;TrkSel;MultNTrkPV", {HistType::kTHnSparseF, {qVecAxisCent, qVecAxisMulti, qVecAxisMulti, qVecAxisMulti}}); + registry.add("qvecQA/QxQy", ";QxAll;QyAll;QxNeg;QyNeg;QxPos;QyPos", {HistType::kTHnSparseF, {qVecAxis, qVecAxis, qVecAxis, qVecAxis, qVecAxis, qVecAxis}}); + } + if (doprocessq2) { registry.add("mq2/eventcounter", "", HistType::kTH1F, {{10, 0, 10}}); registry.add("mq2/h2_cent_q2_etapos", ";Centrality;#it{q}_{2}^{#eta pos};", HistType::kTH2D, {{100, 0, 100}, {600, 0, 6}}); @@ -884,28 +917,108 @@ struct FlowFlucGfwPp { int nMid; }; - template - float computeqnVec(T const& col, bool useNegativeEtaHalf) + struct InTaskTPCQVector { + float qVectTPCPos[2] = {0.f, 0.f}; // Always computed + float qVectTPCNeg[2] = {0.f, 0.f}; // Always computed + float qVectTPCAll[2] = {0.f, 0.f}; // Always computed + + int nTrkTPCPos = 0; + int nTrkTPCNeg = 0; + int nTrkTPCAll = 0; + + std::vector trkTPCPosLabel{}; + std::vector trkTPCNegLabel{}; + std::vector trkTPCAllLabel{}; + }; + + template + bool selTrack(const TrackType track) { - if (col.qvecTPCposReVec().empty() || col.qvecTPCposImVec().empty() || - col.qvecTPCnegReVec().empty() || col.qvecTPCnegImVec().empty()) { - return -1.f; + if (track.pt() < cfgMinPtOnTPC) + return false; + if (track.pt() > cfgMaxPtOnTPC) + return false; + if (!track.passedITSNCls()) + return false; + if (!track.passedITSChi2NDF()) + return false; + if (!track.passedITSHits()) + return false; + if (!track.passedTPCCrossedRowsOverNCls()) + return false; + if (!track.passedTPCChi2NDF()) + return false; + if (!track.passedDCAxy()) + return false; + if (!track.passedDCAz()) + return false; + + return true; + } + + template + InTaskTPCQVector calcQVec(const Nmode nMode, const TrackType& tracks, const TCollision& collision) + { + InTaskTPCQVector qvec; + + double nTrkSel = 0.; + for (auto const& trk : tracks) { + if (!selTrack(trk)) { + continue; + } + if (trk.eta() > cfgEtaMax) { + continue; + } + if (trk.eta() < cfgEtaMin) { + continue; + } + qvec.qVectTPCAll[0] += trk.pt() * std::cos(trk.phi() * nMode); + qvec.qVectTPCAll[1] += trk.pt() * std::sin(trk.phi() * nMode); + qvec.trkTPCAllLabel.push_back(trk.globalIndex()); + qvec.nTrkTPCAll++; + nTrkSel++; + if (std::abs(trk.eta()) < kTPCSubeventEtaGapHalfWidth) { + continue; + } + if (cfgQvecQA) { + registry.fill(HIST("qvecQA/ChTracks"), trk.pt(), trk.eta(), trk.phi()); + } + + if (trk.eta() > 0) { + // In qVectorsTable this branch is additionally guarded by useDetector["QvectorTPCposs"] || useDetector["QvectorBPoss"]. + // Here TPCpos is always computed because the downstream ESE selector can require it. + qvec.qVectTPCPos[0] += trk.pt() * std::cos(trk.phi() * nMode); + qvec.qVectTPCPos[1] += trk.pt() * std::sin(trk.phi() * nMode); + qvec.trkTPCPosLabel.push_back(trk.globalIndex()); + qvec.nTrkTPCPos++; + } else if (trk.eta() < 0) { + // In qVectorsTable this branch is additionally guarded by useDetector["QvectorTPCnegs"] || useDetector["QvectorBNegs"]. + // Here TPCneg is always computed because the downstream ESE selector can require it. + qvec.qVectTPCNeg[0] += trk.pt() * std::cos(trk.phi() * nMode); + qvec.qVectTPCNeg[1] += trk.pt() * std::sin(trk.phi() * nMode); + qvec.trkTPCNegLabel.push_back(trk.globalIndex()); + qvec.nTrkTPCNeg++; + } } - if (col.nTrkTPCpos() <= 0 || col.nTrkTPCneg() <= 0) - return -1.f; + if (cfgQvecQA) { + registry.fill(HIST("qvecQA/CountEvt"), collision.centFT0C(), tracks.size(), nTrkSel, collision.multNTracksPV()); + registry.fill(HIST("qvecQA/QxQy"), qvec.qVectTPCAll[0], qvec.qVectTPCAll[1], qvec.qVectTPCNeg[0], qvec.qVectTPCNeg[1], qvec.qVectTPCPos[0], qvec.qVectTPCPos[1]); + } - const auto qvecPos = - std::sqrt(col.qvecTPCposReVec()[0] * col.qvecTPCposReVec()[0] + - col.qvecTPCposImVec()[0] * col.qvecTPCposImVec()[0]) * - std::sqrt(col.nTrkTPCpos()); + return qvec; + } - const auto qvecNeg = - std::sqrt(col.qvecTPCnegReVec()[0] * col.qvecTPCnegReVec()[0] + - col.qvecTPCnegImVec()[0] * col.qvecTPCnegImVec()[0]) * - std::sqrt(col.nTrkTPCneg()); + float computeqnVec(InTaskTPCQVector const& qvec, bool useNegativeEtaHalf) + { + if (qvec.nTrkTPCPos <= 0 || qvec.nTrkTPCNeg <= 0) { + return -1.f; + } - return useNegativeEtaHalf ? qvecNeg : qvecPos; + if (useNegativeEtaHalf) { + return std::hypot(qvec.qVectTPCNeg[0], qvec.qVectTPCNeg[1]) / std::sqrt(static_cast(qvec.nTrkTPCNeg)); + } + return std::hypot(qvec.qVectTPCPos[0], qvec.qVectTPCPos[1]) / std::sqrt(static_cast(qvec.nTrkTPCPos)); } /// \return the 1-d qn-vector separator to 2-d @@ -1230,8 +1343,7 @@ struct FlowFlucGfwPp { void processData(soa::Filtered>::iterator const& collision, + aod::CentMFTs>>::iterator const& collision, aod::BCsWithTimestamps const&, GFWTracks const& tracks) { auto bc = collision.bc_as(); @@ -1275,7 +1387,7 @@ struct FlowFlucGfwPp { const XAxis xaxis{ getCentrality(collision), - tracks.size(), + collision.multNTracksPV(), (cfgTimeDependent) ? getTimeSinceStartOfFill(bc.timestamp(), *firstRunOfCurrentFill) : -1.0}; if (cfgTimeDependent && run == *firstRunOfCurrentFill && @@ -1294,7 +1406,8 @@ struct FlowFlucGfwPp { if (cfgFillQA) fillEventQA(collision, xaxis); - float qn = computeqnVec(collision, cfgUseNegativeEtaHalfForq2); + const auto qvecTPC = calcQVec(2, tracks, collision); + float qn = computeqnVec(qvecTPC, cfgUseNegativeEtaHalfForq2); if (qn < 0) return; @@ -1307,7 +1420,7 @@ struct FlowFlucGfwPp { } PROCESS_SWITCH(FlowFlucGfwPp, processData, "Process analysis for non-derived data", false); - void processq2(soa::Filtered>::iterator const& collision, aod::BCsWithTimestamps const&, GFWTracks const& tracks) + void processq2(soa::Filtered>::iterator const& collision, aod::BCsWithTimestamps const&, GFWTracks const& tracks) { float count{0.5}; registry.fill(HIST("mq2/eventcounter"), count++); @@ -1323,14 +1436,15 @@ struct FlowFlucGfwPp { } registry.fill(HIST("mq2/eventcounter"), count++); - const XAxis xaxis{getCentrality(collision), tracks.size(), -1.0}; + const XAxis xaxis{getCentrality(collision), collision.multNTracksPV(), -1.0}; if (!eventSelected(collision, xaxis.multiplicity, xaxis.centrality, -1)) return; const auto centr = xaxis.centrality; const auto multi = xaxis.multiplicity; - const auto qvecPos = computeqnVec(collision, false); - const auto qvecNeg = computeqnVec(collision, true); + const auto qvecTPC = calcQVec(2, tracks, collision); + const auto qvecPos = computeqnVec(qvecTPC, false); + const auto qvecNeg = computeqnVec(qvecTPC, true); if (!std::isfinite(qvecPos) || !std::isfinite(qvecNeg) || qvecPos < 0 || qvecNeg < 0) { return;