From 26fc012f458db037814f43ae65d0e76c820e7cf6 Mon Sep 17 00:00:00 2001 From: wenyaCern Date: Mon, 22 Jun 2026 16:28:52 +0200 Subject: [PATCH 1/4] Add cutoff of eta for q-vec and MC process --- PWGCF/Flow/Tasks/flowFlucGfwPp.cxx | 177 ++- PWGCF/Flow/flowFlucGfwPp.cxx | 1650 ++++++++++++++++++++++++++++ 2 files changed, 1794 insertions(+), 33 deletions(-) create mode 100644 PWGCF/Flow/flowFlucGfwPp.cxx diff --git a/PWGCF/Flow/Tasks/flowFlucGfwPp.cxx b/PWGCF/Flow/Tasks/flowFlucGfwPp.cxx index b481ac7068b..b1f3d561cfc 100644 --- a/PWGCF/Flow/Tasks/flowFlucGfwPp.cxx +++ b/PWGCF/Flow/Tasks/flowFlucGfwPp.cxx @@ -121,7 +121,7 @@ struct FlowFlucGfwPp { O2_DEFINE_CONFIGURABLE(cfgNbootstrap, int, 10, "Number of subsamples") 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(cfgCentEstimator, int, 0, "0:FT0C; 1:FT0CVariant1; 2:FT0M; 3:FV0A, 4:NTPV, 5:NGlobals") 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") @@ -171,6 +171,7 @@ struct FlowFlucGfwPp { O2_DEFINE_CONFIGURABLE(cfgUseNegativeEtaHalfForq2, bool, true, "If true, use -eta half for qn selection; otherwise use +eta half"); O2_DEFINE_CONFIGURABLE(cfgQnSelectionHarmonic, int, 2, "Harmonic n used to build the reduced q_n vector for event shape selection, use 2 for q2 and 3 for q3"); O2_DEFINE_CONFIGURABLE(cfgQnHistMax, float, 6., "Upper range for q_n calibration histograms"); + O2_DEFINE_CONFIGURABLE(cfgQnTrkAbsEtaMax, float, 0.5, "Upper range for abs eta of tracks contributing to q_n"); O2_DEFINE_CONFIGURABLE(cfgBypassQnSelection, bool, false, "Bypass q_n event shape selection and fill one integral q-bin"); 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"); @@ -277,8 +278,7 @@ struct FlowFlucGfwPp { kCentFT0M, kCentFV0A, kCentNTPV, - kCentNGlobal, - kCentMFT + kCentNGlobal }; enum EventSelFlags { kFilteredEvent = 1, @@ -337,11 +337,11 @@ struct FlowFlucGfwPp { o2::framework::expressions::Filter collisionFilter = nabs(aod::collision::posZ) < cfgVtxZ; o2::framework::expressions::Filter trackFilter = nabs(aod::track::eta) < cfgEta && aod::track::pt > cfgPtmin&& aod::track::pt < cfgPtmax && (aod::track::itsChi2NCl < cfgChi2PrITSCls) && (aod::track::tpcChi2NCl < cfgChi2PrTPCCls) && nabs(aod::track::dcaZ) < cfgDCAz; - Preslice perCollision = aod::track::collisionId; 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 GFWTracksMC = soa::Filtered>; void init(InitContext const&) { @@ -415,8 +415,7 @@ struct FlowFlucGfwPp { {kCentFT0M, "FT0M"}, {kCentFV0A, "FV0A"}, {kCentNTPV, "NTPV"}, - {kCentNGlobal, "NGlobals"}, - {kCentMFT, "MFT"}}; + {kCentNGlobal, "NGlobals"}}; sCentralityEstimator = centEstimatorMap.at(cfgCentEstimator); sCentralityEstimator += " centrality (%)"; AxisSpec centAxis = {o2::analysis::gfwflowflucpp::centbinning, sCentralityEstimator.c_str()}; @@ -503,7 +502,7 @@ struct FlowFlucGfwPp { } } - if (doprocessData) { + if (doprocessData || doprocessMC) { registry.add("trackQA/before/phi_eta_vtxZ", "", {HistType::kTH3D, {phiAxis, etaAxis, vtxAxis}}); registry.add("trackQA/before/pt_dcaXY_dcaZ", "", {HistType::kTH3D, {ptAxis, dcaXYAXis, dcaZAXis}}); registry.add("trackQA/before/nch_pt", "#it{p}_{T} vs multiplicity; N_{ch}; #it{p}_{T}", {HistType::kTH2D, {nchAxis, ptAxis}}); @@ -530,24 +529,20 @@ struct FlowFlucGfwPp { registry.add("eventQA/before/globalTracks_multV0A", "", {HistType::kTH2D, {v0aAxis, nchAxis}}); registry.add("eventQA/before/multV0A_multT0A", "", {HistType::kTH2D, {t0aAxis, v0aAxis}}); - if (doprocessData) { - registry.add("eventQA/before/centrality", "", {HistType::kTH1D, {centAxis}}); - registry.add("eventQA/before/globalTracks_centT0C", "", {HistType::kTH2D, {centAxis, nchAxis}}); - registry.add("eventQA/before/PVTracks_centT0C", "", {HistType::kTH2D, {centAxis, multpvAxis}}); - registry.add("eventQA/before/multT0C_centT0C", "", {HistType::kTH2D, {centAxis, t0cAxis}}); - - registry.add("eventQA/before/centT0M_centT0C", "", {HistType::kTH2D, {centAxis, centAxis}}); - registry.add("eventQA/before/centV0A_centT0C", "", {HistType::kTH2D, {centAxis, centAxis}}); - registry.add("eventQA/before/centGlobal_centT0C", "", {HistType::kTH2D, {centAxis, centAxis}}); - registry.add("eventQA/before/centNTPV_centT0C", "", {HistType::kTH2D, {centAxis, centAxis}}); - registry.add("eventQA/before/centMFT_centT0C", "", {HistType::kTH2D, {centAxis, centAxis}}); - - if (cfgIsMC) { - registry.add("MCGen/trackQA/phi_eta_vtxZ", "", {HistType::kTH3D, {phiAxis, etaAxis, vtxAxis}}); - registry.add("MCGen/trackQA/nch_pt", "#it{p}_{T} vs multiplicity; N_{ch}; #it{p}_{T}", {HistType::kTH2D, {nchAxis, ptAxis}}); - registry.add("MCGen/trackQA/pt_ref", "", {HistType::kTH1D, {{100, o2::analysis::gfwflowflucpp::ptreflow, o2::analysis::gfwflowflucpp::ptrefup}}}); - registry.add("MCGen/trackQA/pt_poi", "", {HistType::kTH1D, {{100, o2::analysis::gfwflowflucpp::ptpoilow, o2::analysis::gfwflowflucpp::ptpoiup}}}); - } + registry.add("eventQA/before/centrality", "", {HistType::kTH1D, {centAxis}}); + registry.add("eventQA/before/globalTracks_centT0C", "", {HistType::kTH2D, {centAxis, nchAxis}}); + registry.add("eventQA/before/PVTracks_centT0C", "", {HistType::kTH2D, {centAxis, multpvAxis}}); + registry.add("eventQA/before/multT0C_centT0C", "", {HistType::kTH2D, {centAxis, t0cAxis}}); + + registry.add("eventQA/before/centT0M_centT0C", "", {HistType::kTH2D, {centAxis, centAxis}}); + registry.add("eventQA/before/centV0A_centT0C", "", {HistType::kTH2D, {centAxis, centAxis}}); + registry.add("eventQA/before/centGlobal_centT0C", "", {HistType::kTH2D, {centAxis, centAxis}}); + registry.add("eventQA/before/centNTPV_centT0C", "", {HistType::kTH2D, {centAxis, centAxis}}); + if (cfgIsMC || doprocessMC) { + registry.add("MCGen/trackQA/phi_eta_vtxZ", "", {HistType::kTH3D, {phiAxis, etaAxis, vtxAxis}}); + registry.add("MCGen/trackQA/nch_pt", "#it{p}_{T} vs multiplicity; N_{ch}; #it{p}_{T}", {HistType::kTH2D, {nchAxis, ptAxis}}); + registry.add("MCGen/trackQA/pt_ref", "", {HistType::kTH1D, {{100, o2::analysis::gfwflowflucpp::ptreflow, o2::analysis::gfwflowflucpp::ptrefup}}}); + registry.add("MCGen/trackQA/pt_poi", "", {HistType::kTH1D, {{100, o2::analysis::gfwflowflucpp::ptpoilow, o2::analysis::gfwflowflucpp::ptpoiup}}}); } registry.addClone("eventQA/before/", "eventQA/after/"); @@ -898,6 +893,7 @@ struct FlowFlucGfwPp { std::string name = Form("%s_%d_%s", shapeSel.c_str(), jese, it->Head.c_str()); std::string title = it->Head + std::string("_ese"); oba->Add(new TNamed(name.c_str(), title.c_str())); + } } } @@ -1020,14 +1016,14 @@ struct FlowFlucGfwPp { registry.fill(HIST("qvecQA/ChTracks"), trk.pt(), trk.eta(), trk.phi()); } - if (trk.eta() > 0) { + if (trk.eta() > 0 && fabs(trk.eta())< cfgQnTrkAbsEtaMax) { // 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) { + } else if (trk.eta() < 0 && fabs(trk.eta())< cfgQnTrkAbsEtaMax) { // 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); @@ -1172,6 +1168,48 @@ struct FlowFlucGfwPp { lRandom, qPtmp, run); } + template + void processGenCollision(TCollision collision, TParticles particles, const int& mcCollisionId, const XAxis& xaxis, const int& run, const int& qPtmp) + { + if (xaxis.multiplicity < cfgFixedMultMin || xaxis.multiplicity > cfgFixedMultMax) + return; + + if (cfgFillQA && xaxis.centrality >= 0) + registry.fill(HIST("eventQA/after/centrality"), xaxis.centrality); + if (cfgFillQA) + registry.fill(HIST("eventQA/after/multiplicity"), xaxis.multiplicity); + + fGFW->Clear(); + float lRandom = fRndm->Rndm(); + float vtxz = collision.posZ(); + + AcceptedTracks acceptedTracks{0, 0, 0, 0}; + for (const auto& particle : particles) { + if (particle.mcCollisionId() != mcCollisionId) + continue; + processTrack(particle, vtxz, xaxis.multiplicity, run, acceptedTracks); + } + + if (cfgConsistentEventFlag & kRequireBothEtaSides) + if (!acceptedTracks.nPos || !acceptedTracks.nNeg) + return; + if (cfgConsistentEventFlag & kRequireFullFourParticleTracks) + if (acceptedTracks.nFull < kMinTracksForFourParticleCorrelation) + return; + if (cfgConsistentEventFlag & kRequireTwoTracksInBothEtaSides) + if (acceptedTracks.nPos < kMinTracksPerEtaSideForGapCorrelation || + acceptedTracks.nNeg < kMinTracksPerEtaSideForGapCorrelation) + return; + if (cfgConsistentEventFlag & kRequireTwoTracksInThreeEtaRegions) + if (acceptedTracks.nPos < kMinTracksPerEtaRegionForThreeSubevents || + acceptedTracks.nMid < kMinTracksPerEtaRegionForThreeSubevents || + acceptedTracks.nNeg < kMinTracksPerEtaRegionForThreeSubevents) + return; + + fillOutputContainers(cfgUseNch ? static_cast(xaxis.multiplicity) : xaxis.centrality, + lRandom, qPtmp, run); + } + bool isStable(int pdg) { if (std::abs(pdg) == PDG_t::kPiPlus) @@ -1334,8 +1372,6 @@ struct FlowFlucGfwPp { return collision.centNTPV(); case kCentNGlobal: return collision.centNGlobal(); - case kCentMFT: - return collision.centMFT(); default: return collision.centFT0C(); } @@ -1352,7 +1388,6 @@ struct FlowFlucGfwPp { registry.fill(HIST("eventQA/") + HIST(FillTimeName[ft]) + HIST("centV0A_centT0C"), collision.centFT0C(), collision.centFV0A()); registry.fill(HIST("eventQA/") + HIST(FillTimeName[ft]) + HIST("centGlobal_centT0C"), collision.centFT0C(), collision.centNGlobal()); registry.fill(HIST("eventQA/") + HIST(FillTimeName[ft]) + HIST("centNTPV_centT0C"), collision.centFT0C(), collision.centNTPV()); - registry.fill(HIST("eventQA/") + HIST(FillTimeName[ft]) + HIST("centMFT_centT0C"), collision.centFT0C(), collision.centMFT()); } registry.fill(HIST("eventQA/") + HIST(FillTimeName[ft]) + HIST("globalTracks_PVTracks"), collision.multNTracksPV(), xaxis.multiplicity); registry.fill(HIST("eventQA/") + HIST(FillTimeName[ft]) + HIST("globalTracks_multT0A"), collision.multFT0A(), xaxis.multiplicity); @@ -1406,8 +1441,7 @@ struct FlowFlucGfwPp { void processData(soa::Filtered>::iterator const& collision, + aod::CentFV0As, aod::CentNTPVs, aod::CentNGlobals>>::iterator const& collision, aod::BCsWithTimestamps const&, GFWTracks const& tracks) { auto bc = collision.bc_as(); @@ -1487,7 +1521,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}; fillQnEventCounter(count++); @@ -1521,6 +1555,83 @@ struct FlowFlucGfwPp { fillQnCalibrationHistograms(centr, multi, qvecPos, qvecNeg); } PROCESS_SWITCH(FlowFlucGfwPp, processq2, "Process analysis for filling q_n-vector calibration histograms", true); + + void processMC(soa::Filtered>::iterator const& collision, + aod::BCsWithTimestamps const&, GFWTracksMC const& tracks, aod::McCollisions const&, aod::McParticles const& mcParticles) + { + auto bc = collision.bc_as(); + int run = bc.runNumber(); + if (run != lastRun) { + lastRun = run; + LOGF(info, "run = %d", run); + if (cfgRunByRun) { + if (std::find(runNumbers.begin(), runNumbers.end(), run) == runNumbers.end()) { + LOGF(info, "Creating histograms for run %d", run); + createRunByRunHistograms(run); + runNumbers.push_back(run); + } + if (!cfgFillWeights) + loadCorrections(bc); + } + } + if (!cfgFillWeights && !cfgRunByRun) + loadCorrections(bc); + + registry.fill(HIST("eventQA/eventSel"), kFilteredEvent); + if (cfgRunByRun) + th1sList[run][hEventSel]->Fill(kFilteredEvent); + + if (!collision.sel8()) + return; + + registry.fill(HIST("eventQA/eventSel"), kSel8); + if (cfgRunByRun) + th1sList[run][hEventSel]->Fill(kSel8); + + if (cfgDoOccupancySel) { + int occupancy = collision.trackOccupancyInTimeRange(); + if (occupancy < 0 || occupancy > cfgOccupancySelection) + return; + } + + registry.fill(HIST("eventQA/eventSel"), kOccupancy); + if (cfgRunByRun) + th1sList[run][hEventSel]->Fill(kOccupancy); + + const XAxis xaxis{ + getCentrality(collision), + collision.multNTracksPV(), + (cfgTimeDependent) ? getTimeSinceStartOfFill(bc.timestamp(), *firstRunOfCurrentFill) : -1.0}; + + if (cfgTimeDependent && run == *firstRunOfCurrentFill && + firstRunOfCurrentFill != o2::analysis::gfwflowflucpp::firstRunsOfFill.end() - 1) + ++firstRunOfCurrentFill; + + if (cfgFillQA) + fillEventQA(collision, xaxis); + + registry.fill(HIST("eventQA/before/centrality"), xaxis.centrality); + registry.fill(HIST("eventQA/before/multiplicity"), xaxis.multiplicity); + + if (!eventSelected(collision, xaxis.multiplicity, xaxis.centrality, run)) + return; + + if (cfgFillQA) + fillEventQA(collision, xaxis); + + processCollision(collision, tracks, xaxis, run, 0); + + if (!collision.has_mcCollision()) + return; + + const auto genCollision = collision.template mcCollision_as(); + processGenCollision(genCollision, mcParticles, collision.mcCollisionId(), xaxis, run, 0); + } + PROCESS_SWITCH(FlowFlucGfwPp, processMC, "Process analysis for Monte-Carlo data", false); + }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) diff --git a/PWGCF/Flow/flowFlucGfwPp.cxx b/PWGCF/Flow/flowFlucGfwPp.cxx new file mode 100644 index 00000000000..de048182729 --- /dev/null +++ b/PWGCF/Flow/flowFlucGfwPp.cxx @@ -0,0 +1,1650 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file flowFlucGfwPp.cxx +/// \brief GFW task for Event Shape Engineering studies in pp collisions +/// \author Wenya Wu, TUM, wenya.wu@cern.ch, David Krylenkov, TUM + +#include "PWGCF/GenericFramework/Core/FlowContainer.h" +#include "PWGCF/GenericFramework/Core/GFW.h" +#include "PWGCF/GenericFramework/Core/GFWConfig.h" +#include "PWGCF/GenericFramework/Core/GFWWeights.h" + +#include "Common/CCDB/EventSelectionParams.h" +#include "Common/CCDB/TriggerAliases.h" +#include "Common/DataModel/Centrality.h" +#include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/Multiplicity.h" +#include "Common/DataModel/TrackSelectionTables.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace o2; +using namespace o2::framework; +using namespace o2::analysis::genericframework; + +#define O2_DEFINE_CONFIGURABLE(NAME, TYPE, DEFAULT, HELP) Configurable NAME{#NAME, DEFAULT, HELP}; + +namespace o2::analysis::gfwflowflucpp +{ +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}; +float ptpoilow = 0.2, ptpoiup = 10.0; +float ptreflow = 0.2, ptrefup = 3.0; +float ptlow = 0.2, ptup = 10.0; +int etabins = 16; +float etalow = -0.8, etaup = 0.8; +int vtxZbins = 40; +float vtxZlow = -10.0, vtxZup = 10.0; +int phibins = 72; +float philow = 0.0; +float phiup = o2::constants::math::TwoPI; +int nchbins = 300; +float nchlow = 0; +float nchup = 3000; +std::vector centbinning(90); +int nBootstrap = 10; +GFWRegions regions; +GFWCorrConfigs configs; +std::vector multGlobalCorrCutPars; +std::vector multPVCorrCutPars; +std::vector multGlobalPVCorrCutPars; +std::vector multGlobalV0ACutPars; +std::vector multGlobalT0ACutPars; +std::vector firstRunsOfFill; +} // namespace o2::analysis::gfwflowflucpp + +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; + static constexpr int kRequireTwoTracksInBothEtaSides = 4; + static constexpr int kRequireTwoTracksInThreeEtaRegions = 8; + + static constexpr int kMinTracksForFourParticleCorrelation = 4; + static constexpr int kMinTracksPerEtaSideForGapCorrelation = 2; + static constexpr int kMinTracksPerEtaRegionForThreeSubevents = 2; + + static constexpr int EllipticQVectorHarmonic = 2; + static constexpr int TriangularQVectorHarmonic = 3; + + O2_DEFINE_CONFIGURABLE(cfgNbootstrap, int, 10, "Number of subsamples") + 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)") + O2_DEFINE_CONFIGURABLE(cfgTimeDependent, bool, false, "Fill output as function of time (for contamination studies)") + O2_DEFINE_CONFIGURABLE(cfgFirstRunsOfFill, std::vector, {}, "First runs of a fill for time dependent analysis") + O2_DEFINE_CONFIGURABLE(cfgFillQA, bool, false, "Fill QA histograms") + O2_DEFINE_CONFIGURABLE(cfgUseAdditionalEventCut, bool, false, "Use additional event cut on mult correlations") + O2_DEFINE_CONFIGURABLE(cfgEfficiency, std::string, "", "CCDB path to efficiency object") + O2_DEFINE_CONFIGURABLE(cfgAcceptance, std::string, "", "CCDB path to acceptance object") + O2_DEFINE_CONFIGURABLE(cfgDCAxyNSigma, float, 7, "Cut on number of sigma deviations from expected DCA in the transverse direction"); + O2_DEFINE_CONFIGURABLE(cfgDCAxy, std::string, "(0.0026+0.005/(x^1.01))", "Functional form of pt-dependent DCAxy cut"); + O2_DEFINE_CONFIGURABLE(cfgDCAz, float, 2, "Cut on DCA in the longitudinal direction (cm)"); + O2_DEFINE_CONFIGURABLE(cfgNTPCCls, float, 50, "Cut on number of TPC clusters found"); + O2_DEFINE_CONFIGURABLE(cfgNTPCXrows, float, 70, "Cut on number of TPC crossed rows"); + O2_DEFINE_CONFIGURABLE(cfgMinNITSCls, float, 5, "Cut on minimum number of ITS clusters found"); + O2_DEFINE_CONFIGURABLE(cfgChi2PrITSCls, float, 36, "Cut on chi^2 per ITS clusters found"); + O2_DEFINE_CONFIGURABLE(cfgChi2PrTPCCls, float, 2.5, "Cut on chi^2 per TPC clusters found"); + O2_DEFINE_CONFIGURABLE(cfgPtmin, float, 0.2, "minimum pt (GeV/c)"); + O2_DEFINE_CONFIGURABLE(cfgPtmax, float, 10, "maximum pt (GeV/c)"); + O2_DEFINE_CONFIGURABLE(cfgEta, float, 0.8, "eta cut"); + O2_DEFINE_CONFIGURABLE(cfgVtxZ, float, 10, "vertex cut (cm)"); + O2_DEFINE_CONFIGURABLE(cfgNoSameBunchPileupCut, bool, false, "kNoSameBunchPileupCut"); + O2_DEFINE_CONFIGURABLE(cfgIsGoodZvtxFT0vsPV, bool, false, "kIsGoodZvtxFT0vsPV"); + O2_DEFINE_CONFIGURABLE(cfgIsGoodITSLayersAll, bool, false, "kIsGoodITSLayersAll"); + O2_DEFINE_CONFIGURABLE(cfgNoCollInTimeRangeStandard, bool, false, "kNoCollInTimeRangeStandard"); + O2_DEFINE_CONFIGURABLE(cfgDoOccupancySel, bool, false, "Bool for event selection on detector occupancy"); + O2_DEFINE_CONFIGURABLE(cfgOccupancySelection, int, 5000, "Max occupancy selection, -999 to disable"); + O2_DEFINE_CONFIGURABLE(cfgMultCut, bool, false, "Use additional event cut on mult correlations"); + O2_DEFINE_CONFIGURABLE(cfgTVXinTRD, bool, false, "Use kTVXinTRD (reject TRD triggered events)"); + O2_DEFINE_CONFIGURABLE(cfgIsVertexITSTPC, bool, true, "Selects collisions with at least one ITS-TPC track"); + O2_DEFINE_CONFIGURABLE(cfgMagField, float, 99999, "Configurable magnetic field; default CCDB will be queried"); + O2_DEFINE_CONFIGURABLE(cfgFixedMultMin, int, 1, "Minimum for fixed nch range"); + O2_DEFINE_CONFIGURABLE(cfgUseMultiplicityFlowWeights, bool, true, "Enable or disable the use of multiplicity-based event weighting"); + O2_DEFINE_CONFIGURABLE(cfgFixedMultMax, int, 3000, "Maximum for fixed nch range"); + O2_DEFINE_CONFIGURABLE(cfgConsistentEventFlag, int, 0, "Flag to select consistent events - 0: off, 1: v2{2} gap calculable, 2: v2{4} full calculable, 4: v2{4} gap calculable, 8: v2{4} 3sub calculable"); + Configurable> cfgMultGlobalCutPars{"cfgMultGlobalCutPars", std::vector{2272.16, -76.6932, 1.01204, -0.00631545, 1.59868e-05, 136.336, -4.97006, 0.121199, -0.0015921, 7.66197e-06}, "Global vs FT0C multiplicity cut parameter values"}; + Configurable> cfgMultPVCutPars{"cfgMultPVCutPars", std::vector{3074.43, -106.192, 1.46176, -0.00968364, 2.61923e-05, 182.128, -7.43492, 0.193901, -0.00256715, 1.22594e-05}, "PV vs FT0C multiplicity cut parameter values"}; + Configurable> cfgMultGlobalPVCutPars{"cfgMultGlobalPVCutPars", std::vector{-0.223013, 0.715849, 0.664242, 0.0829653, -0.000503733, 1.21185e-06}, "Global vs PV multiplicity cut parameter values"}; + O2_DEFINE_CONFIGURABLE(cfgMultCorrHighCutFunction, std::string, "[0] + [1]*x + [2]*x*x + [3]*x*x*x + [4]*x*x*x*x + 3.*([5] + [6]*x + [7]*x*x + [8]*x*x*x + [9]*x*x*x*x)", "Functional for multiplicity correlation cut"); + O2_DEFINE_CONFIGURABLE(cfgMultCorrLowCutFunction, std::string, "[0] + [1]*x + [2]*x*x + [3]*x*x*x + [4]*x*x*x*x - 3.*([5] + [6]*x + [7]*x*x + [8]*x*x*x + [9]*x*x*x*x)", "Functional for multiplicity correlation cut"); + O2_DEFINE_CONFIGURABLE(cfgMultGlobalPVCorrCutFunction, std::string, "[0] + [1]*x + 3*([2] + [3]*x + [4]*x*x + [5]*x*x*x)", "Functional for global vs pv multiplicity correlation cut"); + O2_DEFINE_CONFIGURABLE(cfgNumQnBins, int, 10, "Number of qn bins"); + 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 qn selection; otherwise use +eta half"); + O2_DEFINE_CONFIGURABLE(cfgQnSelectionHarmonic, int, 2, "Harmonic n used to build the reduced q_n vector for event shape selection, use 2 for q2 and 3 for q3"); + O2_DEFINE_CONFIGURABLE(cfgQnHistMax, float, 6., "Upper range for q_n calibration histograms"); + O2_DEFINE_CONFIGURABLE(cfgQnTrkAbsEtaMax, float, 0.5, "Upper range for abs eta of tracks contributing to q_n"); + O2_DEFINE_CONFIGURABLE(cfgBypassQnSelection, bool, false, "Bypass q_n event shape selection and fill one integral q-bin"); + 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 { + O2_DEFINE_CONFIGURABLE(cfgMultGlobalASideCorrCutFunction, std::string, "[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x + [10]*([5]+[6]*x+[7]*x*x+[8]*x*x*x+[9]*x*x*x*x)", "Functional for global vs V0A multiplicity low correlation cut"); + Configurable> cfgMultGlobalV0ACutPars{"cfgMultGlobalV0ACutPars", std::vector{567.785, 172.715, 0.77888, -0.00693466, 1.40564e-05, 679.853, 66.8068, -0.444332, 0.00115002, -4.92064e-07}, "Global vs FV0A multiplicity cut parameter values"}; + Configurable> cfgMultGlobalT0ACutPars{"cfgMultGlobalT0ACutPars", std::vector{241.618, 61.8402, 0.348049, -0.00306078, 6.20357e-06, 315.235, 29.1491, -0.188639, 0.00044528, -9.08912e-08}, "Global vs FT0A multiplicity cut parameter values"}; + O2_DEFINE_CONFIGURABLE(cfgGlobalV0ALowSigma, float, -3, "Number of sigma deviations below expected value in global vs V0A correlation"); + O2_DEFINE_CONFIGURABLE(cfgGlobalV0AHighSigma, float, 4, "Number of sigma deviations above expected value in global vs V0A correlation"); + O2_DEFINE_CONFIGURABLE(cfgGlobalT0ALowSigma, float, -3., "Number of sigma deviations below expected value in global vs T0A correlation"); + O2_DEFINE_CONFIGURABLE(cfgGlobalT0AHighSigma, float, 4, "Number of sigma deviations above expected value in global vs T0A correlation"); + } cfgGlobalAsideCorrCuts; + + Configurable cfgGfwBinning{"cfgGfwBinning", + {40, 16, 72, 300, 0, 3000, 0.2, 10.0, 0.2, 3.0, {0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1, 1.25, 1.5, 1.75, 2, 2.25, 2.5, 2.75, 3, 3.25, 3.5, 3.75, 4, 4.5, 5, 5.5, 6, 7, 8, 9, 10}, {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90}}, + "Configuration for binning"}; + + Configurable cfgRegions{"cfgRegions", + { + {"refN", "refP", "refFull"}, + {-0.8, 0.4, -0.8}, + {-0.4, 0.8, 0.8}, + {0, 0, 0}, // pT bins + {1, 1, 1} // bitmask + }, + "Configurations for GFW regions"}; + + Configurable cfgCorrConfig{"cfgCorrConfig", + {{"refN {2 -2}", + "refN {2 2 -2 -2}", + "refN {2 2 2 -2 -2 -2}", + "refN {2 2 2 2 -2 -2 -2 -2}", + "refP {2 -2}", + "refP {2 2 -2 -2}", + "refP {2 2 2 -2 -2 -2}", + "refP {2 2 2 2 -2 -2 -2 -2}", + "refN {2} refP {-2}", + "refN {2 2} refP {-2 -2}", + "refFull {2 -2}", + "refFull {2 2 -2 -2}", + "refFull {2 2 2 -2 -2 -2}", + "refFull {2 2 2 2 -2 -2 -2 -2}"}, + {"ChNeg22", + "ChNeg24", + "ChNeg26", + "ChNeg28", + "ChPos22", + "ChPos24", + "ChPos26", + "ChPos28", + "ChGap22", + "ChGap24", + "ChFull22", + "ChFull24", + "ChFull26", + "ChFull28"}, + {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}}, + "Configurations for pp ESE v2 cumulants"}; + + // Connect to ccdb + Service ccdb; + + struct Config { + TH1D* mEfficiency = nullptr; + GFWWeights* mAcceptance; + bool correctionsLoaded = false; + } cfg; + + // Define output + OutputObj fFC{FlowContainer("FlowContainer")}; + OutputObj fFCgen{FlowContainer("FlowContainer_gen")}; + HistogramRegistry registry{"registry"}; + + // QA outputs + std::map>> th1sList; + std::map>> tpfsList; + std::map>> th3sList; + enum OutputTH1Names { + hPhi = 0, + hEta, + hVtxZ, + hMult, + hCent, + hEventSel, + kCount_TH1Names + }; + enum OutputTProfileNames { + pfCorr22 = 0, + kCount_TProfileNames + }; + // NUA outputs + enum OutputTH3Names { + hNUAref = 0, + kCount_TH3Names + }; + enum CentEstimators { + kCentFT0C = 0, + kCentFT0CVariant1, + kCentFT0M, + kCentFV0A, + kCentNTPV, + kCentNGlobal, + kCentMFT + }; + enum EventSelFlags { + kFilteredEvent = 1, + kSel8, + kOccupancy, + kTVXTRD, + kNoSamebunchPU, + kZVtxFT0PV, + kNoCollTRStd, + kVtxITSTPC, + kGoodITSLayers, + kMultCuts, + kTrackCent + }; + + // Define global variables + // Generic Framework + GFW* fGFW = new GFW(); + std::vector corrconfigs; + + TRandom3* fRndm = new TRandom3(0); + TAxis* fSecondAxis; + int lastRun = -1; + std::vector::iterator firstRunOfCurrentFill; + std::vector runNumbers; + + std::string getShapeSel() const + { + return "ese"; + } + + int getNQnOutputBins() + { + return static_cast(cfgBypassQnSelection) ? 1 : static_cast(cfgNumQnBins); + } + + // region indices for consistency flag + int posRegionIndex = -1; + int negRegionIndex = -1; + int fullRegionIndex = -1; + int midRegionIndex = -1; + + // Event selection cuts - Alex + TF1* fMultPVCutLow = nullptr; + TF1* fMultPVCutHigh = nullptr; + TF1* fMultCutLow = nullptr; + TF1* fMultCutHigh = nullptr; + TF1* fMultPVGlobalCutHigh = nullptr; + TF1* fMultGlobalV0ACutLow = nullptr; + TF1* fMultGlobalV0ACutHigh = nullptr; + TF1* fMultGlobalT0ACutLow = nullptr; + TF1* fMultGlobalT0ACutHigh = nullptr; + + TF1* fPtDepDCAxy = nullptr; + + o2::framework::expressions::Filter collisionFilter = nabs(aod::collision::posZ) < cfgVtxZ; + o2::framework::expressions::Filter trackFilter = nabs(aod::track::eta) < cfgEta && aod::track::pt > cfgPtmin&& aod::track::pt < cfgPtmax && (aod::track::itsChi2NCl < cfgChi2PrITSCls) && (aod::track::tpcChi2NCl < cfgChi2PrTPCCls) && nabs(aod::track::dcaZ) < cfgDCAz; + + Preslice perCollision = aod::track::collisionId; + 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 GFWTracksMC = 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()); + o2::analysis::gfwflowflucpp::regions.SetpTDifs(cfgRegions->GetpTDifs()); + o2::analysis::gfwflowflucpp::regions.SetBitmasks(cfgRegions->GetBitmasks()); + o2::analysis::gfwflowflucpp::configs.SetCorrs(cfgCorrConfig->GetCorrs()); + o2::analysis::gfwflowflucpp::configs.SetHeads(cfgCorrConfig->GetHeads()); + o2::analysis::gfwflowflucpp::configs.SetpTDifs(cfgCorrConfig->GetpTDifs()); + o2::analysis::gfwflowflucpp::configs.SetpTCorrMasks(cfgCorrConfig->GetpTCorrMasks()); + + o2::analysis::gfwflowflucpp::regions.Print(); + o2::analysis::gfwflowflucpp::configs.Print(); + + o2::analysis::gfwflowflucpp::ptbinning = cfgGfwBinning->GetPtBinning(); + // o2::analysis::gfwflowflucpp::ptpoilow = cfgGfwBinning->GetPtPOImin(); + // o2::analysis::gfwflowflucpp::ptpoiup = cfgGfwBinning->GetPtPOImax(); + o2::analysis::gfwflowflucpp::ptreflow = cfgGfwBinning->GetPtRefMin(); + o2::analysis::gfwflowflucpp::ptrefup = cfgGfwBinning->GetPtRefMax(); + o2::analysis::gfwflowflucpp::ptlow = cfgPtmin; + o2::analysis::gfwflowflucpp::ptup = cfgPtmax; + o2::analysis::gfwflowflucpp::etabins = cfgGfwBinning->GetEtaBins(); + o2::analysis::gfwflowflucpp::vtxZbins = cfgGfwBinning->GetVtxZbins(); + o2::analysis::gfwflowflucpp::phibins = cfgGfwBinning->GetPhiBins(); + o2::analysis::gfwflowflucpp::philow = 0.0f; + o2::analysis::gfwflowflucpp::phiup = o2::constants::math::TwoPI; + o2::analysis::gfwflowflucpp::nchbins = cfgGfwBinning->GetNchBins(); + o2::analysis::gfwflowflucpp::nchlow = cfgGfwBinning->GetNchMin(); + o2::analysis::gfwflowflucpp::nchup = cfgGfwBinning->GetNchMax(); + o2::analysis::gfwflowflucpp::centbinning = cfgGfwBinning->GetCentBinning(); + cfgGfwBinning->Print(); + + o2::analysis::gfwflowflucpp::multGlobalCorrCutPars = cfgMultGlobalCutPars; + o2::analysis::gfwflowflucpp::multPVCorrCutPars = cfgMultPVCutPars; + o2::analysis::gfwflowflucpp::multGlobalPVCorrCutPars = cfgMultGlobalPVCutPars; + o2::analysis::gfwflowflucpp::multGlobalV0ACutPars = cfgGlobalAsideCorrCuts.cfgMultGlobalV0ACutPars; + o2::analysis::gfwflowflucpp::multGlobalT0ACutPars = cfgGlobalAsideCorrCuts.cfgMultGlobalT0ACutPars; + o2::analysis::gfwflowflucpp::firstRunsOfFill = cfgFirstRunsOfFill; + if (cfgTimeDependent && !std::is_sorted(o2::analysis::gfwflowflucpp::firstRunsOfFill.begin(), o2::analysis::gfwflowflucpp::firstRunsOfFill.end())) { + std::sort(o2::analysis::gfwflowflucpp::firstRunsOfFill.begin(), o2::analysis::gfwflowflucpp::firstRunsOfFill.end()); + } + firstRunOfCurrentFill = o2::analysis::gfwflowflucpp::firstRunsOfFill.begin(); + + AxisSpec phiAxis = {o2::analysis::gfwflowflucpp::phibins, o2::analysis::gfwflowflucpp::philow, o2::analysis::gfwflowflucpp::phiup, "#phi"}; + AxisSpec etaAxis = {o2::analysis::gfwflowflucpp::etabins, -cfgEta, cfgEta, "#eta"}; + AxisSpec vtxAxis = {o2::analysis::gfwflowflucpp::vtxZbins, -cfgVtxZ, cfgVtxZ, "Vtx_{z} (cm)"}; + AxisSpec ptAxis = {o2::analysis::gfwflowflucpp::ptbinning, "#it{p}_{T} GeV/#it{c}"}; + + std::string sCentralityEstimator; + std::map centEstimatorMap = { + {kCentFT0C, "FT0C"}, + {kCentFT0CVariant1, "FT0C variant 1"}, + {kCentFT0M, "FT0M"}, + {kCentFV0A, "FV0A"}, + {kCentNTPV, "NTPV"}, + {kCentNGlobal, "NGlobals"}, + {kCentMFT, "MFT"}}; + sCentralityEstimator = centEstimatorMap.at(cfgCentEstimator); + sCentralityEstimator += " centrality (%)"; + AxisSpec centAxis = {o2::analysis::gfwflowflucpp::centbinning, sCentralityEstimator.c_str()}; + std::vector nchbinning; + int nchskip = (o2::analysis::gfwflowflucpp::nchup - o2::analysis::gfwflowflucpp::nchlow) / o2::analysis::gfwflowflucpp::nchbins; + for (int i = 0; i <= o2::analysis::gfwflowflucpp::nchbins; ++i) { + nchbinning.push_back(nchskip * i + o2::analysis::gfwflowflucpp::nchlow + 0.5); + } + AxisSpec nchAxis = {nchbinning, "N_{ch}"}; + std::vector bbinning(201); + std::generate(bbinning.begin(), bbinning.end(), [n = -0.1, step = 0.1]() mutable { + n += step; + return n; + }); + AxisSpec bAxis = {bbinning, "#it{b}"}; + AxisSpec t0cAxis = {1000, 0, 10000, "N_{ch} (T0C)"}; + AxisSpec t0aAxis = {300, 0, 30000, "N_{ch} (T0A)"}; + AxisSpec v0aAxis = {800, 0, 80000, "N_{ch} (V0A)"}; + AxisSpec multpvAxis = {600, 0, 600, "N_{ch} (PV)"}; + AxisSpec dcaZAXis = {200, -2, 2, "DCA_{z} (cm)"}; + AxisSpec dcaXYAXis = {200, -0.5, 0.5, "DCA_{xy} (cm)"}; + std::vector timebinning(289); + std::generate(timebinning.begin(), timebinning.end(), [n = -24 / 288., step = 24 / 288.]() mutable { + n += step; + return n; + }); + AxisSpec timeAxis = {timebinning, "time (hrs)"}; + + AxisSpec multAxis = (cfgTimeDependent) ? timeAxis : (cfgUseNch) ? nchAxis + : centAxis; + + ccdb->setURL("http://alice-ccdb.cern.ch"); + ccdb->setCaching(true); + ccdb->setLocalObjectValidityChecking(); + + int64_t now = std::chrono::duration_cast(std::chrono::system_clock::now().time_since_epoch()).count(); + ccdb->setCreatedNotAfter(now); + + 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 (static_cast(cfgBypassQnSelection)) { + LOGF(info, "Bypassing q_n event shape selection, all accepted events will be filled into the integral bin ese_0"); + if (static_cast(cfgNumQnBins) != 1) { + LOGF(warning, "cfgBypassQnSelection is on, cfgNumQnBins=%d will be ignored and only one output q-bin will be made", static_cast(cfgNumQnBins)); + } + } else { + LOGF(info, "Event-shape selector uses q_%d from the %s eta half", + static_cast(cfgQnSelectionHarmonic), + static_cast(cfgUseNegativeEtaHalfForq2) ? "negative" : "positive"); + } + + 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) { + const int qnHarmonic = static_cast(cfgQnSelectionHarmonic); + const double qnHistMax = static_cast(cfgQnHistMax); + + if (qnHarmonic == EllipticQVectorHarmonic) { + 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, qnHistMax}}); + registry.add("mq2/h2_cent_q2_etaneg", ";Centrality;#it{q}_{2}^{#eta neg};", HistType::kTH2D, {{100, 0, 100}, {600, 0, qnHistMax}}); + registry.add("mq2/h2_mult_q2_etapos", ";Multiplicity;#it{q}_{2}^{#eta pos};", HistType::kTH2D, {{150, 0, 150}, {600, 0, qnHistMax}}); + registry.add("mq2/h2_mult_q2_etaneg", ";Multiplicity;#it{q}_{2}^{#eta neg};", HistType::kTH2D, {{150, 0, 150}, {600, 0, qnHistMax}}); + } + + if (qnHarmonic == TriangularQVectorHarmonic) { + registry.add("mq3/eventcounter", "", HistType::kTH1F, {{10, 0, 10}}); + registry.add("mq3/h2_cent_q3_etapos", ";Centrality;#it{q}_{3}^{#eta pos};", HistType::kTH2D, {{100, 0, 100}, {600, 0, qnHistMax}}); + registry.add("mq3/h2_cent_q3_etaneg", ";Centrality;#it{q}_{3}^{#eta neg};", HistType::kTH2D, {{100, 0, 100}, {600, 0, qnHistMax}}); + registry.add("mq3/h2_mult_q3_etapos", ";Multiplicity;#it{q}_{3}^{#eta pos};", HistType::kTH2D, {{150, 0, 150}, {600, 0, qnHistMax}}); + registry.add("mq3/h2_mult_q3_etaneg", ";Multiplicity;#it{q}_{3}^{#eta neg};", HistType::kTH2D, {{150, 0, 150}, {600, 0, qnHistMax}}); + } + } + + if (doprocessData || doprocessMC) { + registry.add("trackQA/before/phi_eta_vtxZ", "", {HistType::kTH3D, {phiAxis, etaAxis, vtxAxis}}); + registry.add("trackQA/before/pt_dcaXY_dcaZ", "", {HistType::kTH3D, {ptAxis, dcaXYAXis, dcaZAXis}}); + registry.add("trackQA/before/nch_pt", "#it{p}_{T} vs multiplicity; N_{ch}; #it{p}_{T}", {HistType::kTH2D, {nchAxis, ptAxis}}); + registry.add("trackQA/before/chi2prTPCcls", "#chi^{2}/cluster for the TPC track segment; #chi^{2}/TPC cluster", {HistType::kTH1D, {{100, 0., 5.}}}); + registry.add("trackQA/before/chi2prITScls", "#chi^{2}/cluster for the ITS track; #chi^{2}/ITS cluster", {HistType::kTH1D, {{100, 0., 50.}}}); + registry.add("trackQA/before/nTPCClusters", "Number of found TPC clusters; TPC N_{cls}; Counts", {HistType::kTH1D, {{100, 40, 180}}}); + registry.add("trackQA/before/nITSClusters", "Number of found ITS clusters; ITS N_{cls}; Counts", {HistType::kTH1D, {{100, 0, 20}}}); + registry.add("trackQA/before/nTPCCrossedRows", "Number of crossed TPC Rows; TPC X-rows; Counts", {HistType::kTH1D, {{100, 40, 180}}}); + + registry.addClone("trackQA/before/", "trackQA/after/"); + registry.add("trackQA/after/pt_ref", "", {HistType::kTH1D, {{100, o2::analysis::gfwflowflucpp::ptreflow, o2::analysis::gfwflowflucpp::ptrefup}}}); + registry.add("trackQA/after/pt_poi", "", {HistType::kTH1D, {{100, o2::analysis::gfwflowflucpp::ptpoilow, o2::analysis::gfwflowflucpp::ptpoiup}}}); + + registry.add("eventQA/before/multiplicity", "", {HistType::kTH1D, {nchAxis}}); + if (cfgTimeDependent) { + registry.add("eventQA/before/multiplicity_time", "Multiplicity vs time; time (hrs); N_{ch}", {HistType::kTH2D, {timeAxis, nchAxis}}); + registry.add("eventQA/before/multT0C_time", "T0C Multiplicity vs time; time (hrs); N_{ch} (T0C)", {HistType::kTH2D, {timeAxis, t0cAxis}}); + registry.add("eventQA/before/multT0A_time", "T0A Multiplicity vs time; time (hrs); N_{ch} (T0A)", {HistType::kTH2D, {timeAxis, t0aAxis}}); + registry.add("eventQA/before/multV0A_time", "V0A Multiplicity vs time; time (hrs); N_{ch} (V0A)", {HistType::kTH2D, {timeAxis, v0aAxis}}); + registry.add("eventQA/before/multPV_time", "PV Multiplicity vs time; time (hrs); N_{ch} (PV)", {HistType::kTH2D, {timeAxis, multpvAxis}}); + } + registry.add("eventQA/before/globalTracks_PVTracks", "", {HistType::kTH2D, {multpvAxis, nchAxis}}); + registry.add("eventQA/before/globalTracks_multT0A", "", {HistType::kTH2D, {t0aAxis, nchAxis}}); + registry.add("eventQA/before/globalTracks_multV0A", "", {HistType::kTH2D, {v0aAxis, nchAxis}}); + registry.add("eventQA/before/multV0A_multT0A", "", {HistType::kTH2D, {t0aAxis, v0aAxis}}); + + registry.add("eventQA/before/centrality", "", {HistType::kTH1D, {centAxis}}); + registry.add("eventQA/before/globalTracks_centT0C", "", {HistType::kTH2D, {centAxis, nchAxis}}); + registry.add("eventQA/before/PVTracks_centT0C", "", {HistType::kTH2D, {centAxis, multpvAxis}}); + registry.add("eventQA/before/multT0C_centT0C", "", {HistType::kTH2D, {centAxis, t0cAxis}}); + + registry.add("eventQA/before/centT0M_centT0C", "", {HistType::kTH2D, {centAxis, centAxis}}); + registry.add("eventQA/before/centV0A_centT0C", "", {HistType::kTH2D, {centAxis, centAxis}}); + registry.add("eventQA/before/centGlobal_centT0C", "", {HistType::kTH2D, {centAxis, centAxis}}); + registry.add("eventQA/before/centNTPV_centT0C", "", {HistType::kTH2D, {centAxis, centAxis}}); + registry.add("eventQA/before/centMFT_centT0C", "", {HistType::kTH2D, {centAxis, centAxis}}); + + if (cfgIsMC || doprocessMC) { + registry.add("MCGen/trackQA/phi_eta_vtxZ", "", {HistType::kTH3D, {phiAxis, etaAxis, vtxAxis}}); + registry.add("MCGen/trackQA/nch_pt", "#it{p}_{T} vs multiplicity; N_{ch}; #it{p}_{T}", {HistType::kTH2D, {nchAxis, ptAxis}}); + registry.add("MCGen/trackQA/pt_ref", "", {HistType::kTH1D, {{100, o2::analysis::gfwflowflucpp::ptreflow, o2::analysis::gfwflowflucpp::ptrefup}}}); + registry.add("MCGen/trackQA/pt_poi", "", {HistType::kTH1D, {{100, o2::analysis::gfwflowflucpp::ptpoilow, o2::analysis::gfwflowflucpp::ptpoiup}}}); + } + + registry.addClone("eventQA/before/", "eventQA/after/"); + registry.add("eventQA/eventSel", "Number of Events;; Counts", {HistType::kTH1D, {{11, 0.5, 11.5}}}); + registry.get(HIST("eventQA/eventSel"))->GetXaxis()->SetBinLabel(kFilteredEvent, "Filtered event"); + registry.get(HIST("eventQA/eventSel"))->GetXaxis()->SetBinLabel(kSel8, "sel8"); + registry.get(HIST("eventQA/eventSel"))->GetXaxis()->SetBinLabel(kOccupancy, "occupancy"); + registry.get(HIST("eventQA/eventSel"))->GetXaxis()->SetBinLabel(kTVXTRD, "kTVXinTRD"); + registry.get(HIST("eventQA/eventSel"))->GetXaxis()->SetBinLabel(kNoSamebunchPU, "kNoSameBunchPileup"); + registry.get(HIST("eventQA/eventSel"))->GetXaxis()->SetBinLabel(kZVtxFT0PV, "kIsGoodZvtxFT0vsPV"); + registry.get(HIST("eventQA/eventSel"))->GetXaxis()->SetBinLabel(kNoCollTRStd, "kNoCollInTimeRangeStandard"); + registry.get(HIST("eventQA/eventSel"))->GetXaxis()->SetBinLabel(kVtxITSTPC, "kIsVertexITSTPC"); + registry.get(HIST("eventQA/eventSel"))->GetXaxis()->SetBinLabel(kGoodITSLayers, "kIsGoodITSLayersAll"); + registry.get(HIST("eventQA/eventSel"))->GetXaxis()->SetBinLabel(kMultCuts, "after Mult cuts"); + registry.get(HIST("eventQA/eventSel"))->GetXaxis()->SetBinLabel(kTrackCent, "has track + within cent"); + if (!cfgRunByRun && cfgFillWeights) { + registry.add("phi_eta_vtxz_ref", "", {HistType::kTH3D, {phiAxis, etaAxis, vtxAxis}}); + } + } + + if (o2::analysis::gfwflowflucpp::regions.GetSize() < 0) + LOGF(error, "Configuration contains vectors of different size - check the GFWRegions configurable"); + for (auto i(0); i < o2::analysis::gfwflowflucpp::regions.GetSize(); ++i) { + fGFW->AddRegion(o2::analysis::gfwflowflucpp::regions.GetNames()[i], o2::analysis::gfwflowflucpp::regions.GetEtaMin()[i], o2::analysis::gfwflowflucpp::regions.GetEtaMax()[i], (o2::analysis::gfwflowflucpp::regions.GetpTDifs()[i]) ? ptbins + 1 : 1, o2::analysis::gfwflowflucpp::regions.GetBitmasks()[i]); + } + for (auto i = 0; i < o2::analysis::gfwflowflucpp::configs.GetSize(); ++i) { + corrconfigs.push_back(fGFW->GetCorrelatorConfig(o2::analysis::gfwflowflucpp::configs.GetCorrs()[i], o2::analysis::gfwflowflucpp::configs.GetHeads()[i], o2::analysis::gfwflowflucpp::configs.GetpTDifs()[i])); + } + if (corrconfigs.empty()) + LOGF(error, "Configuration contains vectors of different size - check the GFWCorrConfig configurable"); + fGFW->CreateRegions(); + TObjArray* oba = new TObjArray(); + addConfigObjectsToObjArray(oba, corrconfigs); + fFC->SetName("FlowContainer"); + fFC->SetXAxis(fSecondAxis); + fFC->Initialize(oba, multAxis, cfgNbootstrap); + + fFCgen->SetName("FlowContainer_gen"); + fFCgen->SetXAxis(fSecondAxis); + fFCgen->Initialize(oba, multAxis, cfgNbootstrap); + + delete oba; + + fPtDepDCAxy = new TF1("ptDepDCAxy", Form("[0]*%s", cfgDCAxy->c_str()), 0.001, 100); + fPtDepDCAxy->SetParameter(0, cfgDCAxyNSigma); + LOGF(info, "DCAxy pt-dependence function: %s", Form("[0]*%s", cfgDCAxy->c_str())); + if (cfgUseAdditionalEventCut) { + fMultPVCutLow = new TF1("fMultPVCutLow", cfgMultCorrLowCutFunction->c_str(), 0, 100); + fMultPVCutLow->SetParameters(&(o2::analysis::gfwflowflucpp::multPVCorrCutPars[0])); + fMultPVCutHigh = new TF1("fMultPVCutHigh", cfgMultCorrHighCutFunction->c_str(), 0, 100); + fMultPVCutHigh->SetParameters(&(o2::analysis::gfwflowflucpp::multPVCorrCutPars[0])); + fMultCutLow = new TF1("fMultCutLow", cfgMultCorrLowCutFunction->c_str(), 0, 100); + fMultCutLow->SetParameters(&(o2::analysis::gfwflowflucpp::multGlobalCorrCutPars[0])); + fMultCutHigh = new TF1("fMultCutHigh", cfgMultCorrHighCutFunction->c_str(), 0, 100); + fMultCutHigh->SetParameters(&(o2::analysis::gfwflowflucpp::multGlobalCorrCutPars[0])); + fMultPVGlobalCutHigh = new TF1("fMultPVGlobalCutHigh", cfgMultGlobalPVCorrCutFunction->c_str(), 0, nchbinning.back()); + fMultPVGlobalCutHigh->SetParameters(&(o2::analysis::gfwflowflucpp::multGlobalPVCorrCutPars[0])); + + LOGF(info, "Global V0A function: %s in range 0-%g", cfgGlobalAsideCorrCuts.cfgMultGlobalASideCorrCutFunction->c_str(), v0aAxis.binEdges.back()); + fMultGlobalV0ACutLow = new TF1("fMultGlobalV0ACutLow", cfgGlobalAsideCorrCuts.cfgMultGlobalASideCorrCutFunction->c_str(), 0, v0aAxis.binEdges.back()); + for (std::size_t i = 0; i < o2::analysis::gfwflowflucpp::multGlobalV0ACutPars.size(); ++i) + fMultGlobalV0ACutLow->SetParameter(i, o2::analysis::gfwflowflucpp::multGlobalV0ACutPars[i]); + fMultGlobalV0ACutLow->SetParameter(o2::analysis::gfwflowflucpp::multGlobalV0ACutPars.size(), cfgGlobalAsideCorrCuts.cfgGlobalV0ALowSigma); + for (int i = 0; i < fMultGlobalV0ACutLow->GetNpar(); ++i) + LOGF(info, "fMultGlobalV0ACutLow par %d = %g", i, fMultGlobalV0ACutLow->GetParameter(i)); + + fMultGlobalV0ACutHigh = new TF1("fMultGlobalV0ACutHigh", cfgGlobalAsideCorrCuts.cfgMultGlobalASideCorrCutFunction->c_str(), 0, v0aAxis.binEdges.back()); + for (std::size_t i = 0; i < o2::analysis::gfwflowflucpp::multGlobalV0ACutPars.size(); ++i) + fMultGlobalV0ACutHigh->SetParameter(i, o2::analysis::gfwflowflucpp::multGlobalV0ACutPars[i]); + fMultGlobalV0ACutHigh->SetParameter(o2::analysis::gfwflowflucpp::multGlobalV0ACutPars.size(), cfgGlobalAsideCorrCuts.cfgGlobalV0AHighSigma); + for (int i = 0; i < fMultGlobalV0ACutHigh->GetNpar(); ++i) + LOGF(info, "fMultGlobalV0ACutHigh par %d = %g", i, fMultGlobalV0ACutHigh->GetParameter(i)); + + LOGF(info, "Global T0A function: %s", cfgGlobalAsideCorrCuts.cfgMultGlobalASideCorrCutFunction->c_str()); + fMultGlobalT0ACutLow = new TF1("fMultGlobalT0ACutLow", cfgGlobalAsideCorrCuts.cfgMultGlobalASideCorrCutFunction->c_str(), 0, t0aAxis.binEdges.back()); + for (std::size_t i = 0; i < o2::analysis::gfwflowflucpp::multGlobalT0ACutPars.size(); ++i) + fMultGlobalT0ACutLow->SetParameter(i, o2::analysis::gfwflowflucpp::multGlobalT0ACutPars[i]); + fMultGlobalT0ACutLow->SetParameter(o2::analysis::gfwflowflucpp::multGlobalT0ACutPars.size(), cfgGlobalAsideCorrCuts.cfgGlobalT0ALowSigma); + for (int i = 0; i < fMultGlobalT0ACutLow->GetNpar(); ++i) + LOGF(info, "fMultGlobalT0ACutLow par %d = %g", i, fMultGlobalT0ACutLow->GetParameter(i)); + + fMultGlobalT0ACutHigh = new TF1("fMultGlobalT0ACutHigh", cfgGlobalAsideCorrCuts.cfgMultGlobalASideCorrCutFunction->c_str(), 0, t0aAxis.binEdges.back()); + for (std::size_t i = 0; i < o2::analysis::gfwflowflucpp::multGlobalT0ACutPars.size(); ++i) + fMultGlobalT0ACutHigh->SetParameter(i, o2::analysis::gfwflowflucpp::multGlobalT0ACutPars[i]); + fMultGlobalT0ACutHigh->SetParameter(o2::analysis::gfwflowflucpp::multGlobalT0ACutPars.size(), cfgGlobalAsideCorrCuts.cfgGlobalT0AHighSigma); + for (int i = 0; i < fMultGlobalT0ACutHigh->GetNpar(); ++i) + LOGF(info, "fMultGlobalT0ACutHigh par %d = %g", i, fMultGlobalT0ACutHigh->GetParameter(i)); + } + + if (cfgConsistentEventFlag) { + const auto& names = cfgRegions->GetNames(); + auto findRegionIndex = [&](const std::string& name) { + auto it = std::find(names.begin(), names.end(), name); + return (it != names.end()) ? std::distance(names.begin(), it) : -1; + }; + posRegionIndex = findRegionIndex("refP"); + negRegionIndex = findRegionIndex("refN"); + fullRegionIndex = findRegionIndex("refFull"); + midRegionIndex = findRegionIndex("refMid"); + } + } + + static constexpr std::string_view FillTimeName[] = {"before/", "after/"}; + + enum QAFillTime { + kBefore, + kAfter + }; + + int getMagneticField(uint64_t timestamp) + { + static o2::parameters::GRPMagField* grpo = nullptr; + if (grpo == nullptr) { + // grpo = ccdb->getForTimeStamp("GLO/GRP/GRP", timestamp); + grpo = ccdb->getForTimeStamp("GLO/Config/GRPMagField", timestamp); + if (grpo == nullptr) { + LOGF(fatal, "GRP object not found for timestamp %llu", timestamp); + return 0; + } + LOGF(info, "Retrieved GRP for timestamp %llu with magnetic field of %d kG", timestamp, grpo->getNominalL3Field()); + } + return grpo->getNominalL3Field(); + } + + void loadCorrections(aod::BCsWithTimestamps::iterator const& bc) + { + uint64_t timestamp = bc.timestamp(); + if (!cfgRunByRun && cfg.correctionsLoaded) + return; + if (!cfgAcceptance.value.empty()) { + std::string runstr = (cfgRunByRun) ? "RunByRun/" : ""; + cfg.mAcceptance = ccdb->getForTimeStamp(cfgAcceptance.value + runstr, timestamp); + } + if (!cfgEfficiency.value.empty()) { + cfg.mEfficiency = ccdb->getForTimeStamp(cfgEfficiency, timestamp); + if (cfg.mEfficiency == nullptr) { + LOGF(fatal, "Could not load efficiency histogram from %s", cfgEfficiency.value.c_str()); + } + LOGF(info, "Loaded efficiency histogram from %s (%p)", cfgEfficiency.value.c_str(), (void*)cfg.mEfficiency); + } + cfg.correctionsLoaded = true; + } + + template + double getAcceptance(TTrack track, const double& vtxz) + { + double wacc = 1; + if (cfg.mAcceptance) + wacc = cfg.mAcceptance->getNUA(track.phi(), track.eta(), vtxz); + return wacc; + } + + template + double getEfficiency(TTrack track) + { + double eff = 1.; + if (cfg.mEfficiency) + eff = cfg.mEfficiency->GetBinContent(cfg.mEfficiency->FindBin(track.pt())); + if (eff == 0) + return -1.; + else + return 1. / eff; + } + + template + bool eventSelected(TCollision collision, const int& multTrk, const float& centrality, const int& run) + { + if (cfgTVXinTRD) { + if (collision.alias_bit(kTVXinTRD)) { + // TRD triggered + // "CMTVX-B-NOPF-TRD,minbias_TVX" + return 0; + } + if (cfgFillQA) { + registry.fill(HIST("eventQA/eventSel"), kTVXTRD); + } + if (cfgRunByRun && run != -1) + th1sList[run][hEventSel]->Fill(kTVXTRD); + } + if (cfgNoSameBunchPileupCut) { + if (!collision.selection_bit(o2::aod::evsel::kNoSameBunchPileup)) { + // rejects collisions which are associated with the same "found-by-T0" bunch crossing + // https://indico.cern.ch/event/1396220/#1-event-selection-with-its-rof + return 0; + } + if (cfgFillQA) { + registry.fill(HIST("eventQA/eventSel"), kNoSamebunchPU); + } + if (cfgRunByRun && run != -1) + th1sList[run][hEventSel]->Fill(kNoSamebunchPU); + } + if (cfgIsGoodZvtxFT0vsPV) { + if (!collision.selection_bit(o2::aod::evsel::kIsGoodZvtxFT0vsPV)) { + // removes collisions with large differences between z of PV by tracks and z of PV from FT0 A-C time difference + // use this cut at low multiplicities with caution + return 0; + } + if (cfgFillQA) { + registry.fill(HIST("eventQA/eventSel"), kZVtxFT0PV); + } + if (cfgRunByRun && run != -1) + th1sList[run][hEventSel]->Fill(kZVtxFT0PV); + } + if (cfgNoCollInTimeRangeStandard) { + if (!collision.selection_bit(o2::aod::evsel::kNoCollInTimeRangeStandard)) { + // Rejection of the collisions which have other events nearby + return 0; + } + if (cfgFillQA) { + registry.fill(HIST("eventQA/eventSel"), kNoCollTRStd); + } + if (cfgRunByRun && run != -1) + th1sList[run][hEventSel]->Fill(kNoCollTRStd); + } + + if (cfgIsVertexITSTPC) { + if (!collision.selection_bit(o2::aod::evsel::kIsVertexITSTPC)) { + // selects collisions with at least one ITS-TPC track, and thus rejects vertices built from ITS-only tracks + return 0; + } + if (cfgFillQA) { + registry.fill(HIST("eventQA/eventSel"), kVtxITSTPC); + } + if (cfgRunByRun && run != -1) + th1sList[run][hEventSel]->Fill(kVtxITSTPC); + } + + if (cfgIsGoodITSLayersAll) { + if (!collision.selection_bit(o2::aod::evsel::kIsGoodITSLayersAll)) { + return 0; + } + if (cfgFillQA) { + registry.fill(HIST("eventQA/eventSel"), kGoodITSLayers); + } + if (cfgRunByRun && run != -1) + th1sList[run][hEventSel]->Fill(kGoodITSLayers); + } + + float vtxz = -999; + if (collision.numContrib() > 1) { + vtxz = collision.posZ(); + float zRes = std::sqrt(collision.covZZ()); + float minZRes = 0.25; + int minNContrib = 20; + if (zRes > minZRes && collision.numContrib() < minNContrib) + vtxz = -999; + } + auto multNTracksPV = collision.multNTracksPV(); + + if (vtxz > o2::analysis::gfwflowflucpp::vtxZup || vtxz < o2::analysis::gfwflowflucpp::vtxZlow) + return 0; + + if (cfgMultCut && cfgUseAdditionalEventCut) { + if (multNTracksPV < fMultPVCutLow->Eval(centrality)) + return 0; + if (multNTracksPV > fMultPVCutHigh->Eval(centrality)) + return 0; + if (multTrk < fMultCutLow->Eval(centrality)) + return 0; + if (multTrk > fMultCutHigh->Eval(centrality)) + return 0; + if (multTrk > fMultPVGlobalCutHigh->Eval(collision.multNTracksPV())) + return 0; + + if (!(cfgGlobalAsideCorrCuts.cfgMultGlobalASideCorrCutFunction->empty()) && static_cast(collision.multFV0A()) < fMultGlobalV0ACutLow->Eval(multTrk)) + return 0; + if (!(cfgGlobalAsideCorrCuts.cfgMultGlobalASideCorrCutFunction->empty()) && static_cast(collision.multFV0A()) > fMultGlobalV0ACutHigh->Eval(multTrk)) + return 0; + if (!(cfgGlobalAsideCorrCuts.cfgMultGlobalASideCorrCutFunction->empty()) && static_cast(collision.multFT0A()) < fMultGlobalT0ACutLow->Eval(multTrk)) + return 0; + if (!(cfgGlobalAsideCorrCuts.cfgMultGlobalASideCorrCutFunction->empty()) && static_cast(collision.multFT0A()) > fMultGlobalT0ACutHigh->Eval(multTrk)) + return 0; + if (cfgRunByRun && run != -1) + th1sList[run][hEventSel]->Fill(kMultCuts); + } + if (cfgFillQA) { + registry.fill(HIST("eventQA/eventSel"), kMultCuts); + } + return 1; + } + + template + bool trackSelected(TTrack track) + { + if (cfgDCAxyNSigma && (std::fabs(track.dcaXY()) > fPtDepDCAxy->Eval(track.pt()))) + return false; + return ((track.tpcNClsCrossedRows() >= cfgNTPCXrows) && (track.tpcNClsFound() >= cfgNTPCCls) && (track.itsNCls() >= cfgMinNITSCls)); + } + + enum DataType { + kReco, + kGen + }; + + template + void fillWeights(const TTrack track, const double vtxz, const int& run) + { + if (cfgRunByRun) + th3sList[run][hNUAref]->Fill(track.phi(), track.eta(), vtxz); + else + registry.fill(HIST("phi_eta_vtxz_ref"), track.phi(), track.eta(), vtxz); + return; + } + + void createRunByRunHistograms(const int& run) + { + AxisSpec phiAxis = {o2::analysis::gfwflowflucpp::phibins, o2::analysis::gfwflowflucpp::philow, o2::analysis::gfwflowflucpp::phiup, "#phi"}; + AxisSpec etaAxis = {o2::analysis::gfwflowflucpp::etabins, -cfgEta, cfgEta, "#eta"}; + AxisSpec vtxAxis = {o2::analysis::gfwflowflucpp::vtxZbins, -cfgVtxZ, cfgVtxZ, "Vtx_{z} (cm)"}; + AxisSpec nchAxis = {o2::analysis::gfwflowflucpp::nchbins, o2::analysis::gfwflowflucpp::nchlow, o2::analysis::gfwflowflucpp::nchup, "N_{ch}"}; + AxisSpec centAxis = {o2::analysis::gfwflowflucpp::centbinning, "Centrality (%)"}; + std::vector> histos(kCount_TH1Names); + histos[hPhi] = registry.add(Form("%d/phi", run), "", {HistType::kTH1D, {phiAxis}}); + histos[hEta] = registry.add(Form("%d/eta", run), "", {HistType::kTH1D, {etaAxis}}); + histos[hVtxZ] = registry.add(Form("%d/vtxz", run), "", {HistType::kTH1D, {vtxAxis}}); + histos[hMult] = registry.add(Form("%d/mult", run), "", {HistType::kTH1D, {nchAxis}}); + histos[hCent] = registry.add(Form("%d/cent", run), "", {HistType::kTH1D, {centAxis}}); + if (cfgFillFlowRunByRun) { + std::vector> profiles(kCount_TProfileNames); + profiles[pfCorr22] = registry.add(Form("%d/corr22", run), "", {HistType::kTProfile, {(cfgUseNch) ? nchAxis : centAxis}}); + tpfsList.insert(std::make_pair(run, profiles)); + } + histos[hEventSel] = registry.add(Form("%d/eventSel", run), "Number of Events;; Counts", {HistType::kTH1D, {{11, 0.5, 11.5}}}); + histos[hEventSel]->GetXaxis()->SetBinLabel(kFilteredEvent, "Filtered event"); + histos[hEventSel]->GetXaxis()->SetBinLabel(kSel8, "sel8"); + histos[hEventSel]->GetXaxis()->SetBinLabel(kOccupancy, "occupancy"); + histos[hEventSel]->GetXaxis()->SetBinLabel(kTVXTRD, "kTVXinTRD"); + histos[hEventSel]->GetXaxis()->SetBinLabel(kNoSamebunchPU, "kNoSameBunchPileup"); + histos[hEventSel]->GetXaxis()->SetBinLabel(kZVtxFT0PV, "kIsGoodZvtxFT0vsPV"); + histos[hEventSel]->GetXaxis()->SetBinLabel(kNoCollTRStd, "kNoCollInTimeRangeStandard"); + histos[hEventSel]->GetXaxis()->SetBinLabel(kVtxITSTPC, "kIsVertexITSTPC"); + histos[hEventSel]->GetXaxis()->SetBinLabel(kGoodITSLayers, "kIsGoodITSLayersAll"); + histos[hEventSel]->GetXaxis()->SetBinLabel(kMultCuts, "after Mult cuts"); + histos[hEventSel]->GetXaxis()->SetBinLabel(kTrackCent, "has track + within cent"); + th1sList.insert(std::make_pair(run, histos)); + std::vector> histos3d(kCount_TH3Names); + histos3d[hNUAref] = registry.add(Form("%d/phi_eta_vtxz_ref", run), "", {HistType::kTH3D, {phiAxis, etaAxis, vtxAxis}}); + th3sList.insert(std::make_pair(run, histos3d)); + return; + } + + void addConfigObjectsToObjArray(TObjArray* oba, const std::vector& configs) + { + const auto shapeSel = getShapeSel(); + const int nQnOutputBins = getNQnOutputBins(); + for (auto it = configs.begin(); it != configs.end(); ++it) { + for (int jese = 0; jese < nQnOutputBins; ++jese) { + std::string name = Form("%s_%d_%s", shapeSel.c_str(), jese, it->Head.c_str()); + std::string title = it->Head + std::string("_ese"); + oba->Add(new TNamed(name.c_str(), title.c_str())); + } + } + } + + template + void fillOutputContainers(const float& centmult, const double& rndm, const int& qPtmp, const int& run = 0) + { + const auto shapeSel = getShapeSel(); + + for (uint l_ind = 0; l_ind < corrconfigs.size(); ++l_ind) { + auto dnx = fGFW->Calculate(corrconfigs.at(l_ind), 0, kTRUE).real(); + if (dnx == 0) + continue; + + auto val = fGFW->Calculate(corrconfigs.at(l_ind), 0, kFALSE).real() / dnx; + if (std::abs(val) >= 1.) + continue; + + std::string profName = Form("%s_%d_%s", shapeSel.c_str(), qPtmp, corrconfigs.at(l_ind).Head.c_str()); + + if constexpr (dt == kGen) { + fFCgen->FillProfile(profName.c_str(), + centmult, + val, + cfgUseMultiplicityFlowWeights ? dnx : 1.0, + rndm); + } else { + fFC->FillProfile(profName.c_str(), + centmult, + val, + cfgUseMultiplicityFlowWeights ? dnx : 1.0, + rndm); + + if (cfgRunByRun && cfgFillFlowRunByRun && l_ind == 0) { + tpfsList[run][pfCorr22]->Fill(centmult, + val, + cfgUseMultiplicityFlowWeights ? dnx : 1.0); + } + } + } + } + + struct XAxis { + float centrality; + int64_t multiplicity; + double time; + }; + + struct AcceptedTracks { + int nPos; + int nNeg; + int nFull; + int nMid; + }; + + 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 (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 && fabs(trk.eta())< cfgQnTrkAbsEtaMax) { + // 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 && fabs(trk.eta())< cfgQnTrkAbsEtaMax) { + // 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 (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]); + } + + return qvec; + } + + float computeqnVec(InTaskTPCQVector const& qvec, bool useNegativeEtaHalf) + { + if (qvec.nTrkTPCPos <= 0 || qvec.nTrkTPCNeg <= 0) { + return -1.f; + } + + 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 + std::vector> getQnBinSeparator2D(std::vector flat, const int numQnBins = 10) + { + size_t nBins = numQnBins + 1; + + if (flat.empty() || flat.size() % nBins != 0) { + LOGP(error, "ConfQnBinSeparator size = {} is not divisible by {}", + flat.size(), nBins); + return {{kInvalidQnSeparator, kInvalidQnSeparator}}; + } + + size_t nCent = flat.size() / nBins; + std::vector> res(nCent, std::vector(nBins)); + + for (size_t i = 0; i < nCent; ++i) { + for (size_t j = 0; j < nBins; ++j) { + res[i][j] = flat[i * nBins + j]; + } + } + return res; + } + + /// Get the bin number of qn-vector(FT0C) of an event + /// \param centBinWidth centrality bin width, example: per 1%, per 10% ... + /// \return bin number of qn-vector of the event + // add a param : bool doFillHisto ? + int myqnBin(float centrality, float centMax, float qn, std::vector qnBinSprt, const int numQnBins, float centBinWidth = 1.f) + { + auto twoDSeparator = getQnBinSeparator2D(qnBinSprt, numQnBins); + if (twoDSeparator.empty() || twoDSeparator[0][0] == kInvalidQnSeparator) { + LOGP(warning, "ConfQnBinSeparator not set, using default fallback!"); + return kInvalidQnBin; // safe fallback + } + + int qnBin = kInvalidQnBin; + int mycentBin = static_cast(centrality / centBinWidth); + if (mycentBin >= static_cast(centMax / centBinWidth)) + return qnBin; + + if (mycentBin > static_cast(twoDSeparator.size()) - 1) + return qnBin; + + for (int iqn(0); iqn < static_cast(twoDSeparator[mycentBin].size()) - 1; ++iqn) { + if (qn > twoDSeparator[mycentBin][iqn] && qn <= twoDSeparator[mycentBin][iqn + 1]) { + qnBin = iqn; + break; + } else { + continue; + } + } + + return qnBin; + } + + template + void processCollision(TCollision collision, TTracks tracks, const XAxis& xaxis, const int& run, const int& qPtmp) + { + if (tracks.size() < 1) + return; + if (dt != kGen && xaxis.centrality >= 0 && + (xaxis.centrality < o2::analysis::gfwflowflucpp::centbinning.front() || + xaxis.centrality > o2::analysis::gfwflowflucpp::centbinning.back())) + return; + if (xaxis.multiplicity < cfgFixedMultMin || xaxis.multiplicity > cfgFixedMultMax) + return; + + if (dt != kGen) { + if (cfgFillQA) { + registry.fill(HIST("eventQA/eventSel"), kTrackCent); + } + if (cfgRunByRun) + th1sList[run][hEventSel]->Fill(kTrackCent); + } + + if (cfgFillQA && xaxis.centrality >= 0) + registry.fill(HIST("eventQA/after/centrality"), xaxis.centrality); + if (cfgFillQA) + registry.fill(HIST("eventQA/after/multiplicity"), xaxis.multiplicity); + + float vtxz = collision.posZ(); + if (dt != kGen && cfgRunByRun) { + th1sList[run][hVtxZ]->Fill(vtxz); + th1sList[run][hMult]->Fill(xaxis.multiplicity); + th1sList[run][hCent]->Fill(xaxis.centrality); + } + + fGFW->Clear(); + float lRandom = fRndm->Rndm(); + + AcceptedTracks acceptedTracks{0, 0, 0, 0}; + + for (const auto& track : tracks) { + processTrack(track, vtxz, xaxis.multiplicity, run, acceptedTracks); + } + + if (cfgConsistentEventFlag & kRequireBothEtaSides) + if (!acceptedTracks.nPos || !acceptedTracks.nNeg) + return; + if (cfgConsistentEventFlag & kRequireFullFourParticleTracks) + if (acceptedTracks.nFull < kMinTracksForFourParticleCorrelation) + return; + if (cfgConsistentEventFlag & kRequireTwoTracksInBothEtaSides) + if (acceptedTracks.nPos < kMinTracksPerEtaSideForGapCorrelation || + acceptedTracks.nNeg < kMinTracksPerEtaSideForGapCorrelation) + return; + if (cfgConsistentEventFlag & kRequireTwoTracksInThreeEtaRegions) + if (acceptedTracks.nPos < kMinTracksPerEtaRegionForThreeSubevents || + acceptedTracks.nMid < kMinTracksPerEtaRegionForThreeSubevents || + acceptedTracks.nNeg < kMinTracksPerEtaRegionForThreeSubevents) + return; + + fillOutputContainers
(cfgUseNch ? static_cast(xaxis.multiplicity) : xaxis.centrality, + lRandom, qPtmp, run); + } + + template + void processGenCollision(TCollision collision, TParticles particles, const int& mcCollisionId, const XAxis& xaxis, const int& run, const int& qPtmp) + { + if (xaxis.multiplicity < cfgFixedMultMin || xaxis.multiplicity > cfgFixedMultMax) + return; + + if (cfgFillQA && xaxis.centrality >= 0) + registry.fill(HIST("eventQA/after/centrality"), xaxis.centrality); + if (cfgFillQA) + registry.fill(HIST("eventQA/after/multiplicity"), xaxis.multiplicity); + + fGFW->Clear(); + float lRandom = fRndm->Rndm(); + float vtxz = collision.posZ(); + + AcceptedTracks acceptedTracks{0, 0, 0, 0}; + for (const auto& particle : particles) { + if (particle.mcCollisionId() != mcCollisionId) + continue; + processTrack(particle, vtxz, xaxis.multiplicity, run, acceptedTracks); + } + + if (cfgConsistentEventFlag & kRequireBothEtaSides) + if (!acceptedTracks.nPos || !acceptedTracks.nNeg) + return; + if (cfgConsistentEventFlag & kRequireFullFourParticleTracks) + if (acceptedTracks.nFull < kMinTracksForFourParticleCorrelation) + return; + if (cfgConsistentEventFlag & kRequireTwoTracksInBothEtaSides) + if (acceptedTracks.nPos < kMinTracksPerEtaSideForGapCorrelation || + acceptedTracks.nNeg < kMinTracksPerEtaSideForGapCorrelation) + return; + if (cfgConsistentEventFlag & kRequireTwoTracksInThreeEtaRegions) + if (acceptedTracks.nPos < kMinTracksPerEtaRegionForThreeSubevents || + acceptedTracks.nMid < kMinTracksPerEtaRegionForThreeSubevents || + acceptedTracks.nNeg < kMinTracksPerEtaRegionForThreeSubevents) + return; + + fillOutputContainers(cfgUseNch ? static_cast(xaxis.multiplicity) : xaxis.centrality, + lRandom, qPtmp, run); + } + + bool isStable(int pdg) + { + if (std::abs(pdg) == PDG_t::kPiPlus) + return true; + if (std::abs(pdg) == PDG_t::kKPlus) + return true; + if (std::abs(pdg) == PDG_t::kProton) + return true; + if (std::abs(pdg) == PDG_t::kElectron) + return true; + if (std::abs(pdg) == PDG_t::kMuonMinus) + return true; + return false; + } + + template + void fillAcceptedTracks(TTrack track, AcceptedTracks& acceptedTracks) + { + if (posRegionIndex >= 0 && track.eta() > o2::analysis::gfwflowflucpp::regions.GetEtaMin()[posRegionIndex] && track.eta() < o2::analysis::gfwflowflucpp::regions.GetEtaMax()[posRegionIndex]) + ++acceptedTracks.nPos; + if (negRegionIndex >= 0 && track.eta() > o2::analysis::gfwflowflucpp::regions.GetEtaMin()[negRegionIndex] && track.eta() < o2::analysis::gfwflowflucpp::regions.GetEtaMax()[negRegionIndex]) + ++acceptedTracks.nNeg; + if (fullRegionIndex >= 0 && track.eta() > o2::analysis::gfwflowflucpp::regions.GetEtaMin()[fullRegionIndex] && track.eta() < o2::analysis::gfwflowflucpp::regions.GetEtaMax()[fullRegionIndex]) + ++acceptedTracks.nFull; + if (midRegionIndex >= 0 && track.eta() > o2::analysis::gfwflowflucpp::regions.GetEtaMin()[midRegionIndex] && track.eta() < o2::analysis::gfwflowflucpp::regions.GetEtaMax()[midRegionIndex]) + ++acceptedTracks.nMid; + } + + template + inline void processTrack(TTrack const& track, const float& vtxz, const int& multiplicity, const int& run, AcceptedTracks& acceptedTracks) + { + if constexpr (framework::has_type_v) { + if (track.mcParticleId() < 0 || !(track.has_mcParticle())) + return; + + auto mcParticle = track.mcParticle(); + if (!mcParticle.isPhysicalPrimary()) + return; + if (!isStable(mcParticle.pdgCode())) + return; + if (cfgFillQA && cfgIsMC) { + fillTrackQA(track, vtxz); + registry.fill(HIST("trackQA/before/nch_pt"), multiplicity, track.pt()); + } + if (!trackSelected(track)) + return; + + if (cfgFillWeights) { + fillWeights(track, vtxz, run); + } else { + fillGFW(track, vtxz); + fillAcceptedTracks(track, acceptedTracks); + } + + if (cfgFillQA) { + fillTrackQA(track, vtxz); + registry.fill(HIST("trackQA/after/nch_pt"), multiplicity, track.pt()); + if (cfgRunByRun) { + th1sList[run][hPhi]->Fill(track.phi()); + th1sList[run][hEta]->Fill(track.eta()); + } + } + + } else if constexpr (framework::has_type_v) { + if (!track.isPhysicalPrimary() || !isStable(track.pdgCode())) + return; + + fillGFW(track, vtxz); + fillAcceptedTracks(track, acceptedTracks); + if (cfgFillQA && cfgIsMC) { + fillTrackQA(track, vtxz); + registry.fill(HIST("MCGen/trackQA/nch_pt"), multiplicity, track.pt()); + } + } else { + if (cfgFillQA) { + fillTrackQA(track, vtxz); + registry.fill(HIST("trackQA/before/nch_pt"), multiplicity, track.pt()); + } + if (!trackSelected(track)) + return; + + if (cfgFillWeights) { + fillWeights(track, vtxz, run); + } else { + fillGFW(track, vtxz); + fillAcceptedTracks(track, acceptedTracks); + } + if (cfgFillQA) { + fillTrackQA(track, vtxz); + registry.fill(HIST("trackQA/after/nch_pt"), multiplicity, track.pt()); + if (cfgRunByRun) { + th1sList[run][hPhi]->Fill(track.phi()); + th1sList[run][hEta]->Fill(track.eta()); + th3sList[run][hNUAref]->Fill(track.phi(), track.eta(), vtxz, getAcceptance(track, vtxz)); + } + } + } + return; + } + + template + inline void fillGFW(TTrack track, const double& vtxz) + { + bool withinPtRef = (track.pt() > o2::analysis::gfwflowflucpp::ptreflow && track.pt() < o2::analysis::gfwflowflucpp::ptrefup); + bool withinPtPOI = (track.pt() > o2::analysis::gfwflowflucpp::ptpoilow && track.pt() < o2::analysis::gfwflowflucpp::ptpoiup); + if (!withinPtPOI && !withinPtRef) + return; + double weff = (dt == kGen) ? 1. : getEfficiency(track); + if (weff < 0) + return; + + double wacc = (dt == kGen) ? 1. : getAcceptance(track, vtxz); + if (withinPtRef) + fGFW->Fill(track.eta(), fSecondAxis->FindBin(track.pt()) - 1, track.phi(), weff * wacc, 1); + if (withinPtPOI) + fGFW->Fill(track.eta(), fSecondAxis->FindBin(track.pt()) - 1, track.phi(), weff * wacc, 2); + if (withinPtRef && withinPtPOI) + fGFW->Fill(track.eta(), fSecondAxis->FindBin(track.pt()) - 1, track.phi(), weff * wacc, 4); + return; + } + + template + inline void fillTrackQA(TTrack track, const float vtxz) + { + if constexpr (dt == kGen) { + registry.fill(HIST("MCGen/trackQA/phi_eta_vtxZ"), track.phi(), track.eta(), vtxz); + registry.fill(HIST("MCGen/trackQA/pt_ref"), track.pt()); + registry.fill(HIST("MCGen/trackQA/pt_poi"), track.pt()); + } else { + double wacc = getAcceptance(track, vtxz); + registry.fill(HIST("trackQA/") + HIST(FillTimeName[ft]) + HIST("phi_eta_vtxZ"), track.phi(), track.eta(), vtxz, (ft == kAfter) ? wacc : 1.0); + registry.fill(HIST("trackQA/") + HIST(FillTimeName[ft]) + HIST("pt_dcaXY_dcaZ"), track.pt(), track.dcaXY(), track.dcaZ()); + + registry.fill(HIST("trackQA/") + HIST(FillTimeName[ft]) + HIST("chi2prTPCcls"), track.tpcChi2NCl()); + registry.fill(HIST("trackQA/") + HIST(FillTimeName[ft]) + HIST("chi2prITScls"), track.itsChi2NCl()); + registry.fill(HIST("trackQA/") + HIST(FillTimeName[ft]) + HIST("nTPCClusters"), track.tpcNClsFound()); + registry.fill(HIST("trackQA/") + HIST(FillTimeName[ft]) + HIST("nITSClusters"), track.itsNCls()); + registry.fill(HIST("trackQA/") + HIST(FillTimeName[ft]) + HIST("nTPCCrossedRows"), track.tpcNClsCrossedRows()); + + if (ft == kAfter) { + registry.fill(HIST("trackQA/") + HIST(FillTimeName[ft]) + HIST("pt_ref"), track.pt()); + registry.fill(HIST("trackQA/") + HIST(FillTimeName[ft]) + HIST("pt_poi"), track.pt()); + } + } + } + + template + float getCentrality(TCollision collision) + { + switch (cfgCentEstimator) { + case kCentFT0C: + return collision.centFT0C(); + case kCentFT0CVariant1: + return collision.centFT0CVariant1(); + case kCentFT0M: + return collision.centFT0M(); + case kCentFV0A: + return collision.centFV0A(); + case kCentNTPV: + return collision.centNTPV(); + case kCentNGlobal: + return collision.centNGlobal(); + case kCentMFT: + return collision.centMFT(); + default: + return collision.centFT0C(); + } + } + + template + inline void fillEventQA(TCollision collision, XAxis xaxis) + { + if constexpr (framework::has_type_v) { + registry.fill(HIST("eventQA/") + HIST(FillTimeName[ft]) + HIST("globalTracks_centT0C"), collision.centFT0C(), xaxis.multiplicity); + registry.fill(HIST("eventQA/") + HIST(FillTimeName[ft]) + HIST("PVTracks_centT0C"), collision.centFT0C(), collision.multNTracksPV()); + registry.fill(HIST("eventQA/") + HIST(FillTimeName[ft]) + HIST("multT0C_centT0C"), collision.centFT0C(), collision.multFT0C()); + registry.fill(HIST("eventQA/") + HIST(FillTimeName[ft]) + HIST("centT0M_centT0C"), collision.centFT0C(), collision.centFT0M()); + registry.fill(HIST("eventQA/") + HIST(FillTimeName[ft]) + HIST("centV0A_centT0C"), collision.centFT0C(), collision.centFV0A()); + registry.fill(HIST("eventQA/") + HIST(FillTimeName[ft]) + HIST("centGlobal_centT0C"), collision.centFT0C(), collision.centNGlobal()); + registry.fill(HIST("eventQA/") + HIST(FillTimeName[ft]) + HIST("centNTPV_centT0C"), collision.centFT0C(), collision.centNTPV()); + registry.fill(HIST("eventQA/") + HIST(FillTimeName[ft]) + HIST("centMFT_centT0C"), collision.centFT0C(), collision.centMFT()); + } + registry.fill(HIST("eventQA/") + HIST(FillTimeName[ft]) + HIST("globalTracks_PVTracks"), collision.multNTracksPV(), xaxis.multiplicity); + registry.fill(HIST("eventQA/") + HIST(FillTimeName[ft]) + HIST("globalTracks_multT0A"), collision.multFT0A(), xaxis.multiplicity); + registry.fill(HIST("eventQA/") + HIST(FillTimeName[ft]) + HIST("globalTracks_multV0A"), collision.multFV0A(), xaxis.multiplicity); + registry.fill(HIST("eventQA/") + HIST(FillTimeName[ft]) + HIST("multV0A_multT0A"), collision.multFT0A(), collision.multFV0A()); + if (cfgTimeDependent) { + registry.fill(HIST("eventQA/") + HIST(FillTimeName[ft]) + HIST("multiplicity_time"), xaxis.time, xaxis.multiplicity); + registry.fill(HIST("eventQA/") + HIST(FillTimeName[ft]) + HIST("multT0C_time"), xaxis.time, collision.multFT0C()); + registry.fill(HIST("eventQA/") + HIST(FillTimeName[ft]) + HIST("multT0A_time"), xaxis.time, collision.multFT0A()); + registry.fill(HIST("eventQA/") + HIST(FillTimeName[ft]) + HIST("multV0A_time"), xaxis.time, collision.multFV0A()); + registry.fill(HIST("eventQA/") + HIST(FillTimeName[ft]) + HIST("multPV_time"), xaxis.time, collision.multNTracksPV()); + } + return; + } + + double getTimeSinceStartOfFill(uint64_t timestamp, int firstRun) + { + auto runDuration = ccdb->getRunDuration(firstRun); + uint64_t tsSOF = runDuration.first; + uint64_t diff = timestamp - tsSOF; + return static_cast(diff) / 3600000.0; + } + + void fillQnEventCounter(float count) + { + const int qnHarmonic = static_cast(cfgQnSelectionHarmonic); + + if (qnHarmonic == EllipticQVectorHarmonic) { + registry.fill(HIST("mq2/eventcounter"), count); + } else if (qnHarmonic == TriangularQVectorHarmonic) { + registry.fill(HIST("mq3/eventcounter"), count); + } + } + + void fillQnCalibrationHistograms(float centrality, float multiplicity, float qnPos, float qnNeg) + { + const int qnHarmonic = static_cast(cfgQnSelectionHarmonic); + + if (qnHarmonic == EllipticQVectorHarmonic) { + registry.fill(HIST("mq2/h2_cent_q2_etapos"), centrality, qnPos); + registry.fill(HIST("mq2/h2_cent_q2_etaneg"), centrality, qnNeg); + registry.fill(HIST("mq2/h2_mult_q2_etapos"), multiplicity, qnPos); + registry.fill(HIST("mq2/h2_mult_q2_etaneg"), multiplicity, qnNeg); + } else if (qnHarmonic == TriangularQVectorHarmonic) { + registry.fill(HIST("mq3/h2_cent_q3_etapos"), centrality, qnPos); + registry.fill(HIST("mq3/h2_cent_q3_etaneg"), centrality, qnNeg); + registry.fill(HIST("mq3/h2_mult_q3_etapos"), multiplicity, qnPos); + registry.fill(HIST("mq3/h2_mult_q3_etaneg"), multiplicity, qnNeg); + } + } + + void processData(soa::Filtered>::iterator const& collision, + aod::BCsWithTimestamps const&, GFWTracks const& tracks) + { + auto bc = collision.bc_as(); + int run = bc.runNumber(); + if (run != lastRun) { + lastRun = run; + LOGF(info, "run = %d", run); + if (cfgRunByRun) { + if (std::find(runNumbers.begin(), runNumbers.end(), run) == runNumbers.end()) { + LOGF(info, "Creating histograms for run %d", run); + createRunByRunHistograms(run); + runNumbers.push_back(run); + } + if (!cfgFillWeights) + loadCorrections(bc); + } + } + if (!cfgFillWeights && !cfgRunByRun) + loadCorrections(bc); + + registry.fill(HIST("eventQA/eventSel"), kFilteredEvent); + if (cfgRunByRun) + th1sList[run][hEventSel]->Fill(kFilteredEvent); + + if (!collision.sel8()) + return; + + registry.fill(HIST("eventQA/eventSel"), kSel8); + if (cfgRunByRun) + th1sList[run][hEventSel]->Fill(kSel8); + + if (cfgDoOccupancySel) { + int occupancy = collision.trackOccupancyInTimeRange(); + if (occupancy < 0 || occupancy > cfgOccupancySelection) + return; + } + + registry.fill(HIST("eventQA/eventSel"), kOccupancy); + if (cfgRunByRun) + th1sList[run][hEventSel]->Fill(kOccupancy); + + const XAxis xaxis{ + getCentrality(collision), + collision.multNTracksPV(), + (cfgTimeDependent) ? getTimeSinceStartOfFill(bc.timestamp(), *firstRunOfCurrentFill) : -1.0}; + + if (cfgTimeDependent && run == *firstRunOfCurrentFill && + firstRunOfCurrentFill != o2::analysis::gfwflowflucpp::firstRunsOfFill.end() - 1) + ++firstRunOfCurrentFill; + + if (cfgFillQA) + fillEventQA(collision, xaxis); + + registry.fill(HIST("eventQA/before/centrality"), xaxis.centrality); + registry.fill(HIST("eventQA/before/multiplicity"), xaxis.multiplicity); + + if (!eventSelected(collision, xaxis.multiplicity, xaxis.centrality, run)) + return; + + if (cfgFillQA) + fillEventQA(collision, xaxis); + + int qPtmp = 0; + if (!static_cast(cfgBypassQnSelection)) { + const auto qvecTPC = calcQVec(static_cast(cfgQnSelectionHarmonic), tracks, collision); + float qn = computeqnVec(qvecTPC, cfgUseNegativeEtaHalfForq2); + if (qn < 0) + return; + + qPtmp = myqnBin(cfgEvtSelCent ? xaxis.centrality : xaxis.multiplicity, + cfgCentMax, qn, qnBinSeparator, cfgNumQnBins); + if (qPtmp < 0) + return; + } + + processCollision(collision, tracks, xaxis, run, qPtmp); + } + PROCESS_SWITCH(FlowFlucGfwPp, processData, "Process analysis for non-derived data", false); + + void processq2(soa::Filtered>::iterator const& collision, aod::BCsWithTimestamps const&, GFWTracks const& tracks) + { + float count{0.5}; + fillQnEventCounter(count++); + if (!collision.sel8()) { + return; + } + fillQnEventCounter(count++); + if (cfgDoOccupancySel) { + int occupancy = collision.trackOccupancyInTimeRange(); + if (occupancy < 0 || occupancy > cfgOccupancySelection) { + return; + } + } + fillQnEventCounter(count++); + + 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 qvecTPC = calcQVec(static_cast(cfgQnSelectionHarmonic), 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; + } + + fillQnEventCounter(count++); + fillQnCalibrationHistograms(centr, multi, qvecPos, qvecNeg); + } + PROCESS_SWITCH(FlowFlucGfwPp, processq2, "Process analysis for filling q_n-vector calibration histograms", true); + + void processMC(soa::Filtered>::iterator const& collision, + aod::BCsWithTimestamps const&, GFWTracksMC const& tracks, aod::McCollisions const&, aod::McParticles const& mcParticles) + { + auto bc = collision.bc_as(); + int run = bc.runNumber(); + if (run != lastRun) { + lastRun = run; + LOGF(info, "run = %d", run); + if (cfgRunByRun) { + if (std::find(runNumbers.begin(), runNumbers.end(), run) == runNumbers.end()) { + LOGF(info, "Creating histograms for run %d", run); + createRunByRunHistograms(run); + runNumbers.push_back(run); + } + if (!cfgFillWeights) + loadCorrections(bc); + } + } + if (!cfgFillWeights && !cfgRunByRun) + loadCorrections(bc); + + registry.fill(HIST("eventQA/eventSel"), kFilteredEvent); + if (cfgRunByRun) + th1sList[run][hEventSel]->Fill(kFilteredEvent); + + if (!collision.sel8()) + return; + + registry.fill(HIST("eventQA/eventSel"), kSel8); + if (cfgRunByRun) + th1sList[run][hEventSel]->Fill(kSel8); + + if (cfgDoOccupancySel) { + int occupancy = collision.trackOccupancyInTimeRange(); + if (occupancy < 0 || occupancy > cfgOccupancySelection) + return; + } + + registry.fill(HIST("eventQA/eventSel"), kOccupancy); + if (cfgRunByRun) + th1sList[run][hEventSel]->Fill(kOccupancy); + + const XAxis xaxis{ + getCentrality(collision), + collision.multNTracksPV(), + (cfgTimeDependent) ? getTimeSinceStartOfFill(bc.timestamp(), *firstRunOfCurrentFill) : -1.0}; + + if (cfgTimeDependent && run == *firstRunOfCurrentFill && + firstRunOfCurrentFill != o2::analysis::gfwflowflucpp::firstRunsOfFill.end() - 1) + ++firstRunOfCurrentFill; + + if (cfgFillQA) + fillEventQA(collision, xaxis); + + registry.fill(HIST("eventQA/before/centrality"), xaxis.centrality); + registry.fill(HIST("eventQA/before/multiplicity"), xaxis.multiplicity); + + if (!eventSelected(collision, xaxis.multiplicity, xaxis.centrality, run)) + return; + + if (cfgFillQA) + fillEventQA(collision, xaxis); + + processCollision(collision, tracks, xaxis, run, 0); + + if (!collision.has_mcCollision()) + return; + + const auto genCollision = collision.template mcCollision_as(); + processGenCollision(genCollision, mcParticles, collision.mcCollisionId(), xaxis, run, 0); + } + PROCESS_SWITCH(FlowFlucGfwPp, processMC, "Process analysis for Monte-Carlo data", false); + +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{ + adaptAnalysisTask(cfgc), + }; +} From 28c53e5a6df7369f2e767bf9ed674e6a4cff6aff Mon Sep 17 00:00:00 2001 From: wenyaCern Date: Mon, 22 Jun 2026 16:52:12 +0200 Subject: [PATCH 2/4] Add cutoff of eta for q-vec and MC process --- PWGCF/Flow/flowFlucGfwPp.cxx | 1650 ---------------------------------- 1 file changed, 1650 deletions(-) delete mode 100644 PWGCF/Flow/flowFlucGfwPp.cxx diff --git a/PWGCF/Flow/flowFlucGfwPp.cxx b/PWGCF/Flow/flowFlucGfwPp.cxx deleted file mode 100644 index de048182729..00000000000 --- a/PWGCF/Flow/flowFlucGfwPp.cxx +++ /dev/null @@ -1,1650 +0,0 @@ -// Copyright 2019-2020 CERN and copyright holders of ALICE O2. -// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. -// All rights not expressly granted are reserved. -// -// This software is distributed under the terms of the GNU General Public -// License v3 (GPL Version 3), copied verbatim in the file "COPYING". -// -// In applying this license CERN does not waive the privileges and immunities -// granted to it by virtue of its status as an Intergovernmental Organization -// or submit itself to any jurisdiction. - -/// \file flowFlucGfwPp.cxx -/// \brief GFW task for Event Shape Engineering studies in pp collisions -/// \author Wenya Wu, TUM, wenya.wu@cern.ch, David Krylenkov, TUM - -#include "PWGCF/GenericFramework/Core/FlowContainer.h" -#include "PWGCF/GenericFramework/Core/GFW.h" -#include "PWGCF/GenericFramework/Core/GFWConfig.h" -#include "PWGCF/GenericFramework/Core/GFWWeights.h" - -#include "Common/CCDB/EventSelectionParams.h" -#include "Common/CCDB/TriggerAliases.h" -#include "Common/DataModel/Centrality.h" -#include "Common/DataModel/EventSelection.h" -#include "Common/DataModel/Multiplicity.h" -#include "Common/DataModel/TrackSelectionTables.h" - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include - -#include - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -using namespace o2; -using namespace o2::framework; -using namespace o2::analysis::genericframework; - -#define O2_DEFINE_CONFIGURABLE(NAME, TYPE, DEFAULT, HELP) Configurable NAME{#NAME, DEFAULT, HELP}; - -namespace o2::analysis::gfwflowflucpp -{ -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}; -float ptpoilow = 0.2, ptpoiup = 10.0; -float ptreflow = 0.2, ptrefup = 3.0; -float ptlow = 0.2, ptup = 10.0; -int etabins = 16; -float etalow = -0.8, etaup = 0.8; -int vtxZbins = 40; -float vtxZlow = -10.0, vtxZup = 10.0; -int phibins = 72; -float philow = 0.0; -float phiup = o2::constants::math::TwoPI; -int nchbins = 300; -float nchlow = 0; -float nchup = 3000; -std::vector centbinning(90); -int nBootstrap = 10; -GFWRegions regions; -GFWCorrConfigs configs; -std::vector multGlobalCorrCutPars; -std::vector multPVCorrCutPars; -std::vector multGlobalPVCorrCutPars; -std::vector multGlobalV0ACutPars; -std::vector multGlobalT0ACutPars; -std::vector firstRunsOfFill; -} // namespace o2::analysis::gfwflowflucpp - -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; - static constexpr int kRequireTwoTracksInBothEtaSides = 4; - static constexpr int kRequireTwoTracksInThreeEtaRegions = 8; - - static constexpr int kMinTracksForFourParticleCorrelation = 4; - static constexpr int kMinTracksPerEtaSideForGapCorrelation = 2; - static constexpr int kMinTracksPerEtaRegionForThreeSubevents = 2; - - static constexpr int EllipticQVectorHarmonic = 2; - static constexpr int TriangularQVectorHarmonic = 3; - - O2_DEFINE_CONFIGURABLE(cfgNbootstrap, int, 10, "Number of subsamples") - 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)") - O2_DEFINE_CONFIGURABLE(cfgTimeDependent, bool, false, "Fill output as function of time (for contamination studies)") - O2_DEFINE_CONFIGURABLE(cfgFirstRunsOfFill, std::vector, {}, "First runs of a fill for time dependent analysis") - O2_DEFINE_CONFIGURABLE(cfgFillQA, bool, false, "Fill QA histograms") - O2_DEFINE_CONFIGURABLE(cfgUseAdditionalEventCut, bool, false, "Use additional event cut on mult correlations") - O2_DEFINE_CONFIGURABLE(cfgEfficiency, std::string, "", "CCDB path to efficiency object") - O2_DEFINE_CONFIGURABLE(cfgAcceptance, std::string, "", "CCDB path to acceptance object") - O2_DEFINE_CONFIGURABLE(cfgDCAxyNSigma, float, 7, "Cut on number of sigma deviations from expected DCA in the transverse direction"); - O2_DEFINE_CONFIGURABLE(cfgDCAxy, std::string, "(0.0026+0.005/(x^1.01))", "Functional form of pt-dependent DCAxy cut"); - O2_DEFINE_CONFIGURABLE(cfgDCAz, float, 2, "Cut on DCA in the longitudinal direction (cm)"); - O2_DEFINE_CONFIGURABLE(cfgNTPCCls, float, 50, "Cut on number of TPC clusters found"); - O2_DEFINE_CONFIGURABLE(cfgNTPCXrows, float, 70, "Cut on number of TPC crossed rows"); - O2_DEFINE_CONFIGURABLE(cfgMinNITSCls, float, 5, "Cut on minimum number of ITS clusters found"); - O2_DEFINE_CONFIGURABLE(cfgChi2PrITSCls, float, 36, "Cut on chi^2 per ITS clusters found"); - O2_DEFINE_CONFIGURABLE(cfgChi2PrTPCCls, float, 2.5, "Cut on chi^2 per TPC clusters found"); - O2_DEFINE_CONFIGURABLE(cfgPtmin, float, 0.2, "minimum pt (GeV/c)"); - O2_DEFINE_CONFIGURABLE(cfgPtmax, float, 10, "maximum pt (GeV/c)"); - O2_DEFINE_CONFIGURABLE(cfgEta, float, 0.8, "eta cut"); - O2_DEFINE_CONFIGURABLE(cfgVtxZ, float, 10, "vertex cut (cm)"); - O2_DEFINE_CONFIGURABLE(cfgNoSameBunchPileupCut, bool, false, "kNoSameBunchPileupCut"); - O2_DEFINE_CONFIGURABLE(cfgIsGoodZvtxFT0vsPV, bool, false, "kIsGoodZvtxFT0vsPV"); - O2_DEFINE_CONFIGURABLE(cfgIsGoodITSLayersAll, bool, false, "kIsGoodITSLayersAll"); - O2_DEFINE_CONFIGURABLE(cfgNoCollInTimeRangeStandard, bool, false, "kNoCollInTimeRangeStandard"); - O2_DEFINE_CONFIGURABLE(cfgDoOccupancySel, bool, false, "Bool for event selection on detector occupancy"); - O2_DEFINE_CONFIGURABLE(cfgOccupancySelection, int, 5000, "Max occupancy selection, -999 to disable"); - O2_DEFINE_CONFIGURABLE(cfgMultCut, bool, false, "Use additional event cut on mult correlations"); - O2_DEFINE_CONFIGURABLE(cfgTVXinTRD, bool, false, "Use kTVXinTRD (reject TRD triggered events)"); - O2_DEFINE_CONFIGURABLE(cfgIsVertexITSTPC, bool, true, "Selects collisions with at least one ITS-TPC track"); - O2_DEFINE_CONFIGURABLE(cfgMagField, float, 99999, "Configurable magnetic field; default CCDB will be queried"); - O2_DEFINE_CONFIGURABLE(cfgFixedMultMin, int, 1, "Minimum for fixed nch range"); - O2_DEFINE_CONFIGURABLE(cfgUseMultiplicityFlowWeights, bool, true, "Enable or disable the use of multiplicity-based event weighting"); - O2_DEFINE_CONFIGURABLE(cfgFixedMultMax, int, 3000, "Maximum for fixed nch range"); - O2_DEFINE_CONFIGURABLE(cfgConsistentEventFlag, int, 0, "Flag to select consistent events - 0: off, 1: v2{2} gap calculable, 2: v2{4} full calculable, 4: v2{4} gap calculable, 8: v2{4} 3sub calculable"); - Configurable> cfgMultGlobalCutPars{"cfgMultGlobalCutPars", std::vector{2272.16, -76.6932, 1.01204, -0.00631545, 1.59868e-05, 136.336, -4.97006, 0.121199, -0.0015921, 7.66197e-06}, "Global vs FT0C multiplicity cut parameter values"}; - Configurable> cfgMultPVCutPars{"cfgMultPVCutPars", std::vector{3074.43, -106.192, 1.46176, -0.00968364, 2.61923e-05, 182.128, -7.43492, 0.193901, -0.00256715, 1.22594e-05}, "PV vs FT0C multiplicity cut parameter values"}; - Configurable> cfgMultGlobalPVCutPars{"cfgMultGlobalPVCutPars", std::vector{-0.223013, 0.715849, 0.664242, 0.0829653, -0.000503733, 1.21185e-06}, "Global vs PV multiplicity cut parameter values"}; - O2_DEFINE_CONFIGURABLE(cfgMultCorrHighCutFunction, std::string, "[0] + [1]*x + [2]*x*x + [3]*x*x*x + [4]*x*x*x*x + 3.*([5] + [6]*x + [7]*x*x + [8]*x*x*x + [9]*x*x*x*x)", "Functional for multiplicity correlation cut"); - O2_DEFINE_CONFIGURABLE(cfgMultCorrLowCutFunction, std::string, "[0] + [1]*x + [2]*x*x + [3]*x*x*x + [4]*x*x*x*x - 3.*([5] + [6]*x + [7]*x*x + [8]*x*x*x + [9]*x*x*x*x)", "Functional for multiplicity correlation cut"); - O2_DEFINE_CONFIGURABLE(cfgMultGlobalPVCorrCutFunction, std::string, "[0] + [1]*x + 3*([2] + [3]*x + [4]*x*x + [5]*x*x*x)", "Functional for global vs pv multiplicity correlation cut"); - O2_DEFINE_CONFIGURABLE(cfgNumQnBins, int, 10, "Number of qn bins"); - 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 qn selection; otherwise use +eta half"); - O2_DEFINE_CONFIGURABLE(cfgQnSelectionHarmonic, int, 2, "Harmonic n used to build the reduced q_n vector for event shape selection, use 2 for q2 and 3 for q3"); - O2_DEFINE_CONFIGURABLE(cfgQnHistMax, float, 6., "Upper range for q_n calibration histograms"); - O2_DEFINE_CONFIGURABLE(cfgQnTrkAbsEtaMax, float, 0.5, "Upper range for abs eta of tracks contributing to q_n"); - O2_DEFINE_CONFIGURABLE(cfgBypassQnSelection, bool, false, "Bypass q_n event shape selection and fill one integral q-bin"); - 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 { - O2_DEFINE_CONFIGURABLE(cfgMultGlobalASideCorrCutFunction, std::string, "[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x + [10]*([5]+[6]*x+[7]*x*x+[8]*x*x*x+[9]*x*x*x*x)", "Functional for global vs V0A multiplicity low correlation cut"); - Configurable> cfgMultGlobalV0ACutPars{"cfgMultGlobalV0ACutPars", std::vector{567.785, 172.715, 0.77888, -0.00693466, 1.40564e-05, 679.853, 66.8068, -0.444332, 0.00115002, -4.92064e-07}, "Global vs FV0A multiplicity cut parameter values"}; - Configurable> cfgMultGlobalT0ACutPars{"cfgMultGlobalT0ACutPars", std::vector{241.618, 61.8402, 0.348049, -0.00306078, 6.20357e-06, 315.235, 29.1491, -0.188639, 0.00044528, -9.08912e-08}, "Global vs FT0A multiplicity cut parameter values"}; - O2_DEFINE_CONFIGURABLE(cfgGlobalV0ALowSigma, float, -3, "Number of sigma deviations below expected value in global vs V0A correlation"); - O2_DEFINE_CONFIGURABLE(cfgGlobalV0AHighSigma, float, 4, "Number of sigma deviations above expected value in global vs V0A correlation"); - O2_DEFINE_CONFIGURABLE(cfgGlobalT0ALowSigma, float, -3., "Number of sigma deviations below expected value in global vs T0A correlation"); - O2_DEFINE_CONFIGURABLE(cfgGlobalT0AHighSigma, float, 4, "Number of sigma deviations above expected value in global vs T0A correlation"); - } cfgGlobalAsideCorrCuts; - - Configurable cfgGfwBinning{"cfgGfwBinning", - {40, 16, 72, 300, 0, 3000, 0.2, 10.0, 0.2, 3.0, {0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1, 1.25, 1.5, 1.75, 2, 2.25, 2.5, 2.75, 3, 3.25, 3.5, 3.75, 4, 4.5, 5, 5.5, 6, 7, 8, 9, 10}, {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90}}, - "Configuration for binning"}; - - Configurable cfgRegions{"cfgRegions", - { - {"refN", "refP", "refFull"}, - {-0.8, 0.4, -0.8}, - {-0.4, 0.8, 0.8}, - {0, 0, 0}, // pT bins - {1, 1, 1} // bitmask - }, - "Configurations for GFW regions"}; - - Configurable cfgCorrConfig{"cfgCorrConfig", - {{"refN {2 -2}", - "refN {2 2 -2 -2}", - "refN {2 2 2 -2 -2 -2}", - "refN {2 2 2 2 -2 -2 -2 -2}", - "refP {2 -2}", - "refP {2 2 -2 -2}", - "refP {2 2 2 -2 -2 -2}", - "refP {2 2 2 2 -2 -2 -2 -2}", - "refN {2} refP {-2}", - "refN {2 2} refP {-2 -2}", - "refFull {2 -2}", - "refFull {2 2 -2 -2}", - "refFull {2 2 2 -2 -2 -2}", - "refFull {2 2 2 2 -2 -2 -2 -2}"}, - {"ChNeg22", - "ChNeg24", - "ChNeg26", - "ChNeg28", - "ChPos22", - "ChPos24", - "ChPos26", - "ChPos28", - "ChGap22", - "ChGap24", - "ChFull22", - "ChFull24", - "ChFull26", - "ChFull28"}, - {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, - {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}}, - "Configurations for pp ESE v2 cumulants"}; - - // Connect to ccdb - Service ccdb; - - struct Config { - TH1D* mEfficiency = nullptr; - GFWWeights* mAcceptance; - bool correctionsLoaded = false; - } cfg; - - // Define output - OutputObj fFC{FlowContainer("FlowContainer")}; - OutputObj fFCgen{FlowContainer("FlowContainer_gen")}; - HistogramRegistry registry{"registry"}; - - // QA outputs - std::map>> th1sList; - std::map>> tpfsList; - std::map>> th3sList; - enum OutputTH1Names { - hPhi = 0, - hEta, - hVtxZ, - hMult, - hCent, - hEventSel, - kCount_TH1Names - }; - enum OutputTProfileNames { - pfCorr22 = 0, - kCount_TProfileNames - }; - // NUA outputs - enum OutputTH3Names { - hNUAref = 0, - kCount_TH3Names - }; - enum CentEstimators { - kCentFT0C = 0, - kCentFT0CVariant1, - kCentFT0M, - kCentFV0A, - kCentNTPV, - kCentNGlobal, - kCentMFT - }; - enum EventSelFlags { - kFilteredEvent = 1, - kSel8, - kOccupancy, - kTVXTRD, - kNoSamebunchPU, - kZVtxFT0PV, - kNoCollTRStd, - kVtxITSTPC, - kGoodITSLayers, - kMultCuts, - kTrackCent - }; - - // Define global variables - // Generic Framework - GFW* fGFW = new GFW(); - std::vector corrconfigs; - - TRandom3* fRndm = new TRandom3(0); - TAxis* fSecondAxis; - int lastRun = -1; - std::vector::iterator firstRunOfCurrentFill; - std::vector runNumbers; - - std::string getShapeSel() const - { - return "ese"; - } - - int getNQnOutputBins() - { - return static_cast(cfgBypassQnSelection) ? 1 : static_cast(cfgNumQnBins); - } - - // region indices for consistency flag - int posRegionIndex = -1; - int negRegionIndex = -1; - int fullRegionIndex = -1; - int midRegionIndex = -1; - - // Event selection cuts - Alex - TF1* fMultPVCutLow = nullptr; - TF1* fMultPVCutHigh = nullptr; - TF1* fMultCutLow = nullptr; - TF1* fMultCutHigh = nullptr; - TF1* fMultPVGlobalCutHigh = nullptr; - TF1* fMultGlobalV0ACutLow = nullptr; - TF1* fMultGlobalV0ACutHigh = nullptr; - TF1* fMultGlobalT0ACutLow = nullptr; - TF1* fMultGlobalT0ACutHigh = nullptr; - - TF1* fPtDepDCAxy = nullptr; - - o2::framework::expressions::Filter collisionFilter = nabs(aod::collision::posZ) < cfgVtxZ; - o2::framework::expressions::Filter trackFilter = nabs(aod::track::eta) < cfgEta && aod::track::pt > cfgPtmin&& aod::track::pt < cfgPtmax && (aod::track::itsChi2NCl < cfgChi2PrITSCls) && (aod::track::tpcChi2NCl < cfgChi2PrTPCCls) && nabs(aod::track::dcaZ) < cfgDCAz; - - Preslice perCollision = aod::track::collisionId; - 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 GFWTracksMC = 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()); - o2::analysis::gfwflowflucpp::regions.SetpTDifs(cfgRegions->GetpTDifs()); - o2::analysis::gfwflowflucpp::regions.SetBitmasks(cfgRegions->GetBitmasks()); - o2::analysis::gfwflowflucpp::configs.SetCorrs(cfgCorrConfig->GetCorrs()); - o2::analysis::gfwflowflucpp::configs.SetHeads(cfgCorrConfig->GetHeads()); - o2::analysis::gfwflowflucpp::configs.SetpTDifs(cfgCorrConfig->GetpTDifs()); - o2::analysis::gfwflowflucpp::configs.SetpTCorrMasks(cfgCorrConfig->GetpTCorrMasks()); - - o2::analysis::gfwflowflucpp::regions.Print(); - o2::analysis::gfwflowflucpp::configs.Print(); - - o2::analysis::gfwflowflucpp::ptbinning = cfgGfwBinning->GetPtBinning(); - // o2::analysis::gfwflowflucpp::ptpoilow = cfgGfwBinning->GetPtPOImin(); - // o2::analysis::gfwflowflucpp::ptpoiup = cfgGfwBinning->GetPtPOImax(); - o2::analysis::gfwflowflucpp::ptreflow = cfgGfwBinning->GetPtRefMin(); - o2::analysis::gfwflowflucpp::ptrefup = cfgGfwBinning->GetPtRefMax(); - o2::analysis::gfwflowflucpp::ptlow = cfgPtmin; - o2::analysis::gfwflowflucpp::ptup = cfgPtmax; - o2::analysis::gfwflowflucpp::etabins = cfgGfwBinning->GetEtaBins(); - o2::analysis::gfwflowflucpp::vtxZbins = cfgGfwBinning->GetVtxZbins(); - o2::analysis::gfwflowflucpp::phibins = cfgGfwBinning->GetPhiBins(); - o2::analysis::gfwflowflucpp::philow = 0.0f; - o2::analysis::gfwflowflucpp::phiup = o2::constants::math::TwoPI; - o2::analysis::gfwflowflucpp::nchbins = cfgGfwBinning->GetNchBins(); - o2::analysis::gfwflowflucpp::nchlow = cfgGfwBinning->GetNchMin(); - o2::analysis::gfwflowflucpp::nchup = cfgGfwBinning->GetNchMax(); - o2::analysis::gfwflowflucpp::centbinning = cfgGfwBinning->GetCentBinning(); - cfgGfwBinning->Print(); - - o2::analysis::gfwflowflucpp::multGlobalCorrCutPars = cfgMultGlobalCutPars; - o2::analysis::gfwflowflucpp::multPVCorrCutPars = cfgMultPVCutPars; - o2::analysis::gfwflowflucpp::multGlobalPVCorrCutPars = cfgMultGlobalPVCutPars; - o2::analysis::gfwflowflucpp::multGlobalV0ACutPars = cfgGlobalAsideCorrCuts.cfgMultGlobalV0ACutPars; - o2::analysis::gfwflowflucpp::multGlobalT0ACutPars = cfgGlobalAsideCorrCuts.cfgMultGlobalT0ACutPars; - o2::analysis::gfwflowflucpp::firstRunsOfFill = cfgFirstRunsOfFill; - if (cfgTimeDependent && !std::is_sorted(o2::analysis::gfwflowflucpp::firstRunsOfFill.begin(), o2::analysis::gfwflowflucpp::firstRunsOfFill.end())) { - std::sort(o2::analysis::gfwflowflucpp::firstRunsOfFill.begin(), o2::analysis::gfwflowflucpp::firstRunsOfFill.end()); - } - firstRunOfCurrentFill = o2::analysis::gfwflowflucpp::firstRunsOfFill.begin(); - - AxisSpec phiAxis = {o2::analysis::gfwflowflucpp::phibins, o2::analysis::gfwflowflucpp::philow, o2::analysis::gfwflowflucpp::phiup, "#phi"}; - AxisSpec etaAxis = {o2::analysis::gfwflowflucpp::etabins, -cfgEta, cfgEta, "#eta"}; - AxisSpec vtxAxis = {o2::analysis::gfwflowflucpp::vtxZbins, -cfgVtxZ, cfgVtxZ, "Vtx_{z} (cm)"}; - AxisSpec ptAxis = {o2::analysis::gfwflowflucpp::ptbinning, "#it{p}_{T} GeV/#it{c}"}; - - std::string sCentralityEstimator; - std::map centEstimatorMap = { - {kCentFT0C, "FT0C"}, - {kCentFT0CVariant1, "FT0C variant 1"}, - {kCentFT0M, "FT0M"}, - {kCentFV0A, "FV0A"}, - {kCentNTPV, "NTPV"}, - {kCentNGlobal, "NGlobals"}, - {kCentMFT, "MFT"}}; - sCentralityEstimator = centEstimatorMap.at(cfgCentEstimator); - sCentralityEstimator += " centrality (%)"; - AxisSpec centAxis = {o2::analysis::gfwflowflucpp::centbinning, sCentralityEstimator.c_str()}; - std::vector nchbinning; - int nchskip = (o2::analysis::gfwflowflucpp::nchup - o2::analysis::gfwflowflucpp::nchlow) / o2::analysis::gfwflowflucpp::nchbins; - for (int i = 0; i <= o2::analysis::gfwflowflucpp::nchbins; ++i) { - nchbinning.push_back(nchskip * i + o2::analysis::gfwflowflucpp::nchlow + 0.5); - } - AxisSpec nchAxis = {nchbinning, "N_{ch}"}; - std::vector bbinning(201); - std::generate(bbinning.begin(), bbinning.end(), [n = -0.1, step = 0.1]() mutable { - n += step; - return n; - }); - AxisSpec bAxis = {bbinning, "#it{b}"}; - AxisSpec t0cAxis = {1000, 0, 10000, "N_{ch} (T0C)"}; - AxisSpec t0aAxis = {300, 0, 30000, "N_{ch} (T0A)"}; - AxisSpec v0aAxis = {800, 0, 80000, "N_{ch} (V0A)"}; - AxisSpec multpvAxis = {600, 0, 600, "N_{ch} (PV)"}; - AxisSpec dcaZAXis = {200, -2, 2, "DCA_{z} (cm)"}; - AxisSpec dcaXYAXis = {200, -0.5, 0.5, "DCA_{xy} (cm)"}; - std::vector timebinning(289); - std::generate(timebinning.begin(), timebinning.end(), [n = -24 / 288., step = 24 / 288.]() mutable { - n += step; - return n; - }); - AxisSpec timeAxis = {timebinning, "time (hrs)"}; - - AxisSpec multAxis = (cfgTimeDependent) ? timeAxis : (cfgUseNch) ? nchAxis - : centAxis; - - ccdb->setURL("http://alice-ccdb.cern.ch"); - ccdb->setCaching(true); - ccdb->setLocalObjectValidityChecking(); - - int64_t now = std::chrono::duration_cast(std::chrono::system_clock::now().time_since_epoch()).count(); - ccdb->setCreatedNotAfter(now); - - 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 (static_cast(cfgBypassQnSelection)) { - LOGF(info, "Bypassing q_n event shape selection, all accepted events will be filled into the integral bin ese_0"); - if (static_cast(cfgNumQnBins) != 1) { - LOGF(warning, "cfgBypassQnSelection is on, cfgNumQnBins=%d will be ignored and only one output q-bin will be made", static_cast(cfgNumQnBins)); - } - } else { - LOGF(info, "Event-shape selector uses q_%d from the %s eta half", - static_cast(cfgQnSelectionHarmonic), - static_cast(cfgUseNegativeEtaHalfForq2) ? "negative" : "positive"); - } - - 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) { - const int qnHarmonic = static_cast(cfgQnSelectionHarmonic); - const double qnHistMax = static_cast(cfgQnHistMax); - - if (qnHarmonic == EllipticQVectorHarmonic) { - 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, qnHistMax}}); - registry.add("mq2/h2_cent_q2_etaneg", ";Centrality;#it{q}_{2}^{#eta neg};", HistType::kTH2D, {{100, 0, 100}, {600, 0, qnHistMax}}); - registry.add("mq2/h2_mult_q2_etapos", ";Multiplicity;#it{q}_{2}^{#eta pos};", HistType::kTH2D, {{150, 0, 150}, {600, 0, qnHistMax}}); - registry.add("mq2/h2_mult_q2_etaneg", ";Multiplicity;#it{q}_{2}^{#eta neg};", HistType::kTH2D, {{150, 0, 150}, {600, 0, qnHistMax}}); - } - - if (qnHarmonic == TriangularQVectorHarmonic) { - registry.add("mq3/eventcounter", "", HistType::kTH1F, {{10, 0, 10}}); - registry.add("mq3/h2_cent_q3_etapos", ";Centrality;#it{q}_{3}^{#eta pos};", HistType::kTH2D, {{100, 0, 100}, {600, 0, qnHistMax}}); - registry.add("mq3/h2_cent_q3_etaneg", ";Centrality;#it{q}_{3}^{#eta neg};", HistType::kTH2D, {{100, 0, 100}, {600, 0, qnHistMax}}); - registry.add("mq3/h2_mult_q3_etapos", ";Multiplicity;#it{q}_{3}^{#eta pos};", HistType::kTH2D, {{150, 0, 150}, {600, 0, qnHistMax}}); - registry.add("mq3/h2_mult_q3_etaneg", ";Multiplicity;#it{q}_{3}^{#eta neg};", HistType::kTH2D, {{150, 0, 150}, {600, 0, qnHistMax}}); - } - } - - if (doprocessData || doprocessMC) { - registry.add("trackQA/before/phi_eta_vtxZ", "", {HistType::kTH3D, {phiAxis, etaAxis, vtxAxis}}); - registry.add("trackQA/before/pt_dcaXY_dcaZ", "", {HistType::kTH3D, {ptAxis, dcaXYAXis, dcaZAXis}}); - registry.add("trackQA/before/nch_pt", "#it{p}_{T} vs multiplicity; N_{ch}; #it{p}_{T}", {HistType::kTH2D, {nchAxis, ptAxis}}); - registry.add("trackQA/before/chi2prTPCcls", "#chi^{2}/cluster for the TPC track segment; #chi^{2}/TPC cluster", {HistType::kTH1D, {{100, 0., 5.}}}); - registry.add("trackQA/before/chi2prITScls", "#chi^{2}/cluster for the ITS track; #chi^{2}/ITS cluster", {HistType::kTH1D, {{100, 0., 50.}}}); - registry.add("trackQA/before/nTPCClusters", "Number of found TPC clusters; TPC N_{cls}; Counts", {HistType::kTH1D, {{100, 40, 180}}}); - registry.add("trackQA/before/nITSClusters", "Number of found ITS clusters; ITS N_{cls}; Counts", {HistType::kTH1D, {{100, 0, 20}}}); - registry.add("trackQA/before/nTPCCrossedRows", "Number of crossed TPC Rows; TPC X-rows; Counts", {HistType::kTH1D, {{100, 40, 180}}}); - - registry.addClone("trackQA/before/", "trackQA/after/"); - registry.add("trackQA/after/pt_ref", "", {HistType::kTH1D, {{100, o2::analysis::gfwflowflucpp::ptreflow, o2::analysis::gfwflowflucpp::ptrefup}}}); - registry.add("trackQA/after/pt_poi", "", {HistType::kTH1D, {{100, o2::analysis::gfwflowflucpp::ptpoilow, o2::analysis::gfwflowflucpp::ptpoiup}}}); - - registry.add("eventQA/before/multiplicity", "", {HistType::kTH1D, {nchAxis}}); - if (cfgTimeDependent) { - registry.add("eventQA/before/multiplicity_time", "Multiplicity vs time; time (hrs); N_{ch}", {HistType::kTH2D, {timeAxis, nchAxis}}); - registry.add("eventQA/before/multT0C_time", "T0C Multiplicity vs time; time (hrs); N_{ch} (T0C)", {HistType::kTH2D, {timeAxis, t0cAxis}}); - registry.add("eventQA/before/multT0A_time", "T0A Multiplicity vs time; time (hrs); N_{ch} (T0A)", {HistType::kTH2D, {timeAxis, t0aAxis}}); - registry.add("eventQA/before/multV0A_time", "V0A Multiplicity vs time; time (hrs); N_{ch} (V0A)", {HistType::kTH2D, {timeAxis, v0aAxis}}); - registry.add("eventQA/before/multPV_time", "PV Multiplicity vs time; time (hrs); N_{ch} (PV)", {HistType::kTH2D, {timeAxis, multpvAxis}}); - } - registry.add("eventQA/before/globalTracks_PVTracks", "", {HistType::kTH2D, {multpvAxis, nchAxis}}); - registry.add("eventQA/before/globalTracks_multT0A", "", {HistType::kTH2D, {t0aAxis, nchAxis}}); - registry.add("eventQA/before/globalTracks_multV0A", "", {HistType::kTH2D, {v0aAxis, nchAxis}}); - registry.add("eventQA/before/multV0A_multT0A", "", {HistType::kTH2D, {t0aAxis, v0aAxis}}); - - registry.add("eventQA/before/centrality", "", {HistType::kTH1D, {centAxis}}); - registry.add("eventQA/before/globalTracks_centT0C", "", {HistType::kTH2D, {centAxis, nchAxis}}); - registry.add("eventQA/before/PVTracks_centT0C", "", {HistType::kTH2D, {centAxis, multpvAxis}}); - registry.add("eventQA/before/multT0C_centT0C", "", {HistType::kTH2D, {centAxis, t0cAxis}}); - - registry.add("eventQA/before/centT0M_centT0C", "", {HistType::kTH2D, {centAxis, centAxis}}); - registry.add("eventQA/before/centV0A_centT0C", "", {HistType::kTH2D, {centAxis, centAxis}}); - registry.add("eventQA/before/centGlobal_centT0C", "", {HistType::kTH2D, {centAxis, centAxis}}); - registry.add("eventQA/before/centNTPV_centT0C", "", {HistType::kTH2D, {centAxis, centAxis}}); - registry.add("eventQA/before/centMFT_centT0C", "", {HistType::kTH2D, {centAxis, centAxis}}); - - if (cfgIsMC || doprocessMC) { - registry.add("MCGen/trackQA/phi_eta_vtxZ", "", {HistType::kTH3D, {phiAxis, etaAxis, vtxAxis}}); - registry.add("MCGen/trackQA/nch_pt", "#it{p}_{T} vs multiplicity; N_{ch}; #it{p}_{T}", {HistType::kTH2D, {nchAxis, ptAxis}}); - registry.add("MCGen/trackQA/pt_ref", "", {HistType::kTH1D, {{100, o2::analysis::gfwflowflucpp::ptreflow, o2::analysis::gfwflowflucpp::ptrefup}}}); - registry.add("MCGen/trackQA/pt_poi", "", {HistType::kTH1D, {{100, o2::analysis::gfwflowflucpp::ptpoilow, o2::analysis::gfwflowflucpp::ptpoiup}}}); - } - - registry.addClone("eventQA/before/", "eventQA/after/"); - registry.add("eventQA/eventSel", "Number of Events;; Counts", {HistType::kTH1D, {{11, 0.5, 11.5}}}); - registry.get(HIST("eventQA/eventSel"))->GetXaxis()->SetBinLabel(kFilteredEvent, "Filtered event"); - registry.get(HIST("eventQA/eventSel"))->GetXaxis()->SetBinLabel(kSel8, "sel8"); - registry.get(HIST("eventQA/eventSel"))->GetXaxis()->SetBinLabel(kOccupancy, "occupancy"); - registry.get(HIST("eventQA/eventSel"))->GetXaxis()->SetBinLabel(kTVXTRD, "kTVXinTRD"); - registry.get(HIST("eventQA/eventSel"))->GetXaxis()->SetBinLabel(kNoSamebunchPU, "kNoSameBunchPileup"); - registry.get(HIST("eventQA/eventSel"))->GetXaxis()->SetBinLabel(kZVtxFT0PV, "kIsGoodZvtxFT0vsPV"); - registry.get(HIST("eventQA/eventSel"))->GetXaxis()->SetBinLabel(kNoCollTRStd, "kNoCollInTimeRangeStandard"); - registry.get(HIST("eventQA/eventSel"))->GetXaxis()->SetBinLabel(kVtxITSTPC, "kIsVertexITSTPC"); - registry.get(HIST("eventQA/eventSel"))->GetXaxis()->SetBinLabel(kGoodITSLayers, "kIsGoodITSLayersAll"); - registry.get(HIST("eventQA/eventSel"))->GetXaxis()->SetBinLabel(kMultCuts, "after Mult cuts"); - registry.get(HIST("eventQA/eventSel"))->GetXaxis()->SetBinLabel(kTrackCent, "has track + within cent"); - if (!cfgRunByRun && cfgFillWeights) { - registry.add("phi_eta_vtxz_ref", "", {HistType::kTH3D, {phiAxis, etaAxis, vtxAxis}}); - } - } - - if (o2::analysis::gfwflowflucpp::regions.GetSize() < 0) - LOGF(error, "Configuration contains vectors of different size - check the GFWRegions configurable"); - for (auto i(0); i < o2::analysis::gfwflowflucpp::regions.GetSize(); ++i) { - fGFW->AddRegion(o2::analysis::gfwflowflucpp::regions.GetNames()[i], o2::analysis::gfwflowflucpp::regions.GetEtaMin()[i], o2::analysis::gfwflowflucpp::regions.GetEtaMax()[i], (o2::analysis::gfwflowflucpp::regions.GetpTDifs()[i]) ? ptbins + 1 : 1, o2::analysis::gfwflowflucpp::regions.GetBitmasks()[i]); - } - for (auto i = 0; i < o2::analysis::gfwflowflucpp::configs.GetSize(); ++i) { - corrconfigs.push_back(fGFW->GetCorrelatorConfig(o2::analysis::gfwflowflucpp::configs.GetCorrs()[i], o2::analysis::gfwflowflucpp::configs.GetHeads()[i], o2::analysis::gfwflowflucpp::configs.GetpTDifs()[i])); - } - if (corrconfigs.empty()) - LOGF(error, "Configuration contains vectors of different size - check the GFWCorrConfig configurable"); - fGFW->CreateRegions(); - TObjArray* oba = new TObjArray(); - addConfigObjectsToObjArray(oba, corrconfigs); - fFC->SetName("FlowContainer"); - fFC->SetXAxis(fSecondAxis); - fFC->Initialize(oba, multAxis, cfgNbootstrap); - - fFCgen->SetName("FlowContainer_gen"); - fFCgen->SetXAxis(fSecondAxis); - fFCgen->Initialize(oba, multAxis, cfgNbootstrap); - - delete oba; - - fPtDepDCAxy = new TF1("ptDepDCAxy", Form("[0]*%s", cfgDCAxy->c_str()), 0.001, 100); - fPtDepDCAxy->SetParameter(0, cfgDCAxyNSigma); - LOGF(info, "DCAxy pt-dependence function: %s", Form("[0]*%s", cfgDCAxy->c_str())); - if (cfgUseAdditionalEventCut) { - fMultPVCutLow = new TF1("fMultPVCutLow", cfgMultCorrLowCutFunction->c_str(), 0, 100); - fMultPVCutLow->SetParameters(&(o2::analysis::gfwflowflucpp::multPVCorrCutPars[0])); - fMultPVCutHigh = new TF1("fMultPVCutHigh", cfgMultCorrHighCutFunction->c_str(), 0, 100); - fMultPVCutHigh->SetParameters(&(o2::analysis::gfwflowflucpp::multPVCorrCutPars[0])); - fMultCutLow = new TF1("fMultCutLow", cfgMultCorrLowCutFunction->c_str(), 0, 100); - fMultCutLow->SetParameters(&(o2::analysis::gfwflowflucpp::multGlobalCorrCutPars[0])); - fMultCutHigh = new TF1("fMultCutHigh", cfgMultCorrHighCutFunction->c_str(), 0, 100); - fMultCutHigh->SetParameters(&(o2::analysis::gfwflowflucpp::multGlobalCorrCutPars[0])); - fMultPVGlobalCutHigh = new TF1("fMultPVGlobalCutHigh", cfgMultGlobalPVCorrCutFunction->c_str(), 0, nchbinning.back()); - fMultPVGlobalCutHigh->SetParameters(&(o2::analysis::gfwflowflucpp::multGlobalPVCorrCutPars[0])); - - LOGF(info, "Global V0A function: %s in range 0-%g", cfgGlobalAsideCorrCuts.cfgMultGlobalASideCorrCutFunction->c_str(), v0aAxis.binEdges.back()); - fMultGlobalV0ACutLow = new TF1("fMultGlobalV0ACutLow", cfgGlobalAsideCorrCuts.cfgMultGlobalASideCorrCutFunction->c_str(), 0, v0aAxis.binEdges.back()); - for (std::size_t i = 0; i < o2::analysis::gfwflowflucpp::multGlobalV0ACutPars.size(); ++i) - fMultGlobalV0ACutLow->SetParameter(i, o2::analysis::gfwflowflucpp::multGlobalV0ACutPars[i]); - fMultGlobalV0ACutLow->SetParameter(o2::analysis::gfwflowflucpp::multGlobalV0ACutPars.size(), cfgGlobalAsideCorrCuts.cfgGlobalV0ALowSigma); - for (int i = 0; i < fMultGlobalV0ACutLow->GetNpar(); ++i) - LOGF(info, "fMultGlobalV0ACutLow par %d = %g", i, fMultGlobalV0ACutLow->GetParameter(i)); - - fMultGlobalV0ACutHigh = new TF1("fMultGlobalV0ACutHigh", cfgGlobalAsideCorrCuts.cfgMultGlobalASideCorrCutFunction->c_str(), 0, v0aAxis.binEdges.back()); - for (std::size_t i = 0; i < o2::analysis::gfwflowflucpp::multGlobalV0ACutPars.size(); ++i) - fMultGlobalV0ACutHigh->SetParameter(i, o2::analysis::gfwflowflucpp::multGlobalV0ACutPars[i]); - fMultGlobalV0ACutHigh->SetParameter(o2::analysis::gfwflowflucpp::multGlobalV0ACutPars.size(), cfgGlobalAsideCorrCuts.cfgGlobalV0AHighSigma); - for (int i = 0; i < fMultGlobalV0ACutHigh->GetNpar(); ++i) - LOGF(info, "fMultGlobalV0ACutHigh par %d = %g", i, fMultGlobalV0ACutHigh->GetParameter(i)); - - LOGF(info, "Global T0A function: %s", cfgGlobalAsideCorrCuts.cfgMultGlobalASideCorrCutFunction->c_str()); - fMultGlobalT0ACutLow = new TF1("fMultGlobalT0ACutLow", cfgGlobalAsideCorrCuts.cfgMultGlobalASideCorrCutFunction->c_str(), 0, t0aAxis.binEdges.back()); - for (std::size_t i = 0; i < o2::analysis::gfwflowflucpp::multGlobalT0ACutPars.size(); ++i) - fMultGlobalT0ACutLow->SetParameter(i, o2::analysis::gfwflowflucpp::multGlobalT0ACutPars[i]); - fMultGlobalT0ACutLow->SetParameter(o2::analysis::gfwflowflucpp::multGlobalT0ACutPars.size(), cfgGlobalAsideCorrCuts.cfgGlobalT0ALowSigma); - for (int i = 0; i < fMultGlobalT0ACutLow->GetNpar(); ++i) - LOGF(info, "fMultGlobalT0ACutLow par %d = %g", i, fMultGlobalT0ACutLow->GetParameter(i)); - - fMultGlobalT0ACutHigh = new TF1("fMultGlobalT0ACutHigh", cfgGlobalAsideCorrCuts.cfgMultGlobalASideCorrCutFunction->c_str(), 0, t0aAxis.binEdges.back()); - for (std::size_t i = 0; i < o2::analysis::gfwflowflucpp::multGlobalT0ACutPars.size(); ++i) - fMultGlobalT0ACutHigh->SetParameter(i, o2::analysis::gfwflowflucpp::multGlobalT0ACutPars[i]); - fMultGlobalT0ACutHigh->SetParameter(o2::analysis::gfwflowflucpp::multGlobalT0ACutPars.size(), cfgGlobalAsideCorrCuts.cfgGlobalT0AHighSigma); - for (int i = 0; i < fMultGlobalT0ACutHigh->GetNpar(); ++i) - LOGF(info, "fMultGlobalT0ACutHigh par %d = %g", i, fMultGlobalT0ACutHigh->GetParameter(i)); - } - - if (cfgConsistentEventFlag) { - const auto& names = cfgRegions->GetNames(); - auto findRegionIndex = [&](const std::string& name) { - auto it = std::find(names.begin(), names.end(), name); - return (it != names.end()) ? std::distance(names.begin(), it) : -1; - }; - posRegionIndex = findRegionIndex("refP"); - negRegionIndex = findRegionIndex("refN"); - fullRegionIndex = findRegionIndex("refFull"); - midRegionIndex = findRegionIndex("refMid"); - } - } - - static constexpr std::string_view FillTimeName[] = {"before/", "after/"}; - - enum QAFillTime { - kBefore, - kAfter - }; - - int getMagneticField(uint64_t timestamp) - { - static o2::parameters::GRPMagField* grpo = nullptr; - if (grpo == nullptr) { - // grpo = ccdb->getForTimeStamp("GLO/GRP/GRP", timestamp); - grpo = ccdb->getForTimeStamp("GLO/Config/GRPMagField", timestamp); - if (grpo == nullptr) { - LOGF(fatal, "GRP object not found for timestamp %llu", timestamp); - return 0; - } - LOGF(info, "Retrieved GRP for timestamp %llu with magnetic field of %d kG", timestamp, grpo->getNominalL3Field()); - } - return grpo->getNominalL3Field(); - } - - void loadCorrections(aod::BCsWithTimestamps::iterator const& bc) - { - uint64_t timestamp = bc.timestamp(); - if (!cfgRunByRun && cfg.correctionsLoaded) - return; - if (!cfgAcceptance.value.empty()) { - std::string runstr = (cfgRunByRun) ? "RunByRun/" : ""; - cfg.mAcceptance = ccdb->getForTimeStamp(cfgAcceptance.value + runstr, timestamp); - } - if (!cfgEfficiency.value.empty()) { - cfg.mEfficiency = ccdb->getForTimeStamp(cfgEfficiency, timestamp); - if (cfg.mEfficiency == nullptr) { - LOGF(fatal, "Could not load efficiency histogram from %s", cfgEfficiency.value.c_str()); - } - LOGF(info, "Loaded efficiency histogram from %s (%p)", cfgEfficiency.value.c_str(), (void*)cfg.mEfficiency); - } - cfg.correctionsLoaded = true; - } - - template - double getAcceptance(TTrack track, const double& vtxz) - { - double wacc = 1; - if (cfg.mAcceptance) - wacc = cfg.mAcceptance->getNUA(track.phi(), track.eta(), vtxz); - return wacc; - } - - template - double getEfficiency(TTrack track) - { - double eff = 1.; - if (cfg.mEfficiency) - eff = cfg.mEfficiency->GetBinContent(cfg.mEfficiency->FindBin(track.pt())); - if (eff == 0) - return -1.; - else - return 1. / eff; - } - - template - bool eventSelected(TCollision collision, const int& multTrk, const float& centrality, const int& run) - { - if (cfgTVXinTRD) { - if (collision.alias_bit(kTVXinTRD)) { - // TRD triggered - // "CMTVX-B-NOPF-TRD,minbias_TVX" - return 0; - } - if (cfgFillQA) { - registry.fill(HIST("eventQA/eventSel"), kTVXTRD); - } - if (cfgRunByRun && run != -1) - th1sList[run][hEventSel]->Fill(kTVXTRD); - } - if (cfgNoSameBunchPileupCut) { - if (!collision.selection_bit(o2::aod::evsel::kNoSameBunchPileup)) { - // rejects collisions which are associated with the same "found-by-T0" bunch crossing - // https://indico.cern.ch/event/1396220/#1-event-selection-with-its-rof - return 0; - } - if (cfgFillQA) { - registry.fill(HIST("eventQA/eventSel"), kNoSamebunchPU); - } - if (cfgRunByRun && run != -1) - th1sList[run][hEventSel]->Fill(kNoSamebunchPU); - } - if (cfgIsGoodZvtxFT0vsPV) { - if (!collision.selection_bit(o2::aod::evsel::kIsGoodZvtxFT0vsPV)) { - // removes collisions with large differences between z of PV by tracks and z of PV from FT0 A-C time difference - // use this cut at low multiplicities with caution - return 0; - } - if (cfgFillQA) { - registry.fill(HIST("eventQA/eventSel"), kZVtxFT0PV); - } - if (cfgRunByRun && run != -1) - th1sList[run][hEventSel]->Fill(kZVtxFT0PV); - } - if (cfgNoCollInTimeRangeStandard) { - if (!collision.selection_bit(o2::aod::evsel::kNoCollInTimeRangeStandard)) { - // Rejection of the collisions which have other events nearby - return 0; - } - if (cfgFillQA) { - registry.fill(HIST("eventQA/eventSel"), kNoCollTRStd); - } - if (cfgRunByRun && run != -1) - th1sList[run][hEventSel]->Fill(kNoCollTRStd); - } - - if (cfgIsVertexITSTPC) { - if (!collision.selection_bit(o2::aod::evsel::kIsVertexITSTPC)) { - // selects collisions with at least one ITS-TPC track, and thus rejects vertices built from ITS-only tracks - return 0; - } - if (cfgFillQA) { - registry.fill(HIST("eventQA/eventSel"), kVtxITSTPC); - } - if (cfgRunByRun && run != -1) - th1sList[run][hEventSel]->Fill(kVtxITSTPC); - } - - if (cfgIsGoodITSLayersAll) { - if (!collision.selection_bit(o2::aod::evsel::kIsGoodITSLayersAll)) { - return 0; - } - if (cfgFillQA) { - registry.fill(HIST("eventQA/eventSel"), kGoodITSLayers); - } - if (cfgRunByRun && run != -1) - th1sList[run][hEventSel]->Fill(kGoodITSLayers); - } - - float vtxz = -999; - if (collision.numContrib() > 1) { - vtxz = collision.posZ(); - float zRes = std::sqrt(collision.covZZ()); - float minZRes = 0.25; - int minNContrib = 20; - if (zRes > minZRes && collision.numContrib() < minNContrib) - vtxz = -999; - } - auto multNTracksPV = collision.multNTracksPV(); - - if (vtxz > o2::analysis::gfwflowflucpp::vtxZup || vtxz < o2::analysis::gfwflowflucpp::vtxZlow) - return 0; - - if (cfgMultCut && cfgUseAdditionalEventCut) { - if (multNTracksPV < fMultPVCutLow->Eval(centrality)) - return 0; - if (multNTracksPV > fMultPVCutHigh->Eval(centrality)) - return 0; - if (multTrk < fMultCutLow->Eval(centrality)) - return 0; - if (multTrk > fMultCutHigh->Eval(centrality)) - return 0; - if (multTrk > fMultPVGlobalCutHigh->Eval(collision.multNTracksPV())) - return 0; - - if (!(cfgGlobalAsideCorrCuts.cfgMultGlobalASideCorrCutFunction->empty()) && static_cast(collision.multFV0A()) < fMultGlobalV0ACutLow->Eval(multTrk)) - return 0; - if (!(cfgGlobalAsideCorrCuts.cfgMultGlobalASideCorrCutFunction->empty()) && static_cast(collision.multFV0A()) > fMultGlobalV0ACutHigh->Eval(multTrk)) - return 0; - if (!(cfgGlobalAsideCorrCuts.cfgMultGlobalASideCorrCutFunction->empty()) && static_cast(collision.multFT0A()) < fMultGlobalT0ACutLow->Eval(multTrk)) - return 0; - if (!(cfgGlobalAsideCorrCuts.cfgMultGlobalASideCorrCutFunction->empty()) && static_cast(collision.multFT0A()) > fMultGlobalT0ACutHigh->Eval(multTrk)) - return 0; - if (cfgRunByRun && run != -1) - th1sList[run][hEventSel]->Fill(kMultCuts); - } - if (cfgFillQA) { - registry.fill(HIST("eventQA/eventSel"), kMultCuts); - } - return 1; - } - - template - bool trackSelected(TTrack track) - { - if (cfgDCAxyNSigma && (std::fabs(track.dcaXY()) > fPtDepDCAxy->Eval(track.pt()))) - return false; - return ((track.tpcNClsCrossedRows() >= cfgNTPCXrows) && (track.tpcNClsFound() >= cfgNTPCCls) && (track.itsNCls() >= cfgMinNITSCls)); - } - - enum DataType { - kReco, - kGen - }; - - template - void fillWeights(const TTrack track, const double vtxz, const int& run) - { - if (cfgRunByRun) - th3sList[run][hNUAref]->Fill(track.phi(), track.eta(), vtxz); - else - registry.fill(HIST("phi_eta_vtxz_ref"), track.phi(), track.eta(), vtxz); - return; - } - - void createRunByRunHistograms(const int& run) - { - AxisSpec phiAxis = {o2::analysis::gfwflowflucpp::phibins, o2::analysis::gfwflowflucpp::philow, o2::analysis::gfwflowflucpp::phiup, "#phi"}; - AxisSpec etaAxis = {o2::analysis::gfwflowflucpp::etabins, -cfgEta, cfgEta, "#eta"}; - AxisSpec vtxAxis = {o2::analysis::gfwflowflucpp::vtxZbins, -cfgVtxZ, cfgVtxZ, "Vtx_{z} (cm)"}; - AxisSpec nchAxis = {o2::analysis::gfwflowflucpp::nchbins, o2::analysis::gfwflowflucpp::nchlow, o2::analysis::gfwflowflucpp::nchup, "N_{ch}"}; - AxisSpec centAxis = {o2::analysis::gfwflowflucpp::centbinning, "Centrality (%)"}; - std::vector> histos(kCount_TH1Names); - histos[hPhi] = registry.add(Form("%d/phi", run), "", {HistType::kTH1D, {phiAxis}}); - histos[hEta] = registry.add(Form("%d/eta", run), "", {HistType::kTH1D, {etaAxis}}); - histos[hVtxZ] = registry.add(Form("%d/vtxz", run), "", {HistType::kTH1D, {vtxAxis}}); - histos[hMult] = registry.add(Form("%d/mult", run), "", {HistType::kTH1D, {nchAxis}}); - histos[hCent] = registry.add(Form("%d/cent", run), "", {HistType::kTH1D, {centAxis}}); - if (cfgFillFlowRunByRun) { - std::vector> profiles(kCount_TProfileNames); - profiles[pfCorr22] = registry.add(Form("%d/corr22", run), "", {HistType::kTProfile, {(cfgUseNch) ? nchAxis : centAxis}}); - tpfsList.insert(std::make_pair(run, profiles)); - } - histos[hEventSel] = registry.add(Form("%d/eventSel", run), "Number of Events;; Counts", {HistType::kTH1D, {{11, 0.5, 11.5}}}); - histos[hEventSel]->GetXaxis()->SetBinLabel(kFilteredEvent, "Filtered event"); - histos[hEventSel]->GetXaxis()->SetBinLabel(kSel8, "sel8"); - histos[hEventSel]->GetXaxis()->SetBinLabel(kOccupancy, "occupancy"); - histos[hEventSel]->GetXaxis()->SetBinLabel(kTVXTRD, "kTVXinTRD"); - histos[hEventSel]->GetXaxis()->SetBinLabel(kNoSamebunchPU, "kNoSameBunchPileup"); - histos[hEventSel]->GetXaxis()->SetBinLabel(kZVtxFT0PV, "kIsGoodZvtxFT0vsPV"); - histos[hEventSel]->GetXaxis()->SetBinLabel(kNoCollTRStd, "kNoCollInTimeRangeStandard"); - histos[hEventSel]->GetXaxis()->SetBinLabel(kVtxITSTPC, "kIsVertexITSTPC"); - histos[hEventSel]->GetXaxis()->SetBinLabel(kGoodITSLayers, "kIsGoodITSLayersAll"); - histos[hEventSel]->GetXaxis()->SetBinLabel(kMultCuts, "after Mult cuts"); - histos[hEventSel]->GetXaxis()->SetBinLabel(kTrackCent, "has track + within cent"); - th1sList.insert(std::make_pair(run, histos)); - std::vector> histos3d(kCount_TH3Names); - histos3d[hNUAref] = registry.add(Form("%d/phi_eta_vtxz_ref", run), "", {HistType::kTH3D, {phiAxis, etaAxis, vtxAxis}}); - th3sList.insert(std::make_pair(run, histos3d)); - return; - } - - void addConfigObjectsToObjArray(TObjArray* oba, const std::vector& configs) - { - const auto shapeSel = getShapeSel(); - const int nQnOutputBins = getNQnOutputBins(); - for (auto it = configs.begin(); it != configs.end(); ++it) { - for (int jese = 0; jese < nQnOutputBins; ++jese) { - std::string name = Form("%s_%d_%s", shapeSel.c_str(), jese, it->Head.c_str()); - std::string title = it->Head + std::string("_ese"); - oba->Add(new TNamed(name.c_str(), title.c_str())); - } - } - } - - template - void fillOutputContainers(const float& centmult, const double& rndm, const int& qPtmp, const int& run = 0) - { - const auto shapeSel = getShapeSel(); - - for (uint l_ind = 0; l_ind < corrconfigs.size(); ++l_ind) { - auto dnx = fGFW->Calculate(corrconfigs.at(l_ind), 0, kTRUE).real(); - if (dnx == 0) - continue; - - auto val = fGFW->Calculate(corrconfigs.at(l_ind), 0, kFALSE).real() / dnx; - if (std::abs(val) >= 1.) - continue; - - std::string profName = Form("%s_%d_%s", shapeSel.c_str(), qPtmp, corrconfigs.at(l_ind).Head.c_str()); - - if constexpr (dt == kGen) { - fFCgen->FillProfile(profName.c_str(), - centmult, - val, - cfgUseMultiplicityFlowWeights ? dnx : 1.0, - rndm); - } else { - fFC->FillProfile(profName.c_str(), - centmult, - val, - cfgUseMultiplicityFlowWeights ? dnx : 1.0, - rndm); - - if (cfgRunByRun && cfgFillFlowRunByRun && l_ind == 0) { - tpfsList[run][pfCorr22]->Fill(centmult, - val, - cfgUseMultiplicityFlowWeights ? dnx : 1.0); - } - } - } - } - - struct XAxis { - float centrality; - int64_t multiplicity; - double time; - }; - - struct AcceptedTracks { - int nPos; - int nNeg; - int nFull; - int nMid; - }; - - 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 (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 && fabs(trk.eta())< cfgQnTrkAbsEtaMax) { - // 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 && fabs(trk.eta())< cfgQnTrkAbsEtaMax) { - // 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 (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]); - } - - return qvec; - } - - float computeqnVec(InTaskTPCQVector const& qvec, bool useNegativeEtaHalf) - { - if (qvec.nTrkTPCPos <= 0 || qvec.nTrkTPCNeg <= 0) { - return -1.f; - } - - 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 - std::vector> getQnBinSeparator2D(std::vector flat, const int numQnBins = 10) - { - size_t nBins = numQnBins + 1; - - if (flat.empty() || flat.size() % nBins != 0) { - LOGP(error, "ConfQnBinSeparator size = {} is not divisible by {}", - flat.size(), nBins); - return {{kInvalidQnSeparator, kInvalidQnSeparator}}; - } - - size_t nCent = flat.size() / nBins; - std::vector> res(nCent, std::vector(nBins)); - - for (size_t i = 0; i < nCent; ++i) { - for (size_t j = 0; j < nBins; ++j) { - res[i][j] = flat[i * nBins + j]; - } - } - return res; - } - - /// Get the bin number of qn-vector(FT0C) of an event - /// \param centBinWidth centrality bin width, example: per 1%, per 10% ... - /// \return bin number of qn-vector of the event - // add a param : bool doFillHisto ? - int myqnBin(float centrality, float centMax, float qn, std::vector qnBinSprt, const int numQnBins, float centBinWidth = 1.f) - { - auto twoDSeparator = getQnBinSeparator2D(qnBinSprt, numQnBins); - if (twoDSeparator.empty() || twoDSeparator[0][0] == kInvalidQnSeparator) { - LOGP(warning, "ConfQnBinSeparator not set, using default fallback!"); - return kInvalidQnBin; // safe fallback - } - - int qnBin = kInvalidQnBin; - int mycentBin = static_cast(centrality / centBinWidth); - if (mycentBin >= static_cast(centMax / centBinWidth)) - return qnBin; - - if (mycentBin > static_cast(twoDSeparator.size()) - 1) - return qnBin; - - for (int iqn(0); iqn < static_cast(twoDSeparator[mycentBin].size()) - 1; ++iqn) { - if (qn > twoDSeparator[mycentBin][iqn] && qn <= twoDSeparator[mycentBin][iqn + 1]) { - qnBin = iqn; - break; - } else { - continue; - } - } - - return qnBin; - } - - template - void processCollision(TCollision collision, TTracks tracks, const XAxis& xaxis, const int& run, const int& qPtmp) - { - if (tracks.size() < 1) - return; - if (dt != kGen && xaxis.centrality >= 0 && - (xaxis.centrality < o2::analysis::gfwflowflucpp::centbinning.front() || - xaxis.centrality > o2::analysis::gfwflowflucpp::centbinning.back())) - return; - if (xaxis.multiplicity < cfgFixedMultMin || xaxis.multiplicity > cfgFixedMultMax) - return; - - if (dt != kGen) { - if (cfgFillQA) { - registry.fill(HIST("eventQA/eventSel"), kTrackCent); - } - if (cfgRunByRun) - th1sList[run][hEventSel]->Fill(kTrackCent); - } - - if (cfgFillQA && xaxis.centrality >= 0) - registry.fill(HIST("eventQA/after/centrality"), xaxis.centrality); - if (cfgFillQA) - registry.fill(HIST("eventQA/after/multiplicity"), xaxis.multiplicity); - - float vtxz = collision.posZ(); - if (dt != kGen && cfgRunByRun) { - th1sList[run][hVtxZ]->Fill(vtxz); - th1sList[run][hMult]->Fill(xaxis.multiplicity); - th1sList[run][hCent]->Fill(xaxis.centrality); - } - - fGFW->Clear(); - float lRandom = fRndm->Rndm(); - - AcceptedTracks acceptedTracks{0, 0, 0, 0}; - - for (const auto& track : tracks) { - processTrack(track, vtxz, xaxis.multiplicity, run, acceptedTracks); - } - - if (cfgConsistentEventFlag & kRequireBothEtaSides) - if (!acceptedTracks.nPos || !acceptedTracks.nNeg) - return; - if (cfgConsistentEventFlag & kRequireFullFourParticleTracks) - if (acceptedTracks.nFull < kMinTracksForFourParticleCorrelation) - return; - if (cfgConsistentEventFlag & kRequireTwoTracksInBothEtaSides) - if (acceptedTracks.nPos < kMinTracksPerEtaSideForGapCorrelation || - acceptedTracks.nNeg < kMinTracksPerEtaSideForGapCorrelation) - return; - if (cfgConsistentEventFlag & kRequireTwoTracksInThreeEtaRegions) - if (acceptedTracks.nPos < kMinTracksPerEtaRegionForThreeSubevents || - acceptedTracks.nMid < kMinTracksPerEtaRegionForThreeSubevents || - acceptedTracks.nNeg < kMinTracksPerEtaRegionForThreeSubevents) - return; - - fillOutputContainers
(cfgUseNch ? static_cast(xaxis.multiplicity) : xaxis.centrality, - lRandom, qPtmp, run); - } - - template - void processGenCollision(TCollision collision, TParticles particles, const int& mcCollisionId, const XAxis& xaxis, const int& run, const int& qPtmp) - { - if (xaxis.multiplicity < cfgFixedMultMin || xaxis.multiplicity > cfgFixedMultMax) - return; - - if (cfgFillQA && xaxis.centrality >= 0) - registry.fill(HIST("eventQA/after/centrality"), xaxis.centrality); - if (cfgFillQA) - registry.fill(HIST("eventQA/after/multiplicity"), xaxis.multiplicity); - - fGFW->Clear(); - float lRandom = fRndm->Rndm(); - float vtxz = collision.posZ(); - - AcceptedTracks acceptedTracks{0, 0, 0, 0}; - for (const auto& particle : particles) { - if (particle.mcCollisionId() != mcCollisionId) - continue; - processTrack(particle, vtxz, xaxis.multiplicity, run, acceptedTracks); - } - - if (cfgConsistentEventFlag & kRequireBothEtaSides) - if (!acceptedTracks.nPos || !acceptedTracks.nNeg) - return; - if (cfgConsistentEventFlag & kRequireFullFourParticleTracks) - if (acceptedTracks.nFull < kMinTracksForFourParticleCorrelation) - return; - if (cfgConsistentEventFlag & kRequireTwoTracksInBothEtaSides) - if (acceptedTracks.nPos < kMinTracksPerEtaSideForGapCorrelation || - acceptedTracks.nNeg < kMinTracksPerEtaSideForGapCorrelation) - return; - if (cfgConsistentEventFlag & kRequireTwoTracksInThreeEtaRegions) - if (acceptedTracks.nPos < kMinTracksPerEtaRegionForThreeSubevents || - acceptedTracks.nMid < kMinTracksPerEtaRegionForThreeSubevents || - acceptedTracks.nNeg < kMinTracksPerEtaRegionForThreeSubevents) - return; - - fillOutputContainers(cfgUseNch ? static_cast(xaxis.multiplicity) : xaxis.centrality, - lRandom, qPtmp, run); - } - - bool isStable(int pdg) - { - if (std::abs(pdg) == PDG_t::kPiPlus) - return true; - if (std::abs(pdg) == PDG_t::kKPlus) - return true; - if (std::abs(pdg) == PDG_t::kProton) - return true; - if (std::abs(pdg) == PDG_t::kElectron) - return true; - if (std::abs(pdg) == PDG_t::kMuonMinus) - return true; - return false; - } - - template - void fillAcceptedTracks(TTrack track, AcceptedTracks& acceptedTracks) - { - if (posRegionIndex >= 0 && track.eta() > o2::analysis::gfwflowflucpp::regions.GetEtaMin()[posRegionIndex] && track.eta() < o2::analysis::gfwflowflucpp::regions.GetEtaMax()[posRegionIndex]) - ++acceptedTracks.nPos; - if (negRegionIndex >= 0 && track.eta() > o2::analysis::gfwflowflucpp::regions.GetEtaMin()[negRegionIndex] && track.eta() < o2::analysis::gfwflowflucpp::regions.GetEtaMax()[negRegionIndex]) - ++acceptedTracks.nNeg; - if (fullRegionIndex >= 0 && track.eta() > o2::analysis::gfwflowflucpp::regions.GetEtaMin()[fullRegionIndex] && track.eta() < o2::analysis::gfwflowflucpp::regions.GetEtaMax()[fullRegionIndex]) - ++acceptedTracks.nFull; - if (midRegionIndex >= 0 && track.eta() > o2::analysis::gfwflowflucpp::regions.GetEtaMin()[midRegionIndex] && track.eta() < o2::analysis::gfwflowflucpp::regions.GetEtaMax()[midRegionIndex]) - ++acceptedTracks.nMid; - } - - template - inline void processTrack(TTrack const& track, const float& vtxz, const int& multiplicity, const int& run, AcceptedTracks& acceptedTracks) - { - if constexpr (framework::has_type_v) { - if (track.mcParticleId() < 0 || !(track.has_mcParticle())) - return; - - auto mcParticle = track.mcParticle(); - if (!mcParticle.isPhysicalPrimary()) - return; - if (!isStable(mcParticle.pdgCode())) - return; - if (cfgFillQA && cfgIsMC) { - fillTrackQA(track, vtxz); - registry.fill(HIST("trackQA/before/nch_pt"), multiplicity, track.pt()); - } - if (!trackSelected(track)) - return; - - if (cfgFillWeights) { - fillWeights(track, vtxz, run); - } else { - fillGFW(track, vtxz); - fillAcceptedTracks(track, acceptedTracks); - } - - if (cfgFillQA) { - fillTrackQA(track, vtxz); - registry.fill(HIST("trackQA/after/nch_pt"), multiplicity, track.pt()); - if (cfgRunByRun) { - th1sList[run][hPhi]->Fill(track.phi()); - th1sList[run][hEta]->Fill(track.eta()); - } - } - - } else if constexpr (framework::has_type_v) { - if (!track.isPhysicalPrimary() || !isStable(track.pdgCode())) - return; - - fillGFW(track, vtxz); - fillAcceptedTracks(track, acceptedTracks); - if (cfgFillQA && cfgIsMC) { - fillTrackQA(track, vtxz); - registry.fill(HIST("MCGen/trackQA/nch_pt"), multiplicity, track.pt()); - } - } else { - if (cfgFillQA) { - fillTrackQA(track, vtxz); - registry.fill(HIST("trackQA/before/nch_pt"), multiplicity, track.pt()); - } - if (!trackSelected(track)) - return; - - if (cfgFillWeights) { - fillWeights(track, vtxz, run); - } else { - fillGFW(track, vtxz); - fillAcceptedTracks(track, acceptedTracks); - } - if (cfgFillQA) { - fillTrackQA(track, vtxz); - registry.fill(HIST("trackQA/after/nch_pt"), multiplicity, track.pt()); - if (cfgRunByRun) { - th1sList[run][hPhi]->Fill(track.phi()); - th1sList[run][hEta]->Fill(track.eta()); - th3sList[run][hNUAref]->Fill(track.phi(), track.eta(), vtxz, getAcceptance(track, vtxz)); - } - } - } - return; - } - - template - inline void fillGFW(TTrack track, const double& vtxz) - { - bool withinPtRef = (track.pt() > o2::analysis::gfwflowflucpp::ptreflow && track.pt() < o2::analysis::gfwflowflucpp::ptrefup); - bool withinPtPOI = (track.pt() > o2::analysis::gfwflowflucpp::ptpoilow && track.pt() < o2::analysis::gfwflowflucpp::ptpoiup); - if (!withinPtPOI && !withinPtRef) - return; - double weff = (dt == kGen) ? 1. : getEfficiency(track); - if (weff < 0) - return; - - double wacc = (dt == kGen) ? 1. : getAcceptance(track, vtxz); - if (withinPtRef) - fGFW->Fill(track.eta(), fSecondAxis->FindBin(track.pt()) - 1, track.phi(), weff * wacc, 1); - if (withinPtPOI) - fGFW->Fill(track.eta(), fSecondAxis->FindBin(track.pt()) - 1, track.phi(), weff * wacc, 2); - if (withinPtRef && withinPtPOI) - fGFW->Fill(track.eta(), fSecondAxis->FindBin(track.pt()) - 1, track.phi(), weff * wacc, 4); - return; - } - - template - inline void fillTrackQA(TTrack track, const float vtxz) - { - if constexpr (dt == kGen) { - registry.fill(HIST("MCGen/trackQA/phi_eta_vtxZ"), track.phi(), track.eta(), vtxz); - registry.fill(HIST("MCGen/trackQA/pt_ref"), track.pt()); - registry.fill(HIST("MCGen/trackQA/pt_poi"), track.pt()); - } else { - double wacc = getAcceptance(track, vtxz); - registry.fill(HIST("trackQA/") + HIST(FillTimeName[ft]) + HIST("phi_eta_vtxZ"), track.phi(), track.eta(), vtxz, (ft == kAfter) ? wacc : 1.0); - registry.fill(HIST("trackQA/") + HIST(FillTimeName[ft]) + HIST("pt_dcaXY_dcaZ"), track.pt(), track.dcaXY(), track.dcaZ()); - - registry.fill(HIST("trackQA/") + HIST(FillTimeName[ft]) + HIST("chi2prTPCcls"), track.tpcChi2NCl()); - registry.fill(HIST("trackQA/") + HIST(FillTimeName[ft]) + HIST("chi2prITScls"), track.itsChi2NCl()); - registry.fill(HIST("trackQA/") + HIST(FillTimeName[ft]) + HIST("nTPCClusters"), track.tpcNClsFound()); - registry.fill(HIST("trackQA/") + HIST(FillTimeName[ft]) + HIST("nITSClusters"), track.itsNCls()); - registry.fill(HIST("trackQA/") + HIST(FillTimeName[ft]) + HIST("nTPCCrossedRows"), track.tpcNClsCrossedRows()); - - if (ft == kAfter) { - registry.fill(HIST("trackQA/") + HIST(FillTimeName[ft]) + HIST("pt_ref"), track.pt()); - registry.fill(HIST("trackQA/") + HIST(FillTimeName[ft]) + HIST("pt_poi"), track.pt()); - } - } - } - - template - float getCentrality(TCollision collision) - { - switch (cfgCentEstimator) { - case kCentFT0C: - return collision.centFT0C(); - case kCentFT0CVariant1: - return collision.centFT0CVariant1(); - case kCentFT0M: - return collision.centFT0M(); - case kCentFV0A: - return collision.centFV0A(); - case kCentNTPV: - return collision.centNTPV(); - case kCentNGlobal: - return collision.centNGlobal(); - case kCentMFT: - return collision.centMFT(); - default: - return collision.centFT0C(); - } - } - - template - inline void fillEventQA(TCollision collision, XAxis xaxis) - { - if constexpr (framework::has_type_v) { - registry.fill(HIST("eventQA/") + HIST(FillTimeName[ft]) + HIST("globalTracks_centT0C"), collision.centFT0C(), xaxis.multiplicity); - registry.fill(HIST("eventQA/") + HIST(FillTimeName[ft]) + HIST("PVTracks_centT0C"), collision.centFT0C(), collision.multNTracksPV()); - registry.fill(HIST("eventQA/") + HIST(FillTimeName[ft]) + HIST("multT0C_centT0C"), collision.centFT0C(), collision.multFT0C()); - registry.fill(HIST("eventQA/") + HIST(FillTimeName[ft]) + HIST("centT0M_centT0C"), collision.centFT0C(), collision.centFT0M()); - registry.fill(HIST("eventQA/") + HIST(FillTimeName[ft]) + HIST("centV0A_centT0C"), collision.centFT0C(), collision.centFV0A()); - registry.fill(HIST("eventQA/") + HIST(FillTimeName[ft]) + HIST("centGlobal_centT0C"), collision.centFT0C(), collision.centNGlobal()); - registry.fill(HIST("eventQA/") + HIST(FillTimeName[ft]) + HIST("centNTPV_centT0C"), collision.centFT0C(), collision.centNTPV()); - registry.fill(HIST("eventQA/") + HIST(FillTimeName[ft]) + HIST("centMFT_centT0C"), collision.centFT0C(), collision.centMFT()); - } - registry.fill(HIST("eventQA/") + HIST(FillTimeName[ft]) + HIST("globalTracks_PVTracks"), collision.multNTracksPV(), xaxis.multiplicity); - registry.fill(HIST("eventQA/") + HIST(FillTimeName[ft]) + HIST("globalTracks_multT0A"), collision.multFT0A(), xaxis.multiplicity); - registry.fill(HIST("eventQA/") + HIST(FillTimeName[ft]) + HIST("globalTracks_multV0A"), collision.multFV0A(), xaxis.multiplicity); - registry.fill(HIST("eventQA/") + HIST(FillTimeName[ft]) + HIST("multV0A_multT0A"), collision.multFT0A(), collision.multFV0A()); - if (cfgTimeDependent) { - registry.fill(HIST("eventQA/") + HIST(FillTimeName[ft]) + HIST("multiplicity_time"), xaxis.time, xaxis.multiplicity); - registry.fill(HIST("eventQA/") + HIST(FillTimeName[ft]) + HIST("multT0C_time"), xaxis.time, collision.multFT0C()); - registry.fill(HIST("eventQA/") + HIST(FillTimeName[ft]) + HIST("multT0A_time"), xaxis.time, collision.multFT0A()); - registry.fill(HIST("eventQA/") + HIST(FillTimeName[ft]) + HIST("multV0A_time"), xaxis.time, collision.multFV0A()); - registry.fill(HIST("eventQA/") + HIST(FillTimeName[ft]) + HIST("multPV_time"), xaxis.time, collision.multNTracksPV()); - } - return; - } - - double getTimeSinceStartOfFill(uint64_t timestamp, int firstRun) - { - auto runDuration = ccdb->getRunDuration(firstRun); - uint64_t tsSOF = runDuration.first; - uint64_t diff = timestamp - tsSOF; - return static_cast(diff) / 3600000.0; - } - - void fillQnEventCounter(float count) - { - const int qnHarmonic = static_cast(cfgQnSelectionHarmonic); - - if (qnHarmonic == EllipticQVectorHarmonic) { - registry.fill(HIST("mq2/eventcounter"), count); - } else if (qnHarmonic == TriangularQVectorHarmonic) { - registry.fill(HIST("mq3/eventcounter"), count); - } - } - - void fillQnCalibrationHistograms(float centrality, float multiplicity, float qnPos, float qnNeg) - { - const int qnHarmonic = static_cast(cfgQnSelectionHarmonic); - - if (qnHarmonic == EllipticQVectorHarmonic) { - registry.fill(HIST("mq2/h2_cent_q2_etapos"), centrality, qnPos); - registry.fill(HIST("mq2/h2_cent_q2_etaneg"), centrality, qnNeg); - registry.fill(HIST("mq2/h2_mult_q2_etapos"), multiplicity, qnPos); - registry.fill(HIST("mq2/h2_mult_q2_etaneg"), multiplicity, qnNeg); - } else if (qnHarmonic == TriangularQVectorHarmonic) { - registry.fill(HIST("mq3/h2_cent_q3_etapos"), centrality, qnPos); - registry.fill(HIST("mq3/h2_cent_q3_etaneg"), centrality, qnNeg); - registry.fill(HIST("mq3/h2_mult_q3_etapos"), multiplicity, qnPos); - registry.fill(HIST("mq3/h2_mult_q3_etaneg"), multiplicity, qnNeg); - } - } - - void processData(soa::Filtered>::iterator const& collision, - aod::BCsWithTimestamps const&, GFWTracks const& tracks) - { - auto bc = collision.bc_as(); - int run = bc.runNumber(); - if (run != lastRun) { - lastRun = run; - LOGF(info, "run = %d", run); - if (cfgRunByRun) { - if (std::find(runNumbers.begin(), runNumbers.end(), run) == runNumbers.end()) { - LOGF(info, "Creating histograms for run %d", run); - createRunByRunHistograms(run); - runNumbers.push_back(run); - } - if (!cfgFillWeights) - loadCorrections(bc); - } - } - if (!cfgFillWeights && !cfgRunByRun) - loadCorrections(bc); - - registry.fill(HIST("eventQA/eventSel"), kFilteredEvent); - if (cfgRunByRun) - th1sList[run][hEventSel]->Fill(kFilteredEvent); - - if (!collision.sel8()) - return; - - registry.fill(HIST("eventQA/eventSel"), kSel8); - if (cfgRunByRun) - th1sList[run][hEventSel]->Fill(kSel8); - - if (cfgDoOccupancySel) { - int occupancy = collision.trackOccupancyInTimeRange(); - if (occupancy < 0 || occupancy > cfgOccupancySelection) - return; - } - - registry.fill(HIST("eventQA/eventSel"), kOccupancy); - if (cfgRunByRun) - th1sList[run][hEventSel]->Fill(kOccupancy); - - const XAxis xaxis{ - getCentrality(collision), - collision.multNTracksPV(), - (cfgTimeDependent) ? getTimeSinceStartOfFill(bc.timestamp(), *firstRunOfCurrentFill) : -1.0}; - - if (cfgTimeDependent && run == *firstRunOfCurrentFill && - firstRunOfCurrentFill != o2::analysis::gfwflowflucpp::firstRunsOfFill.end() - 1) - ++firstRunOfCurrentFill; - - if (cfgFillQA) - fillEventQA(collision, xaxis); - - registry.fill(HIST("eventQA/before/centrality"), xaxis.centrality); - registry.fill(HIST("eventQA/before/multiplicity"), xaxis.multiplicity); - - if (!eventSelected(collision, xaxis.multiplicity, xaxis.centrality, run)) - return; - - if (cfgFillQA) - fillEventQA(collision, xaxis); - - int qPtmp = 0; - if (!static_cast(cfgBypassQnSelection)) { - const auto qvecTPC = calcQVec(static_cast(cfgQnSelectionHarmonic), tracks, collision); - float qn = computeqnVec(qvecTPC, cfgUseNegativeEtaHalfForq2); - if (qn < 0) - return; - - qPtmp = myqnBin(cfgEvtSelCent ? xaxis.centrality : xaxis.multiplicity, - cfgCentMax, qn, qnBinSeparator, cfgNumQnBins); - if (qPtmp < 0) - return; - } - - processCollision(collision, tracks, xaxis, run, qPtmp); - } - PROCESS_SWITCH(FlowFlucGfwPp, processData, "Process analysis for non-derived data", false); - - void processq2(soa::Filtered>::iterator const& collision, aod::BCsWithTimestamps const&, GFWTracks const& tracks) - { - float count{0.5}; - fillQnEventCounter(count++); - if (!collision.sel8()) { - return; - } - fillQnEventCounter(count++); - if (cfgDoOccupancySel) { - int occupancy = collision.trackOccupancyInTimeRange(); - if (occupancy < 0 || occupancy > cfgOccupancySelection) { - return; - } - } - fillQnEventCounter(count++); - - 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 qvecTPC = calcQVec(static_cast(cfgQnSelectionHarmonic), 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; - } - - fillQnEventCounter(count++); - fillQnCalibrationHistograms(centr, multi, qvecPos, qvecNeg); - } - PROCESS_SWITCH(FlowFlucGfwPp, processq2, "Process analysis for filling q_n-vector calibration histograms", true); - - void processMC(soa::Filtered>::iterator const& collision, - aod::BCsWithTimestamps const&, GFWTracksMC const& tracks, aod::McCollisions const&, aod::McParticles const& mcParticles) - { - auto bc = collision.bc_as(); - int run = bc.runNumber(); - if (run != lastRun) { - lastRun = run; - LOGF(info, "run = %d", run); - if (cfgRunByRun) { - if (std::find(runNumbers.begin(), runNumbers.end(), run) == runNumbers.end()) { - LOGF(info, "Creating histograms for run %d", run); - createRunByRunHistograms(run); - runNumbers.push_back(run); - } - if (!cfgFillWeights) - loadCorrections(bc); - } - } - if (!cfgFillWeights && !cfgRunByRun) - loadCorrections(bc); - - registry.fill(HIST("eventQA/eventSel"), kFilteredEvent); - if (cfgRunByRun) - th1sList[run][hEventSel]->Fill(kFilteredEvent); - - if (!collision.sel8()) - return; - - registry.fill(HIST("eventQA/eventSel"), kSel8); - if (cfgRunByRun) - th1sList[run][hEventSel]->Fill(kSel8); - - if (cfgDoOccupancySel) { - int occupancy = collision.trackOccupancyInTimeRange(); - if (occupancy < 0 || occupancy > cfgOccupancySelection) - return; - } - - registry.fill(HIST("eventQA/eventSel"), kOccupancy); - if (cfgRunByRun) - th1sList[run][hEventSel]->Fill(kOccupancy); - - const XAxis xaxis{ - getCentrality(collision), - collision.multNTracksPV(), - (cfgTimeDependent) ? getTimeSinceStartOfFill(bc.timestamp(), *firstRunOfCurrentFill) : -1.0}; - - if (cfgTimeDependent && run == *firstRunOfCurrentFill && - firstRunOfCurrentFill != o2::analysis::gfwflowflucpp::firstRunsOfFill.end() - 1) - ++firstRunOfCurrentFill; - - if (cfgFillQA) - fillEventQA(collision, xaxis); - - registry.fill(HIST("eventQA/before/centrality"), xaxis.centrality); - registry.fill(HIST("eventQA/before/multiplicity"), xaxis.multiplicity); - - if (!eventSelected(collision, xaxis.multiplicity, xaxis.centrality, run)) - return; - - if (cfgFillQA) - fillEventQA(collision, xaxis); - - processCollision(collision, tracks, xaxis, run, 0); - - if (!collision.has_mcCollision()) - return; - - const auto genCollision = collision.template mcCollision_as(); - processGenCollision(genCollision, mcParticles, collision.mcCollisionId(), xaxis, run, 0); - } - PROCESS_SWITCH(FlowFlucGfwPp, processMC, "Process analysis for Monte-Carlo data", false); - -}; - -WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) -{ - return WorkflowSpec{ - adaptAnalysisTask(cfgc), - }; -} From 56f4d8a115537d8bf189bd910cc8d7718ce1c67a Mon Sep 17 00:00:00 2001 From: ALICE Action Bot Date: Mon, 22 Jun 2026 15:03:14 +0000 Subject: [PATCH 3/4] Please consider the following formatting changes --- PWGCF/Flow/Tasks/flowFlucGfwPp.cxx | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/PWGCF/Flow/Tasks/flowFlucGfwPp.cxx b/PWGCF/Flow/Tasks/flowFlucGfwPp.cxx index b1f3d561cfc..2e015e6d526 100644 --- a/PWGCF/Flow/Tasks/flowFlucGfwPp.cxx +++ b/PWGCF/Flow/Tasks/flowFlucGfwPp.cxx @@ -893,7 +893,6 @@ struct FlowFlucGfwPp { std::string name = Form("%s_%d_%s", shapeSel.c_str(), jese, it->Head.c_str()); std::string title = it->Head + std::string("_ese"); oba->Add(new TNamed(name.c_str(), title.c_str())); - } } } @@ -1016,14 +1015,14 @@ struct FlowFlucGfwPp { registry.fill(HIST("qvecQA/ChTracks"), trk.pt(), trk.eta(), trk.phi()); } - if (trk.eta() > 0 && fabs(trk.eta())< cfgQnTrkAbsEtaMax) { + if (trk.eta() > 0 && fabs(trk.eta()) < cfgQnTrkAbsEtaMax) { // 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 && fabs(trk.eta())< cfgQnTrkAbsEtaMax) { + } else if (trk.eta() < 0 && fabs(trk.eta()) < cfgQnTrkAbsEtaMax) { // 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); @@ -1557,10 +1556,10 @@ struct FlowFlucGfwPp { PROCESS_SWITCH(FlowFlucGfwPp, processq2, "Process analysis for filling q_n-vector calibration histograms", true); void processMC(soa::Filtered>::iterator const& collision, - aod::BCsWithTimestamps const&, GFWTracksMC const& tracks, aod::McCollisions const&, aod::McParticles const& mcParticles) + aod::CentFT0Cs, aod::CentFT0CVariant1s, aod::CentFT0Ms, + aod::CentFV0As, aod::CentNTPVs, aod::CentNGlobals, + aod::McCollisionLabels>>::iterator const& collision, + aod::BCsWithTimestamps const&, GFWTracksMC const& tracks, aod::McCollisions const&, aod::McParticles const& mcParticles) { auto bc = collision.bc_as(); int run = bc.runNumber(); @@ -1631,7 +1630,6 @@ struct FlowFlucGfwPp { processGenCollision(genCollision, mcParticles, collision.mcCollisionId(), xaxis, run, 0); } PROCESS_SWITCH(FlowFlucGfwPp, processMC, "Process analysis for Monte-Carlo data", false); - }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) From 7a785ee664e4f85fc8b70b7939387b3761b401c1 Mon Sep 17 00:00:00 2001 From: wenyaCern <31894577+wenyaCern@users.noreply.github.com> Date: Tue, 23 Jun 2026 16:55:27 +0200 Subject: [PATCH 4/4] fix the o2-linter error Use std::fabs for absolute value of eta --- PWGCF/Flow/Tasks/flowFlucGfwPp.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/PWGCF/Flow/Tasks/flowFlucGfwPp.cxx b/PWGCF/Flow/Tasks/flowFlucGfwPp.cxx index 2e015e6d526..9a253ab83b8 100644 --- a/PWGCF/Flow/Tasks/flowFlucGfwPp.cxx +++ b/PWGCF/Flow/Tasks/flowFlucGfwPp.cxx @@ -1015,14 +1015,14 @@ struct FlowFlucGfwPp { registry.fill(HIST("qvecQA/ChTracks"), trk.pt(), trk.eta(), trk.phi()); } - if (trk.eta() > 0 && fabs(trk.eta()) < cfgQnTrkAbsEtaMax) { + if (trk.eta() > 0 && std::fabs(trk.eta()) < cfgQnTrkAbsEtaMax) { // 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 && fabs(trk.eta()) < cfgQnTrkAbsEtaMax) { + } else if (trk.eta() < 0 && std::fabs(trk.eta()) < cfgQnTrkAbsEtaMax) { // 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);