diff --git a/PWGCF/EbyEFluctuations/Tasks/radialFlowDecorr.cxx b/PWGCF/EbyEFluctuations/Tasks/radialFlowDecorr.cxx index cded5d14896..ba1aa5fa2f6 100644 --- a/PWGCF/EbyEFluctuations/Tasks/radialFlowDecorr.cxx +++ b/PWGCF/EbyEFluctuations/Tasks/radialFlowDecorr.cxx @@ -18,8 +18,6 @@ #include "Common/DataModel/EventSelection.h" #include "Common/DataModel/FT0Corrected.h" #include "Common/DataModel/Multiplicity.h" -#include "Common/DataModel/PIDResponseTOF.h" -#include "Common/DataModel/PIDResponseTPC.h" #include "Common/DataModel/TrackSelectionTables.h" #include @@ -38,8 +36,6 @@ #include #include -#include -#include #include #include #include @@ -47,6 +43,7 @@ #include #include #include +#include #include #include @@ -68,19 +65,17 @@ using namespace constants::math; struct RadialFlowDecorr { - static constexpr int KPidPionOne = 1; - static constexpr int KPidKaonTwo = 2; - static constexpr int KPidProtonThree = 3; - static constexpr int KConstTen = 10; - + // --- fixed constants --------------------------------------------------------- static constexpr int KnFt0cCell = 96; - static constexpr int KIntM = 3; - static constexpr int KIntK = 3; - static constexpr int KNEta = 9; + static constexpr int KIntM = 3; // pT-moment order used in the sums + static constexpr int KIntK = 3; // weight-power order used in the sums + + // KNEtaMax sizes every fixed-length eta array. The *active* count nEta is a + // runtime value: 17 for the 0.1-wide observable binning (16 bins + the index-0 + // full-range reference bin), 9 for the 0.2-wide binning (8 bins + reference). + static constexpr int KNEtaMax = 17; + static constexpr float KFloatEpsilon = 1e-6f; - static constexpr int KPiPlus = 211; - static constexpr int KKPlus = 321; - static constexpr int KProton = 2212; static constexpr float KBinOffset = 0.5f; static constexpr float KPhiMin = 0.f; static constexpr int KNbinsZvtx = 240; @@ -94,47 +89,18 @@ struct RadialFlowDecorr { static constexpr float KEtaMin = -1.2f; static constexpr float KEtaMax = 1.2f; static constexpr int KNbinsPhi = 64; - static constexpr float KEtaAxisMin = -0.8f; - static constexpr float KEtaAxisMax = 0.8f; static constexpr int KNbinsPtRes = 50; - static constexpr float KPtResMax = 1.f; static constexpr int KNbinsEtaRes = 100; - static constexpr float KEtaResMax = 0.5f; static constexpr int KNbinsVz = 80; static constexpr float KVzMin = -40.f; static constexpr float KVzMax = 40.f; - static constexpr float KVzResMax = 20.f; static constexpr int KNbinsEtaFine = 20; static constexpr float KEtaFineMax = 1.f; - static constexpr int KNbinsDca = 400; - static constexpr float KDcaMax = 0.2f; static constexpr float KCentMax = 90; - enum PIDIdx { - kInclusiveIdx = 0, - kPiMinusIdx, - kPiPlusIdx, - kPiAllIdx, - kKaMinusIdx, - kKaPlusIdx, - kKaAllIdx, - kAntiPrIdx, - kPrIdx, - kPrAllIdx, - KNsp - }; - - inline static const std::vector pidSuffix = {"", "_PiMinus", "_PiPlus", "_PiAll", "_KaMinus", "_KaPlus", "_KaAll", "_AntiPr", "_Pr", "_PrAll"}; - - struct PIDMeanSigmaMap { - static constexpr int MaxCentBins = 100; - double meanTOF[KNsp][MaxCentBins] = {{0.0}}; - double sigmaTOF[KNsp][MaxCentBins] = {{1.0}}; // Default sigma = 1 - double meanTPC[KNsp][MaxCentBins] = {{0.0}}; - double sigmaTPC[KNsp][MaxCentBins] = {{1.0}}; // Default sigma = 1 - }; - - PIDMeanSigmaMap* pidMeanSigmaMap = nullptr; + // Bootstrap: KMaxBoot is the compile-time storage cap; the number actually + // filled is nBoot = min(cfgNBootstrap, KMaxBoot). + static constexpr int KMaxBoot = 64; enum ECentralityEstimator { kCentFT0C = 1, @@ -148,21 +114,43 @@ struct RadialFlowDecorr { kOO = 3, kpp = 4 }; + + // Systematic-variation selector. Base is the main measurement; the others each + // read their own DataMean object for the correlation step and run at 0.2. + enum ESystType { + kSystBase = 0, + kSystDCA, + kSystEff, + kSystFlat, + kSystNEta, + kSystNITS, + kSystNTPC, + kSystPileup, + kSystVz, + kNSystType + }; + // Suffix appended to the DataMean CCDB path per systematic (Base -> no suffix). + inline static const std::vector systSuffix = { + "", "_systDCA", "_systEff", "_systFlat", "_systNEta", + "_systNITS", "_systNTPC", "_systPileup", "_systVz"}; + static constexpr float KinvalidCentrality = -1.0f; - inline static const std::vector etaLw = { - -0.8, - -0.8, -0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6}; - inline static const std::vector etaUp = { - 0.8, - -0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8}; - - Configurable cfgVtxZCut{"cfgVtxZCut", 10.f, "z-vertex range"}; - Configurable cfgPtMin{"cfgPtMin", 0.2f, "min pT"}; - Configurable cfgPtMax{"cfgPtMax", 5.0f, "max pT"}; - Configurable cfgEtaCut{"cfgEtaCut", 0.8f, "|η| cut"}; - Configurable cfgCutVertex{"cfgCutVertex", 10.0f, "Accepted z-vertex range"}; + + // Observable eta binning, filled at init(): index 0 = full-range reference bin. + std::vector etaLw; + std::vector etaUp; + int nEta = 9; // active eta-bin count (set in init()) + + int nBoot = 0; // active bootstrap samples (set in init()) + bool doBoot = false; // bootstrap active for this run (base data fluc only) + + // --- configurables ----------------------------------------------------------- + Configurable cfgVtxZCut{"cfgVtxZCut", 10.f, "|z_{vtx}| acceptance (cm): collision filter + explicit event/particle vertex checks"}; + Configurable cfgPtMin{"cfgPtMin", 0.2f, "min pT for observables"}; + Configurable cfgPtMax{"cfgPtMax", 5.0f, "max pT for observables"}; + Configurable cfgEtaCut{"cfgEtaCut", 0.8f, "|eta| cut"}; Configurable cfgCutTracKDcaMaxZ{"cfgCutTracKDcaMaxZ", 2.0f, "Maximum DcaZ"}; - Configurable cfgCutTracKDcaMaxXY{"cfgCutTracKDcaMaxXY", 0.2f, "Maximum DcaZ"}; + Configurable cfgCutTracKDcaMaxXY{"cfgCutTracKDcaMaxXY", 0.2f, "Maximum DcaXY"}; Configurable cfgPtDepDCAxy{"cfgPtDepDCAxy", false, "Use pt-dependent DCAxy cut"}; Configurable cfgDcaXyP0{"cfgDcaXyP0", 0.0026f, "p0 for DCAxy"}; @@ -177,18 +165,11 @@ struct RadialFlowDecorr { Configurable cfgITScluster{"cfgITScluster", 1, "Minimum Number of ITS cluster"}; Configurable cfgTPCcluster{"cfgTPCcluster", 80, "Minimum Number of TPC cluster"}; Configurable cfgTPCnCrossedRows{"cfgTPCnCrossedRows", 70, "Minimum Number of TPC crossed-rows"}; - Configurable cfgCutPtUpperTPC{"cfgCutPtUpperTPC", 0.6f, "Upper pT cut for PID using TPC only"}; - Configurable cfgnSigmaOtherParticles{"cfgnSigmaOtherParticles", 3.0f, "PID nSigma cut to remove other particles (default:3)"}; - Configurable cfgnSigmaCutTPC{"cfgnSigmaCutTPC", 2.0f, "PID nSigma cut for TPC"}; - Configurable cfgnSigmaCutTOF{"cfgnSigmaCutTOF", 2.0f, "PID nSigma cut for TOF"}; - Configurable cfgnSigmaCutCombTPCTOF{"cfgnSigmaCutCombTPCTOF", 2.0f, "PID nSigma combined cut for TPC and TOF"}; - - Configurable cfgCutPtLower{"cfgCutPtLower", 0.2f, "Lower pT cut"}; - Configurable cfgCutPtUpper{"cfgCutPtUpper", 10.0f, "Higher pT cut for inclusive hadron analysis"}; - Configurable cfgCutPtUpperPID{"cfgCutPtUpperPID", 6.0f, "Higher pT cut for identified particle analysis"}; + + Configurable cfgCutPtLower{"cfgCutPtLower", 0.2f, "Lower pT cut (track selection)"}; + Configurable cfgCutPtUpper{"cfgCutPtUpper", 10.0f, "Higher pT cut (track selection)"}; Configurable cfgCutEta{"cfgCutEta", 0.8f, "absolute Eta cut"}; - Configurable cfgMinTracksPerEtaBin{"cfgMinTracksPerEtaBin", 0, "Min weighted-track sum required in every narrow eta bin for inclusive species (0 = disabled)"}; - Configurable cfgNsubsample{"cfgNsubsample", 10, "Number of subsamples"}; + Configurable cfgMinTracksPerEtaBin{"cfgMinTracksPerEtaBin", 2, "Min weighted-track sum required in every narrow eta bin (systNEta; 0 = disabled)"}; Configurable cfgCentralityChoice{"cfgCentralityChoice", 1, "Which centrality estimator? 1-->FT0C, 2-->FT0M, 3-->FDDM, 4-->FV0A"}; Configurable cfgEvSelNoSameBunchPileup{"cfgEvSelNoSameBunchPileup", true, "Pileup removal"}; Configurable cfgUseGoodITSLayerAllCut{"cfgUseGoodITSLayerAllCut", true, "Remove time interval with dead ITS zone"}; @@ -205,7 +186,11 @@ struct RadialFlowDecorr { Configurable cfgNchPbMax{"cfgNchPbMax", 5000, "Max Nch range for PbPb collisions"}; Configurable cfgNchOMax{"cfgNchOMax", 800, "Max Nch range for OO collisions"}; - Configurable cfgSys{"cfgSys", 1, "Efficiency to be used for which system? 1-->PbPb, 2-->NeNe, 3-->OO, 4-->pp"}; + Configurable cfgSys{"cfgSys", 1, "Which collision system? 1-->PbPb, 2-->NeNe, 3-->OO, 4-->pp"}; + Configurable cfgSystType{"cfgSystType", 0, "Systematic variation: 0=Base,1=systDCA,2=systEff,3=systFlat,4=systNEta,5=systNITS,6=systNTPC,7=systPileup,8=systVz"}; + Configurable cfgNBootstrap{"cfgNBootstrap", 30, "Number of Poisson bootstrap samples (base data run only)"}; + Configurable cfgBootstrapSeed{"cfgBootstrapSeed", 0, "TRandom3 seed for bootstrap (0 = machine-random per job)"}; + Configurable cfgFlat{"cfgFlat", false, "Whether to use flattening weights"}; Configurable cfgEff{"cfgEff", false, "Whether to use Efficiency weights"}; Configurable cfgZDC{"cfgZDC", false, "Whether to use ZDC for pileup histograms"}; @@ -215,362 +200,184 @@ struct RadialFlowDecorr { ConfigurableAxis cfgAxisCent{"cfgAxisCent", {0.0, 1.0, 5.0, 10, 20, 40, 60, 80, 100}, "centrality axis (percentile)"}; - const AxisSpec centAxis{cfgAxisCent, "Centrality (%)"}; - const AxisSpec centAxis1Per{100, 0.0, 100.0, "Centrality (%)"}; + // --- axes -------------------------------------------------------------------- + AxisSpec centAxis{cfgAxisCent, "Centrality (%)"}; + AxisSpec centAxis1Per{100, 0.0, 100.0, "Centrality (%)"}; AxisSpec nChAxis{1, 0., 1., "Nch", "Nch"}; AxisSpec nChAxis2{1, 0., 1., "Nch", "Nch"}; - const AxisSpec vzAxis{5, -12.5, 12.5, "Vz"}; - const AxisSpec chgAxis{3, -1.5, 1.5}; - const AxisSpec pTAxis{{0.0, 0.2, 0.4, 0.6, 0.8, 1, 3, 5, 7, 10}, "pT Axis"}; - const AxisSpec etaAxis{{-0.8, -0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8}, "Eta"}; - const AxisSpec phiAxis{KNbinsPhi, KPhiMin, TwoPI, "#phi"}; - const AxisSpec etaBinAxis{KNEta + 1, -0.5, KNEta + 0.5, "#eta bin Number"}; - const AxisSpec spBinAxis{KNsp + 1, -KBinOffset, static_cast(KNsp) + KBinOffset, "species index Number"}; + AxisSpec vzAxis{5, -12.5, 12.5, "Vz"}; + AxisSpec chgAxis{3, -1.5, 1.5}; + AxisSpec pTAxis{{0.0, 0.2, 0.4, 0.6, 0.8, 1, 3, 5, 7, 10}, "pT Axis"}; + AxisSpec phiAxis{KNbinsPhi, KPhiMin, TwoPI, "#phi"}; - const AxisSpec gapAxis{{-1.5, -1.3, -1.1, -0.9, -0.7, -0.5, -0.3, -0.1, - 0.1, 0.3, 0.5, 0.7, 0.9, 1.1, 1.3, 1.5}, - "Gap"}; + // etaFlatAxis is the *fixed* granularity used only for the flattening map, so + // that a flattening map is reusable across the 0.1/0.2 observable binnings. + AxisSpec etaFlatAxis{{-0.8, -0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8}, "#eta"}; - const AxisSpec sumAxis{{-1.5, -1.3, -1.1, -0.9, -0.7, -0.5, -0.3, -0.1, - 0.1, 0.3, 0.5, 0.7, 0.9, 1.1, 1.3, 1.5}, - "Sum"}; + // etaAxis (physical) and etaBinAxis (integer index) follow the observable + // binning and are rebuilt in init(). Placeholders here. + AxisSpec etaAxis{9, -0.9, 0.9, "#eta"}; + AxisSpec etaBinAxis{10, -0.5, 9.5, "#eta bin Number"}; - Configurable cfgRunMCGetNSig{"cfgRunMCGetNSig", false, "Run MC pass to get mean of Nsig Plots"}; + AxisSpec gapAxis{{-1.5, -1.3, -1.1, -0.9, -0.7, -0.5, -0.3, -0.1, + 0.1, 0.3, 0.5, 0.7, 0.9, 1.1, 1.3, 1.5}, + "Gap"}; + AxisSpec sumAxis{{-1.5, -1.3, -1.1, -0.9, -0.7, -0.5, -0.3, -0.1, + 0.1, 0.3, 0.5, 0.7, 0.9, 1.1, 1.3, 1.5}, + "Sum"}; + + // --- process switches -------------------------------------------------------- Configurable cfgRunGetEff{"cfgRunGetEff", false, "Run MC pass to build efficiency/fake maps"}; - Configurable cfgRunGetMCFlat{"cfgRunGetMCFlat", false, "Run MC to Get Flattening Weights"}; - Configurable cfgRunMCMean{"cfgRunMCMean", false, "Run MC mean(pT) & mean(Et)"}; + Configurable cfgRunGetMCFlat{"cfgRunGetMCFlat", false, "Run MC to get flattening weights"}; + Configurable cfgRunMCMean{"cfgRunMCMean", false, "Run MC mean(pT)"}; Configurable cfgRunMCFluc{"cfgRunMCFluc", false, "Run MC fluctuations (C2, subevent)"}; - Configurable cfgRunDataGetNSig{"cfgRunDataGetNSig", false, "Run MC pass to get mean of Nsig Plots"}; - Configurable cfgRunGetDataFlat{"cfgRunGetDataFlat", false, "Run Data Get Flattening Weights"}; - Configurable cfgRunDataMean{"cfgRunDataMean", false, "Run DATA mean(pT) & mean(Et)"}; + Configurable cfgRunGetDataFlat{"cfgRunGetDataFlat", false, "Run data get flattening weights"}; + Configurable cfgRunDataMean{"cfgRunDataMean", false, "Run DATA mean(pT)"}; Configurable cfgRunDataFluc{"cfgRunDataFluc", false, "Run DATA fluctuations (C2, subevent)"}; - Service ccdb; - Service pdg; + Service ccdb{}; + Service pdg{}; HistogramRegistry histos{"Histos", {}, OutputObjHandlingPolicy::AnalysisObject}; + TRandom3 rng; // bootstrap Poisson weights + + // --- persistent state -------------------------------------------------------- struct InternalState { - std::array hEff{}; - std::array hFake{}; - std::array hFlatWeight{}; + TH3F* hEff = nullptr; + TH3F* hFake = nullptr; + THnSparseF* hFlatWeight = nullptr; std::vector> mLimitsNchCent; float mMinXNchCent = 0, mMaxXNchCent = 0; - TProfile3D* pmeanTruNchEtabinSpbinStep2 = nullptr; - TProfile3D* pmeanRecoNchEtabinSpbinStep2 = nullptr; - TProfile3D* pmeanRecoEffcorrNchEtabinSpbinStep2 = nullptr; + // MC mean maps (per eta bin): truth / reco / reco-eff-corrected + TProfile2D* pmeanTruNchEtabinStep2 = nullptr; + TProfile2D* pmeanRecoNchEtabinStep2 = nullptr; + TProfile2D* pmeanRecoEffcorrNchEtabinStep2 = nullptr; + TProfile2D* pmeanMultTruNchEtabinStep2 = nullptr; + TProfile2D* pmeanMultRecoNchEtabinStep2 = nullptr; + TProfile2D* pmeanMultRecoEffcorrNchEtabinStep2 = nullptr; - TProfile3D* pmeanMultTruNchEtabinSpbinStep2 = nullptr; - TProfile3D* pmeanMultRecoNchEtabinSpbinStep2 = nullptr; - TProfile3D* pmeanMultRecoEffcorrNchEtabinSpbinStep2 = nullptr; - - TProfile3D* pmeanNchEtabinSpbinStep2 = nullptr; - TProfile3D* pmeanMultNchEtabinSpbinStep2 = nullptr; + // Data mean maps (per eta bin) + TProfile2D* pmeanNchEtabinStep2 = nullptr; + TProfile2D* pmeanMultNchEtabinStep2 = nullptr; TProfile* pmeanFT0AmultpvStep2 = nullptr; TProfile* pmeanFT0CmultpvStep2 = nullptr; } state; o2::ft0::Geometry ft0Det; - template - void fillNSigmaBefCut(const T& track, float cent) - { - float pt = track.pt(); - auto sign = track.sign(); - - if (sign > 0) { - histos.fill(HIST("h3DnsigmaTpcVsPtBefCut_Cent_PiPlus"), cent, pt, track.tpcNSigmaPi()); - histos.fill(HIST("h3DnsigmaTofVsPtBefCut_Cent_PiPlus"), cent, pt, track.tofNSigmaPi()); - histos.fill(HIST("h3DnsigmaTpcVsTofBefCut_Cent_PiPlus"), cent, track.tofNSigmaPi(), track.tpcNSigmaPi()); - - histos.fill(HIST("h3DnsigmaTpcVsPtBefCut_Cent_KaPlus"), cent, pt, track.tpcNSigmaKa()); - histos.fill(HIST("h3DnsigmaTofVsPtBefCut_Cent_KaPlus"), cent, pt, track.tofNSigmaKa()); - histos.fill(HIST("h3DnsigmaTpcVsTofBefCut_Cent_KaPlus"), cent, track.tofNSigmaKa(), track.tpcNSigmaKa()); - - histos.fill(HIST("h3DnsigmaTpcVsPtBefCut_Cent_Pr"), cent, pt, track.tpcNSigmaPr()); - histos.fill(HIST("h3DnsigmaTofVsPtBefCut_Cent_Pr"), cent, pt, track.tofNSigmaPr()); - histos.fill(HIST("h3DnsigmaTpcVsTofBefCut_Cent_Pr"), cent, track.tofNSigmaPr(), track.tpcNSigmaPr()); - } else if (sign < 0) { - histos.fill(HIST("h3DnsigmaTpcVsPtBefCut_Cent_PiMinus"), cent, pt, track.tpcNSigmaPi()); - histos.fill(HIST("h3DnsigmaTofVsPtBefCut_Cent_PiMinus"), cent, pt, track.tofNSigmaPi()); - histos.fill(HIST("h3DnsigmaTpcVsTofBefCut_Cent_PiMinus"), cent, track.tofNSigmaPi(), track.tpcNSigmaPi()); - - histos.fill(HIST("h3DnsigmaTpcVsPtBefCut_Cent_KaMinus"), cent, pt, track.tpcNSigmaKa()); - histos.fill(HIST("h3DnsigmaTofVsPtBefCut_Cent_KaMinus"), cent, pt, track.tofNSigmaKa()); - histos.fill(HIST("h3DnsigmaTpcVsTofBefCut_Cent_KaMinus"), cent, track.tofNSigmaKa(), track.tpcNSigmaKa()); - - histos.fill(HIST("h3DnsigmaTpcVsPtBefCut_Cent_AntiPr"), cent, pt, track.tpcNSigmaPr()); - histos.fill(HIST("h3DnsigmaTofVsPtBefCut_Cent_AntiPr"), cent, pt, track.tofNSigmaPr()); - histos.fill(HIST("h3DnsigmaTpcVsTofBefCut_Cent_AntiPr"), cent, track.tofNSigmaPr(), track.tpcNSigmaPr()); - } - histos.fill(HIST("h3DnsigmaTpcVsPtBefCut_Cent_PiAll"), cent, pt, track.tpcNSigmaPi()); - histos.fill(HIST("h3DnsigmaTofVsPtBefCut_Cent_PiAll"), cent, pt, track.tofNSigmaPi()); - histos.fill(HIST("h3DnsigmaTpcVsTofBefCut_Cent_PiAll"), cent, track.tofNSigmaPi(), track.tpcNSigmaPi()); - - histos.fill(HIST("h3DnsigmaTpcVsPtBefCut_Cent_KaAll"), cent, pt, track.tpcNSigmaKa()); - histos.fill(HIST("h3DnsigmaTofVsPtBefCut_Cent_KaAll"), cent, pt, track.tofNSigmaKa()); - histos.fill(HIST("h3DnsigmaTpcVsTofBefCut_Cent_KaAll"), cent, track.tofNSigmaKa(), track.tpcNSigmaKa()); - - histos.fill(HIST("h3DnsigmaTpcVsPtBefCut_Cent_PrAll"), cent, pt, track.tpcNSigmaPr()); - histos.fill(HIST("h3DnsigmaTofVsPtBefCut_Cent_PrAll"), cent, pt, track.tofNSigmaPr()); - histos.fill(HIST("h3DnsigmaTpcVsTofBefCut_Cent_PrAll"), cent, track.tofNSigmaPr(), track.tpcNSigmaPr()); - } - - template - void fillNSigmaAftCut(const T& track, float cent, bool isSpecies[]) - { - float pt = track.pt(); - float tpcPi = track.tpcNSigmaPi(); - float tofPi = track.tofNSigmaPi(); - - float tpcKa = track.tpcNSigmaKa(); - float tofKa = track.tofNSigmaKa(); - - float tpcPr = track.tpcNSigmaPr(); - float tofPr = track.tofNSigmaPr(); - - if (isSpecies[kPiPlusIdx]) { - histos.fill(HIST("h3DnsigmaTpcVsPtAftCut_Cent_PiPlus"), cent, pt, tpcPi); - histos.fill(HIST("h3DnsigmaTofVsPtAftCut_Cent_PiPlus"), cent, pt, tofPi); - histos.fill(HIST("h3DnsigmaTpcVsTofAftCut_Cent_PiPlus"), cent, tofPi, tpcPi); - } - if (isSpecies[kPiMinusIdx]) { - histos.fill(HIST("h3DnsigmaTpcVsPtAftCut_Cent_PiMinus"), cent, pt, tpcPi); - histos.fill(HIST("h3DnsigmaTofVsPtAftCut_Cent_PiMinus"), cent, pt, tofPi); - histos.fill(HIST("h3DnsigmaTpcVsTofAftCut_Cent_PiMinus"), cent, tofPi, tpcPi); - } - if (isSpecies[kPiAllIdx]) { - histos.fill(HIST("h3DnsigmaTpcVsPtAftCut_Cent_PiAll"), cent, pt, tpcPi); - histos.fill(HIST("h3DnsigmaTofVsPtAftCut_Cent_PiAll"), cent, pt, tofPi); - histos.fill(HIST("h3DnsigmaTpcVsTofAftCut_Cent_PiAll"), cent, tofPi, tpcPi); - } - if (isSpecies[kKaPlusIdx]) { - histos.fill(HIST("h3DnsigmaTpcVsPtAftCut_Cent_KaPlus"), cent, pt, tpcKa); - histos.fill(HIST("h3DnsigmaTofVsPtAftCut_Cent_KaPlus"), cent, pt, tofKa); - histos.fill(HIST("h3DnsigmaTpcVsTofAftCut_Cent_KaPlus"), cent, tofKa, tpcKa); - } - if (isSpecies[kKaMinusIdx]) { - histos.fill(HIST("h3DnsigmaTpcVsPtAftCut_Cent_KaMinus"), cent, pt, tpcKa); - histos.fill(HIST("h3DnsigmaTofVsPtAftCut_Cent_KaMinus"), cent, pt, tofKa); - histos.fill(HIST("h3DnsigmaTpcVsTofAftCut_Cent_KaMinus"), cent, tofKa, tpcKa); - } - if (isSpecies[kKaAllIdx]) { - histos.fill(HIST("h3DnsigmaTpcVsPtAftCut_Cent_KaAll"), cent, pt, tpcKa); - histos.fill(HIST("h3DnsigmaTofVsPtAftCut_Cent_KaAll"), cent, pt, tofKa); - histos.fill(HIST("h3DnsigmaTpcVsTofAftCut_Cent_KaAll"), cent, tofKa, tpcKa); - } - if (isSpecies[kPrIdx]) { - histos.fill(HIST("h3DnsigmaTpcVsPtAftCut_Cent_Pr"), cent, pt, tpcPr); - histos.fill(HIST("h3DnsigmaTofVsPtAftCut_Cent_Pr"), cent, pt, tofPr); - histos.fill(HIST("h3DnsigmaTpcVsTofAftCut_Cent_Pr"), cent, tofPr, tpcPr); - } - if (isSpecies[kAntiPrIdx]) { - histos.fill(HIST("h3DnsigmaTpcVsPtAftCut_Cent_AntiPr"), cent, pt, tpcPr); - histos.fill(HIST("h3DnsigmaTofVsPtAftCut_Cent_AntiPr"), cent, pt, tofPr); - histos.fill(HIST("h3DnsigmaTpcVsTofAftCut_Cent_AntiPr"), cent, tofPr, tpcPr); - } - if (isSpecies[kPrAllIdx]) { - histos.fill(HIST("h3DnsigmaTpcVsPtAftCut_Cent_PrAll"), cent, pt, tpcPr); - histos.fill(HIST("h3DnsigmaTofVsPtAftCut_Cent_PrAll"), cent, pt, tofPr); - histos.fill(HIST("h3DnsigmaTpcVsTofAftCut_Cent_PrAll"), cent, tofPr, tpcPr); - } else { - return; - } - } - - // Returns: 0 = Unknown/Reject, 1 = Pion, 2 = Kaon, 3 = Proton - template - int identifyTrack(const T& candidate, int cent) - { - if (!candidate.hasTPC()) - return 0; - - float pt = candidate.pt(); - if (pt <= cfgCutPtLower || pt >= cfgCutPtUpperPID) - return 0; // Out of bounds - - if (!pidMeanSigmaMap) - return 0; - - int centBin = cent + 0.5; - auto charge = candidate.sign(); - int piIdx = (charge > 0) ? kPiPlusIdx : kPiMinusIdx; - int kaIdx = (charge > 0) ? kKaPlusIdx : kKaMinusIdx; - int prIdx = (charge > 0) ? kPrIdx : kAntiPrIdx; - - // TPC - float mPiTpc = pidMeanSigmaMap->meanTPC[piIdx][centBin]; - float sPiTpc = pidMeanSigmaMap->sigmaTPC[piIdx][centBin]; - - float mKaTpc = pidMeanSigmaMap->meanTPC[kaIdx][centBin]; - float sKaTpc = pidMeanSigmaMap->sigmaTPC[kaIdx][centBin]; - - float mPrTpc = pidMeanSigmaMap->meanTPC[prIdx][centBin]; - float sPrTpc = pidMeanSigmaMap->sigmaTPC[prIdx][centBin]; - - // TOF - float mPiTof = pidMeanSigmaMap->meanTOF[piIdx][centBin]; - float sPiTof = pidMeanSigmaMap->sigmaTOF[piIdx][centBin]; - - float mKaTof = pidMeanSigmaMap->meanTOF[kaIdx][centBin]; - float sKaTof = pidMeanSigmaMap->sigmaTOF[kaIdx][centBin]; - - float mPrTof = pidMeanSigmaMap->meanTOF[prIdx][centBin]; - float sPrTof = pidMeanSigmaMap->sigmaTOF[prIdx][centBin]; - - static int debugLogCounter = 0; - if (debugLogCounter < KConstTen) { - LOGF(info, "[PID DEBUG] CentBin: %d, Charge: %d", centBin, charge); - LOGF(info, " -> TPC USED | Pi (\u03bc=%.3f, \u03c3=%.3f) | Ka (\u03bc=%.3f, \u03c3=%.3f) | Pr (\u03bc=%.3f, \u03c3=%.3f)", - mPiTpc, sPiTpc, mKaTpc, sKaTpc, mPrTpc, sPrTpc); - - if (candidate.hasTOF()) { - LOGF(info, " -> TOF USED | Pi (\u03bc=%.3f, \u03c3=%.3f) | Ka (\u03bc=%.3f, \u03c3=%.3f) | Pr (\u03bc=%.3f, \u03c3=%.3f)", - mPiTof, sPiTof, mKaTof, sKaTof, mPrTof, sPrTof); - } else { - LOGF(info, " -> TOF USED | Track has no TOF signal."); - } - debugLogCounter++; - } - // Fetch Raw nSigma Values - float rawTpcPi = candidate.tpcNSigmaPi(); - float rawTpcKa = candidate.tpcNSigmaKa(); - float rawTpcPr = candidate.tpcNSigmaPr(); - - float rawTofPi = 0.f, rawTofKa = 0.f, rawTofPr = 0.f; - if (candidate.hasTOF()) { - rawTofPi = candidate.tofNSigmaPi(); - rawTofKa = candidate.tofNSigmaKa(); - rawTofPr = candidate.tofNSigmaPr(); - } - - // --- Low PT Regime --- - if (pt <= cfgCutPtUpperTPC) { - // Basic TPC passing check: |Raw - Mean| < (Cut * Sigma) - bool inTpcPi = std::abs(rawTpcPi - mPiTpc) < (cfgnSigmaCutTPC * sPiTpc); - bool inTpcKa = std::abs(rawTpcKa - mKaTpc) < (cfgnSigmaCutTPC * sKaTpc); - bool inTpcPr = std::abs(rawTpcPr - mPrTpc) < (cfgnSigmaCutTPC * sPrTpc); - - // Combined passing check (adds TOF if available) - bool passPi = inTpcPi && (!candidate.hasTOF() || std::abs(rawTofPi - mPiTof) < (cfgnSigmaCutTOF * sPiTof)); - bool passKa = inTpcKa && (!candidate.hasTOF() || std::abs(rawTofKa - mKaTof) < (cfgnSigmaCutTOF * sKaTof)); - bool passPr = inTpcPr && (!candidate.hasTOF() || std::abs(rawTofPr - mPrTof) < (cfgnSigmaCutTOF * sPrTof)); - - // Uniqueness check: Must pass target cut, and NOT fall into the TPC range of the others - if (passPi && !passKa && !passPr) - return 1; - if (passKa && !passPi && !passPr) - return 2; - if (passPr && !passPi && !passKa) - return 3; - - return 0; // Ambiguous or failed all cuts - } - - // --- High PT Regime--- - if (candidate.hasTOF() && pt > cfgCutPtUpperTPC) { - // Calculate 2D Normalized Distance (Elliptical distance normalized by sigma) - float dPi = std::hypot((rawTpcPi - mPiTpc) / sPiTpc, (rawTofPi - mPiTof) / sPiTof); - float dKa = std::hypot((rawTpcKa - mKaTpc) / sKaTpc, (rawTofKa - mKaTof) / sKaTof); - float dPr = std::hypot((rawTpcPr - mPrTpc) / sPrTpc, (rawTofPr - mPrTof) / sPrTof); - - // Count how many particles are within the ambiguity radius - int competitors = (dPi < cfgnSigmaOtherParticles) + - (dKa < cfgnSigmaOtherParticles) + - (dPr < cfgnSigmaOtherParticles); - - // If 1 or fewer are in the ambiguity region, pick the absolute best match - if (competitors <= 1) { - if (dPi <= dKa && dPi <= dPr && dPi < cfgnSigmaCutCombTPCTOF) - return 1; - if (dKa <= dPi && dKa <= dPr && dKa < cfgnSigmaCutCombTPCTOF) - return 2; - if (dPr <= dPi && dPr <= dKa && dPr < cfgnSigmaCutCombTPCTOF) - return 3; - } - } - return 0; // Unknown/Reject - } - + // --- bootstrap replica storage (base data fluc only) ------------------------- + struct BootstrapHists { + std::array, KMaxBoot> meanpTCent{}; + std::array, KMaxBoot> meanpTMult{}; + std::array, KMaxBoot> c2Cent{}; + std::array, KMaxBoot> c2Mult{}; + std::array, KMaxBoot> c2SubCent{}; + std::array, KMaxBoot> c2SubMult{}; + std::array, KMaxBoot> covCent{}; + std::array, KMaxBoot> covMult{}; + std::array, KMaxBoot> covFT0ACent{}; + std::array, KMaxBoot> covFT0AMult{}; + std::array, KMaxBoot> covFT0CCent{}; + std::array, KMaxBoot> covFT0CMult{}; + std::array, KMaxBoot> c2Sub2D{}; + std::array, KMaxBoot> gapSum2D{}; + std::array, KMaxBoot> cov2D{}; + std::array, KMaxBoot> covFT0A2D{}; + std::array, KMaxBoot> covFT0C2D{}; + } bs; + + // =========================================================================== + // Selection helpers + // =========================================================================== template bool isEventSelected(const T& col) { histos.fill(HIST("hEvtCount"), 0.5); - - if (!col.sel8()) + if (!col.sel8()) { return false; + } histos.fill(HIST("hEvtCount"), 1.5); - - if (std::abs(col.posZ()) > cfgCutVertex) + if (std::abs(col.posZ()) > cfgVtxZCut) { return false; + } histos.fill(HIST("hEvtCount"), 2.5); - - if (cfgEvSelNoSameBunchPileup && !col.selection_bit(o2::aod::evsel::kNoSameBunchPileup)) + if (cfgEvSelNoSameBunchPileup && !col.selection_bit(o2::aod::evsel::kNoSameBunchPileup)) { return false; + } histos.fill(HIST("hEvtCount"), 3.5); - - if (cfgUseGoodITSLayerAllCut && !col.selection_bit(o2::aod::evsel::kIsGoodITSLayersAll)) + if (cfgUseGoodITSLayerAllCut && !col.selection_bit(o2::aod::evsel::kIsGoodITSLayersAll)) { return false; + } histos.fill(HIST("hEvtCount"), 4.5); - - if (cfgIsGoodZvtxFT0VsPV && !col.selection_bit(o2::aod::evsel::kIsGoodZvtxFT0vsPV)) + if (cfgIsGoodZvtxFT0VsPV && !col.selection_bit(o2::aod::evsel::kIsGoodZvtxFT0vsPV)) { return false; + } histos.fill(HIST("hEvtCount"), 5.5); return true; } - bool isPassAddPileup(int multPV, int trksize, float cent) + bool isPassAddPileup(float multPV, int trksize, float cent) { auto checkLimits = [](float x, float y, const std::vector>& limits, float xM, float xMx) { - if (limits.empty()) + if (limits.empty()) { return true; + } int bin = 1 + static_cast((x - xM) / (xMx - xM) * (limits.size() - 2)); - if (bin < 1 || bin >= static_cast(limits.size() - 1)) + if (bin < 1 || bin >= static_cast(limits.size() - 1)) { return false; + } return (y >= limits[bin].first && y <= limits[bin].second); }; if (cfgApplySigPupCut) { - if (!checkLimits(cent, trksize, state.mLimitsNchCent, state.mMinXNchCent, state.mMaxXNchCent)) + if (!checkLimits(cent, trksize, state.mLimitsNchCent, state.mMinXNchCent, state.mMaxXNchCent)) { return false; + } histos.fill(HIST("hEvtCount"), 6.5); } - if (cfgApplyLinPupCut) { - if (trksize > (cfgLinPupParam0 + cfgLinPupParam1 * multPV)) + if (trksize > (cfgLinPupParam0 + cfgLinPupParam1 * multPV)) { return false; + } histos.fill(HIST("hEvtCount"), 7.5); - if (trksize < (cfgLinPupParam2 + cfgLinPupParam3 * multPV)) + if (trksize < (cfgLinPupParam2 + cfgLinPupParam3 * multPV)) { return false; + } histos.fill(HIST("hEvtCount"), 8.5); } return true; } - template - bool hasMinTracksInAllEtaBins(const double (&sw)[NspT][NetaT][NkT]) + // Minimum weighted-track requirement in every narrow eta bin (systNEta). Loops + // only over the active bins [1, nEta); the array is sized to KNEtaMax. + template + bool hasMinTracksInAllEtaBins(const std::array, NetaT>& sw) { const int minTracks = cfgMinTracksPerEtaBin; - if (minTracks <= 0) + if (minTracks <= 0) { return true; - for (std::size_t ieta = 1; ieta < NetaT; ++ieta) { - if (sw[kInclusiveIdx][ieta][1] < static_cast(minTracks)) + } + for (int ieta = 1; ieta < nEta; ++ieta) { + if (sw[ieta][1] < static_cast(minTracks)) { return false; + } } return true; } - - template - bool hasMinTracksInAllEtaBins(const double (&sw)[NspT][NetaT]) + template + bool hasMinTracksInAllEtaBins(const std::array& sw) { const int minTracks = cfgMinTracksPerEtaBin; - if (minTracks <= 0) + if (minTracks <= 0) { return true; - for (std::size_t ieta = 1; ieta < NetaT; ++ieta) { - if (sw[kInclusiveIdx][ieta] < static_cast(minTracks)) + } + for (int ieta = 1; ieta < nEta; ++ieta) { + if (sw[ieta] < static_cast(minTracks)) { return false; + } } return true; } @@ -579,41 +386,40 @@ struct RadialFlowDecorr { bool isTrackSelected(const T& trk) { histos.fill(HIST("hTrkCount"), 0.5); - - if (trk.sign() == 0) + if (trk.sign() == 0) { return false; + } histos.fill(HIST("hTrkCount"), 1.5); - - if (!trk.has_collision()) + if (!trk.has_collision()) { return false; + } histos.fill(HIST("hTrkCount"), 2.5); - - if (!trk.isPVContributor()) + if (!trk.isPVContributor()) { return false; + } histos.fill(HIST("hTrkCount"), 3.5); - - if (!(trk.itsNCls() > cfgITScluster)) + if (!(trk.itsNCls() > cfgITScluster)) { return false; + } histos.fill(HIST("hTrkCount"), 4.5); - - if (!(trk.tpcNClsFound() >= cfgTPCcluster)) + if (!(trk.tpcNClsFound() >= cfgTPCcluster)) { return false; + } histos.fill(HIST("hTrkCount"), 5.5); - - if (!(trk.tpcNClsCrossedRows() >= cfgTPCnCrossedRows)) + if (!(trk.tpcNClsCrossedRows() >= cfgTPCnCrossedRows)) { return false; + } histos.fill(HIST("hTrkCount"), 6.5); - - if (trk.pt() < cfgCutPtLower || trk.pt() > cfgCutPtUpper || std::abs(trk.eta()) > cfgCutEta) + if (trk.pt() < cfgCutPtLower || trk.pt() > cfgCutPtUpper || std::abs(trk.eta()) > cfgCutEta) { return false; + } histos.fill(HIST("hTrkCount"), 7.5); - - if (!trk.isGlobalTrack()) + if (!trk.isGlobalTrack()) { return false; + } histos.fill(HIST("hTrkCount"), 8.5); if (cfgPtDepDCAxy) { - // Evaluates: P0 + P1 / (pt^P2) float maxDcaXY = cfgDcaXyP0 + cfgDcaXyP1 / std::pow(trk.pt(), cfgDcaXyP2); if (std::abs(trk.dcaXY()) > maxDcaXY) { return false; @@ -626,10 +432,9 @@ struct RadialFlowDecorr { histos.fill(HIST("hTrkCount"), 9.5); } if (cfgPtDepDCAz) { - // Evaluates: P0 + P1 / (pt^P2) float maxDcaZ = cfgDcaZP0 + cfgDcaZP1 / std::pow(trk.pt(), cfgDcaZP2); if (std::abs(trk.dcaZ()) > maxDcaZ) { - return false; // Reject track if DCA is too large + return false; } histos.fill(HIST("hTrkCount"), 10.5); } else { @@ -645,67 +450,74 @@ struct RadialFlowDecorr { bool isParticleSelected(const T& particle) { auto* pd = pdg->GetParticle(particle.pdgCode()); - if (!pd) + if (!pd) { return false; - if (std::abs(pd->Charge()) == 0) + } + if (std::abs(pd->Charge()) == 0) { return false; - if (particle.pt() < cfgCutPtLower || particle.pt() > cfgCutPtUpper || std::abs(particle.eta()) > cfgCutEta) + } + if (particle.pt() < cfgCutPtLower || particle.pt() > cfgCutPtUpper || std::abs(particle.eta()) > cfgCutEta) { return false; - if (std::abs(particle.vz()) > cfgCutVertex) + } + if (std::abs(particle.vz()) > cfgVtxZCut) { return false; + } return true; } float getCentrality(const auto& col) const { - if (cfgCentralityChoice.value == kCentFT0C) + if (cfgCentralityChoice.value == kCentFT0C) { return col.centFT0C(); - if (cfgCentralityChoice.value == kCentFT0M) + } + if (cfgCentralityChoice.value == kCentFT0M) { return col.centFT0M(); - if (cfgCentralityChoice.value == kCentFDDM) + } + if (cfgCentralityChoice.value == kCentFDDM) { return col.centFDDM(); - if (cfgCentralityChoice.value == kCentFV0A) + } + if (cfgCentralityChoice.value == kCentFV0A) { return col.centFV0A(); + } return KinvalidCentrality; } - float getEfficiency(float mult, float pt, float eta, PIDIdx pidType, int effidx, bool cfgEff) const + // Inclusive efficiency/fake lookup (no species dependence). + float getEfficiency(float mult, float pt, float eta, int effidx, bool useEff) const { - if (!cfgEff) { + if (!useEff) { return (effidx == 0) ? 1.0f : 0.0f; } - TH3F* h = (effidx == 0) ? state.hEff[pidType] : state.hFake[pidType]; - - if (!h) + TH3F* h = (effidx == 0) ? state.hEff : state.hFake; + if (!h) { return -1; - + } int ibx = h->GetXaxis()->FindBin(mult); int iby = h->GetYaxis()->FindBin(pt); int ibz = h->GetZaxis()->FindBin(eta); float val = h->GetBinContent(ibx, iby, ibz); - - if (effidx == 0) + if (effidx == 0) { return (val > 0.f) ? val : 1.0f; + } return val; } - float getFlatteningWeight(float vz, float chg, float pt, float eta, float phi, PIDIdx pidType, bool cfgflat) const + float getFlatteningWeight(float vz, float chg, float pt, float eta, float phi, bool useFlat) const { - if (!cfgflat) + if (!useFlat) { return 1.0; - THnSparseF* h = state.hFlatWeight[pidType]; - - if (!h) + } + THnSparseF* h = state.hFlatWeight; + if (!h) { return 0.0; - int bins[5]; + } + std::array bins{}; bins[0] = h->GetAxis(0)->FindBin(vz); bins[1] = h->GetAxis(1)->FindBin(chg); bins[2] = h->GetAxis(2)->FindBin(pt); bins[3] = h->GetAxis(3)->FindBin(eta); bins[4] = h->GetAxis(4)->FindBin(phi); - float val = h->GetBinContent(bins); - - return val; + return h->GetBinContent(bins.data()); } std::vector* offsetFT0 = nullptr; @@ -719,8 +531,6 @@ struct RadialFlowDecorr { auto x = chPos.X() + (*offsetFT0)[i].getX(); auto y = chPos.Y() + (*offsetFT0)[i].getY(); auto z = chPos.Z() + (*offsetFT0)[i].getZ(); - // i == 0 is FT0A (A-side, positive Z) - // i == 1 is FT0C (C-side, negative Z) if (i == 1) { z = -std::abs(z); } else if (i == 0) { @@ -742,13 +552,12 @@ struct RadialFlowDecorr { return; } mLastTimestamp = timestamp; - LOGF(info, "Successfully loaded new alignment parameters for timestamp %llu", timestamp); - LOGF(info, "Offset for FT0A: x = %.3f y = %.3f z = %.3f\n", (*offsetFT0)[0].getX(), (*offsetFT0)[0].getY(), (*offsetFT0)[0].getZ()); - LOGF(info, "Offset for FT0C: x = %.3f y = %.3f z = %.3f\n", (*offsetFT0)[1].getX(), (*offsetFT0)[1].getY(), (*offsetFT0)[1].getZ()); + LOGF(info, "Loaded FT0 alignment for timestamp %llu", timestamp); } - template - std::pair calculateMeanAndC2FromSums(const double sumpmwk[KIntM][KIntK], const double sumwk[KIntK], float referenceMeanPt) const + // Two-particle pT correlator from the per-event power sums. + template + std::pair calculateMeanAndC2FromSums(const std::array, M>& sumpmwk, const std::array& sumwk, float referenceMeanPt) const { if (sumwk[1] == 0.) { return {0.f, 0.f}; @@ -763,16 +572,8 @@ struct RadialFlowDecorr { } double pmk11 = sumpmwk[1][1] / sumwk[1]; - - double pmk12 = 0.f; - if (sumwk[2] != 0.f) { - pmk12 = sumpmwk[1][2] / sumwk[2]; - } - - double pmk22 = 0.f; - if (sumwk[2] != 0.f) { - pmk22 = sumpmwk[2][2] / sumwk[2]; - } + double pmk12 = (sumwk[2] != 0.f) ? sumpmwk[1][2] / sumwk[2] : 0.f; + double pmk22 = (sumwk[2] != 0.f) ? sumpmwk[2][2] / sumwk[2] : 0.f; float calculatedMeanPt = pmk11; @@ -786,6 +587,9 @@ struct RadialFlowDecorr { return {calculatedMeanPt, twopcorr}; } + // =========================================================================== + // Table joins + // =========================================================================== using GeneralCollisions = soa::Join; using UnfilteredTracks = soa::Join< - aod::Tracks, - aod::TracksExtra, - aod::TrackSelection, - aod::TracksDCA, - aod::pidTPCFullPi, aod::pidTPCFullKa, aod::pidTPCFullPr, aod::pidTPCFullEl, - aod::pidTOFFullPi, aod::pidTOFFullKa, aod::pidTOFFullPr, aod::pidTOFFullEl>; - Filter trackFilter = aod::track::pt > KPtMin&& - aod::track::pt < KPtMax&& - requireGlobalTrackInFilter(); + aod::Tracks, aod::TracksExtra, aod::TrackSelection, aod::TracksDCA>; + Filter trackFilter = aod::track::pt > KPtMin&& aod::track::pt < KPtMax&& requireGlobalTrackInFilter(); using AodTracksSel = soa::Filtered; using TCs = soa::Join; using FilteredTCs = soa::Filtered; @@ -816,6 +613,9 @@ struct RadialFlowDecorr { PresliceUnsorted colPerMcCollision = aod::mccollisionlabel::mcCollisionId; + // =========================================================================== + // Histogram declarations + // =========================================================================== void declareCommonQA() { histos.add("hVtxZ_after_sel", ";z_{vtx} (cm)", kTH1F, {{KNbinsZvtx, KZvtxMin, KZvtxMax}}); @@ -823,7 +623,6 @@ struct RadialFlowDecorr { histos.add("hCentrality", ";centrality (%)", kTH1F, {{centAxis1Per}}); histos.add("Hist2D_globalTracks_PVTracks", ";N_{global};N_{PV}", kTH2F, {{nChAxis}, {nChAxis}}); histos.add("Hist2D_cent_nch", ";N_{PV};cent (%)", kTH2F, {{nChAxis}, {centAxis1Per}}); - histos.add("Hist2D_globalTracks_cent", "cent (%);N_{global}", kTH2F, {{centAxis1Per}, {nChAxis}}); histos.add("Hist2D_PVTracks_cent", "cent (%);N_{PV}", kTH2F, {{centAxis1Per}, {nChAxis}}); @@ -859,14 +658,11 @@ struct RadialFlowDecorr { void declareMCCommonHists() { - for (const auto& suf : pidSuffix) { - histos.add("h3_AllPrimary" + suf, ";N_{PV};p_{T};#eta", kTH3F, {{nChAxis2}, {KNbinsPtRes, KPtMin, KPtMax}, {KNbinsEtaFine, -KEtaFineMax, KEtaFineMax}}); - histos.add("h3_RecoMatchedToPrimary" + suf, ";N_{PV};p_{T};#eta", kTH3F, {{nChAxis2}, {KNbinsPtRes, KPtMin, KPtMax}, {KNbinsEtaFine, -KEtaFineMax, KEtaFineMax}}); - histos.add("h3_AllReco" + suf, ";N_{PV};p_{T};#eta", kTH3F, {{nChAxis2}, {KNbinsPtRes, KPtMin, KPtMax}, {KNbinsEtaFine, -KEtaFineMax, KEtaFineMax}}); - histos.add("h3_RecoUnMatchedToPrimary_Secondary" + suf, ";N_{PV};p_{T};#eta", kTH3F, {{nChAxis2}, {KNbinsPtRes, KPtMin, KPtMax}, {KNbinsEtaFine, -KEtaFineMax, KEtaFineMax}}); - histos.add("h3_RecoUnMatchedToPrimary_Fake" + suf, ";N_{PV};p_{T};#eta", kTH3F, {{nChAxis2}, {KNbinsPtRes, KPtMin, KPtMax}, {KNbinsEtaFine, -KEtaFineMax, KEtaFineMax}}); - histos.add("h3_RecoMatchedToPrimary_MisID" + suf, ";N_{PV};p_{T};#eta", kTH3F, {{nChAxis2}, {KNbinsPtRes, KPtMin, KPtMax}, {KNbinsEtaFine, -KEtaFineMax, KEtaFineMax}}); - } + histos.add("h3_AllPrimary", ";N_{PV};p_{T};#eta", kTH3F, {{nChAxis2}, {KNbinsPtRes, KPtMin, KPtMax}, {KNbinsEtaFine, -KEtaFineMax, KEtaFineMax}}); + histos.add("h3_RecoMatchedToPrimary", ";N_{PV};p_{T};#eta", kTH3F, {{nChAxis2}, {KNbinsPtRes, KPtMin, KPtMax}, {KNbinsEtaFine, -KEtaFineMax, KEtaFineMax}}); + histos.add("h3_AllReco", ";N_{PV};p_{T};#eta", kTH3F, {{nChAxis2}, {KNbinsPtRes, KPtMin, KPtMax}, {KNbinsEtaFine, -KEtaFineMax, KEtaFineMax}}); + histos.add("h3_RecoUnMatchedToPrimary_Secondary", ";N_{PV};p_{T};#eta", kTH3F, {{nChAxis2}, {KNbinsPtRes, KPtMin, KPtMax}, {KNbinsEtaFine, -KEtaFineMax, KEtaFineMax}}); + histos.add("h3_RecoUnMatchedToPrimary_Fake", ";N_{PV};p_{T};#eta", kTH3F, {{nChAxis2}, {KNbinsPtRes, KPtMin, KPtMax}, {KNbinsEtaFine, -KEtaFineMax, KEtaFineMax}}); histos.add("ptResolution", ";p_{T}^{MC};(p_{T}^{reco}-p_{T}^{MC})/p_{T}^{MC}", kTH2F, {{KNbinsPtRes, KPtMin, KPtMax}, {100, -0.2, 0.2}}); histos.add("etaResolution", ";#eta^{MC};#eta^{reco}-#eta^{MC}", kTH2F, {{KNbinsEtaRes, -KEtaFineMax, KEtaFineMax}, {100, -0.02, 0.02}}); @@ -875,25 +671,11 @@ struct RadialFlowDecorr { histos.add("vzResolution", ";Vz^{MC};(Vz^{reco}-Vz^{MC})/Vz^{MC}", kTH2F, {{KNbinsVz, KVzMin, KVzMax}, {100, -0.1, 0.1}}); } - void declarenSigHists() - { - for (const auto& suf : pidSuffix) { - histos.add("h3DnsigmaTpcVsPtBefCut_Cent" + suf, "TPC nSigma vs pT Before Cut;cent [%]; p_{T} (GeV/c);n#sigma_{TPC}", kTH3F, {{centAxis1Per}, {KNbinsPtRes, KPtMin, KPtMax}, {200, -10.f, 10.f}}); - histos.add("h3DnsigmaTofVsPtBefCut_Cent" + suf, "TOF nSigma vs pT Before Cut;cent [%]; p_{T} (GeV/c);n#sigma_{TOF}", kTH3F, {{centAxis1Per}, {KNbinsPtRes, KPtMin, KPtMax}, {200, -10.f, 10.f}}); - histos.add("h3DnsigmaTpcVsTofBefCut_Cent" + suf, "TPC vs TOF nSigma Before Cut;cent [%]; n#sigma_{TOF};n#sigma_{TPC}", kTH3F, {{centAxis1Per}, {200, -10.f, 10.f}, {200, -10.f, 10.f}}); - histos.add("h3DnsigmaTpcVsPtAftCut_Cent" + suf, "TPC nSigma vs pT After Cut;cent [%],; p_{T} (GeV/c);n#sigma_{TPC}", kTH3F, {{centAxis1Per}, {KNbinsPtRes, KPtMin, KPtMax}, {200, -10.f, 10.f}}); - histos.add("h3DnsigmaTofVsPtAftCut_Cent" + suf, "TOF nSigma vs pT After Cut;cent [%],; p_{T} (GeV/c);n#sigma_{TOF}", kTH3F, {{centAxis1Per}, {KNbinsPtRes, KPtMin, KPtMax}, {200, -10.f, 10.f}}); - histos.add("h3DnsigmaTpcVsTofAftCut_Cent" + suf, "TPC vs TOF nSigma After Cut;cent [%],; n#sigma_{TOF};n#sigma_{TPC}", kTH3F, {{centAxis1Per}, {200, -10.f, 10.f}, {200, -10.f, 10.f}}); - } - } - void declareMCGetFlatHists() { - for (const auto& suf : pidSuffix) { - histos.add("MCGen/hEtaPhiReco" + suf, ";vz;sign;pt;eta;phi", kTHnSparseF, {{vzAxis}, {chgAxis}, {pTAxis}, {etaAxis}, {phiAxis}}); - histos.add("MCGen/hEtaPhiRecoEffWtd" + suf, ";vz;sign;pt;eta;phi", kTHnSparseF, {{vzAxis}, {chgAxis}, {pTAxis}, {etaAxis}, {phiAxis}}); - histos.add("MCGen/hEtaPhiRecoWtd" + suf, ";vz;sign;pt;eta;phi", kTHnSparseF, {{vzAxis}, {chgAxis}, {pTAxis}, {etaAxis}, {phiAxis}}); - } + histos.add("MCGen/hEtaPhiReco", ";vz;sign;pt;eta;phi", kTHnSparseF, {{vzAxis}, {chgAxis}, {pTAxis}, {etaFlatAxis}, {phiAxis}}); + histos.add("MCGen/hEtaPhiRecoEffWtd", ";vz;sign;pt;eta;phi", kTHnSparseF, {{vzAxis}, {chgAxis}, {pTAxis}, {etaFlatAxis}, {phiAxis}}); + histos.add("MCGen/hEtaPhiRecoWtd", ";vz;sign;pt;eta;phi", kTHnSparseF, {{vzAxis}, {chgAxis}, {pTAxis}, {etaFlatAxis}, {phiAxis}}); } void declareMCMeanHists() @@ -902,12 +684,10 @@ struct RadialFlowDecorr { histos.add("Eff_Ntrk", ";N_{PV}", kTProfile, {nChAxis2}); histos.add("Eff_pT", ";p_{T}", kTProfile, {{KNbinsPtRes, KPtMin, KPtMax}}); histos.add("Eff_eta", ";#eta", kTProfile, {{KNbinsEtaFine, -KEtaFineMax, KEtaFineMax}}); - histos.add("Fake_cent", ";cent", kTProfile, {centAxis1Per}); histos.add("Fake_Ntrk", ";N_{PV}", kTProfile, {nChAxis2}); histos.add("Fake_pT", ";p_{T}", kTProfile, {{KNbinsPtRes, KPtMin, KPtMax}}); histos.add("Fake_eta", ";#eta", kTProfile, {{KNbinsEtaFine, -KEtaFineMax, KEtaFineMax}}); - histos.add("wgt_cent", ";cent", kTProfile, {centAxis1Per}); histos.add("wgt_Ntrk", ";N_{PV}", kTProfile, {nChAxis2}); histos.add("wgt_pT", ";p_{T}", kTProfile, {{KNbinsPtRes, KPtMin, KPtMax}}); @@ -920,71 +700,64 @@ struct RadialFlowDecorr { histos.add("pmean_cent_id_eta_FT0", ";cent;id;#eta", kTProfile3D, {{centAxis1Per}, {200, -0.5, 199.5}, {100, -5.0, 5.0}}); histos.add("h3_cent_id_eta_FT0", ";cent;id;#eta", kTH3F, {{centAxis1Per}, {200, -0.5, 199.5}, {100, -5.0, 5.0}}); - histos.add("MCGen/Prof_Cent_Nsp_Nchrec", ";cent;isp", kTProfile2D, {{centAxis1Per}, {spBinAxis}}); - histos.add("MCGen/Prof_Mult_Nsp_Nchrec", ";mult;isp", kTProfile2D, {{nChAxis}, {spBinAxis}}); - histos.add("MCGen/Prof_Cent_Nsp_MeanpT", ";cent;isp", kTProfile2D, {{centAxis1Per}, {spBinAxis}}); - histos.add("MCGen/Prof_Mult_Nsp_MeanpT", ";mult;isp", kTProfile2D, {{nChAxis}, {spBinAxis}}); - - histos.add("pmeanTru_nch_etabin_spbin", ";mult;eta;isp", kTProfile3D, {{nChAxis}, {etaBinAxis}, {spBinAxis}}); - histos.add("pmeanReco_nch_etabin_spbin", ";mult;eta;isp", kTProfile3D, {{nChAxis}, {etaBinAxis}, {spBinAxis}}); - histos.add("pmeanRecoEffcorr_nch_etabin_spbin", ";mult;eta;isp", kTProfile3D, {{nChAxis}, {etaBinAxis}, {spBinAxis}}); - - histos.add("pmeanMultTru_nch_etabin_spbin", ";mult;eta;isp", kTProfile3D, {{nChAxis}, {etaBinAxis}, {spBinAxis}}); - histos.add("pmeanMultReco_nch_etabin_spbin", ";mult;eta;isp", kTProfile3D, {{nChAxis}, {etaBinAxis}, {spBinAxis}}); - histos.add("pmeanMultRecoEffcorr_nch_etabin_spbin", ";mult;eta;isp", kTProfile3D, {{nChAxis}, {etaBinAxis}, {spBinAxis}}); - - for (const auto& suf : pidSuffix) { - histos.add("hEtaPhiReco" + suf, ";vz;sign;pt;eta;phi", kTHnSparseF, {{vzAxis}, {chgAxis}, {pTAxis}, {etaAxis}, {phiAxis}}); - histos.add("hEtaPhiRecoEffWtd" + suf, ";vz;sign;pt;eta;phi", kTHnSparseF, {{vzAxis}, {chgAxis}, {pTAxis}, {etaAxis}, {phiAxis}}); - histos.add("hEtaPhiRecoWtd" + suf, ";vz;sign;pt;eta;phi", kTHnSparseF, {{vzAxis}, {chgAxis}, {pTAxis}, {etaAxis}, {phiAxis}}); - - histos.add("Prof2D_MeanpTSub_Tru" + suf, ";cent;etaA;etaC", kTProfile3D, {{centAxis1Per}, {etaBinAxis}, {etaBinAxis}}); - histos.add("Prof2D_MeanpTSub_Reco" + suf, ";cent;etaA;etaC", kTProfile3D, {{centAxis1Per}, {etaBinAxis}, {etaBinAxis}}); - histos.add("Prof2D_MeanpTSub_RecoEffCorr" + suf, ";cent;etaA;etaC", kTProfile3D, {{centAxis1Per}, {etaBinAxis}, {etaBinAxis}}); - } + histos.add("MCGen/Prof_Cent_Nchrec", ";cent;#LT N#GT", kTProfile, {centAxis1Per}); + histos.add("MCGen/Prof_Mult_Nchrec", ";mult;#LT N#GT", kTProfile, {nChAxis}); + histos.add("MCGen/Prof_Cent_MeanpT", ";cent;#LT p_{T}#GT", kTProfile, {centAxis1Per}); + histos.add("MCGen/Prof_Mult_MeanpT", ";mult;#LT p_{T}#GT", kTProfile, {nChAxis}); + + histos.add("pmeanTru_nch_etabin", ";mult;eta", kTProfile2D, {{nChAxis}, {etaBinAxis}}); + histos.add("pmeanReco_nch_etabin", ";mult;eta", kTProfile2D, {{nChAxis}, {etaBinAxis}}); + histos.add("pmeanRecoEffcorr_nch_etabin", ";mult;eta", kTProfile2D, {{nChAxis}, {etaBinAxis}}); + histos.add("pmeanMultTru_nch_etabin", ";mult;eta", kTProfile2D, {{nChAxis}, {etaBinAxis}}); + histos.add("pmeanMultReco_nch_etabin", ";mult;eta", kTProfile2D, {{nChAxis}, {etaBinAxis}}); + histos.add("pmeanMultRecoEffcorr_nch_etabin", ";mult;eta", kTProfile2D, {{nChAxis}, {etaBinAxis}}); + + histos.add("MCGen/hEtaPhiReco", ";vz;sign;pt;eta;phi", kTHnSparseF, {{vzAxis}, {chgAxis}, {pTAxis}, {etaFlatAxis}, {phiAxis}}); + histos.add("MCGen/hEtaPhiRecoEffWtd", ";vz;sign;pt;eta;phi", kTHnSparseF, {{vzAxis}, {chgAxis}, {pTAxis}, {etaFlatAxis}, {phiAxis}}); + histos.add("MCGen/hEtaPhiRecoWtd", ";vz;sign;pt;eta;phi", kTHnSparseF, {{vzAxis}, {chgAxis}, {pTAxis}, {etaFlatAxis}, {phiAxis}}); + + histos.add("Prof2D_MeanpTSub_Tru", ";cent;etaA;etaC", kTProfile3D, {{centAxis1Per}, {etaBinAxis}, {etaBinAxis}}); + histos.add("Prof2D_MeanpTSub_Reco", ";cent;etaA;etaC", kTProfile3D, {{centAxis1Per}, {etaBinAxis}, {etaBinAxis}}); + histos.add("Prof2D_MeanpTSub_RecoEffCorr", ";cent;etaA;etaC", kTProfile3D, {{centAxis1Per}, {etaBinAxis}, {etaBinAxis}}); } void declareMCFlucHists() { - histos.add("MCGen/Prof_Cent_NEta_Nsp_Nchrec", ";cent;eta;isp", kTProfile3D, {{centAxis1Per}, {etaBinAxis}, {spBinAxis}}); - histos.add("MCGen/Prof_Mult_NEta_Nsp_Nchrec", ";mult;eta;isp", kTProfile3D, {{nChAxis}, {etaBinAxis}, {spBinAxis}}); - histos.add("MCGen/Prof_Cent_NEta_Nsp_MeanpT", ";cent;eta;isp", kTProfile3D, {{centAxis1Per}, {etaBinAxis}, {spBinAxis}}); - histos.add("MCGen/Prof_Mult_NEta_Nsp_MeanpT", ";mult;eta;isp", kTProfile3D, {{nChAxis}, {etaBinAxis}, {spBinAxis}}); - - histos.add("MCGen/Prof_MeanpT_Cent_etabin_spbin", ";cent;eta;isp", kTProfile3D, {{centAxis1Per}, {etaBinAxis}, {spBinAxis}}); - histos.add("MCGen/Prof_C2_Cent_etabin_spbin", ";cent;eta;isp", kTProfile3D, {{centAxis1Per}, {etaBinAxis}, {spBinAxis}}); - histos.add("MCGen/Prof_C2Sub_Cent_etabin_spbin", ";cent;eta;isp", kTProfile3D, {{centAxis1Per}, {etaBinAxis}, {spBinAxis}}); - histos.add("MCGen/Prof_Cov_Cent_etabin_spbin", ";cent;eta;isp", kTProfile3D, {{centAxis1Per}, {etaBinAxis}, {spBinAxis}}); - histos.add("MCGen/Prof_CovFT0A_Cent_etabin_spbin", ";cent;eta;isp", kTProfile3D, {{centAxis1Per}, {etaBinAxis}, {spBinAxis}}); - histos.add("MCGen/Prof_CovFT0C_Cent_etabin_spbin", ";cent;eta;isp", kTProfile3D, {{centAxis1Per}, {etaBinAxis}, {spBinAxis}}); - - histos.add("MCGen/Prof_MeanpT_Mult_etabin_spbin", ";mult;eta;isp", kTProfile3D, {{nChAxis}, {etaBinAxis}, {spBinAxis}}); - histos.add("MCGen/Prof_C2_Mult_etabin_spbin", ";mult;eta;isp", kTProfile3D, {{nChAxis}, {etaBinAxis}, {spBinAxis}}); - histos.add("MCGen/Prof_C2Sub_Mult_etabin_spbin", ";mult;eta;isp", kTProfile3D, {{nChAxis}, {etaBinAxis}, {spBinAxis}}); - histos.add("MCGen/Prof_Cov_Mult_etabin_spbin", ";mult;eta;isp", kTProfile3D, {{nChAxis}, {etaBinAxis}, {spBinAxis}}); - histos.add("MCGen/Prof_CovFT0A_Mult_etabin_spbin", ";mult;eta;isp", kTProfile3D, {{nChAxis}, {etaBinAxis}, {spBinAxis}}); - histos.add("MCGen/Prof_CovFT0C_Mult_etabin_spbin", ";mult;eta;isp", kTProfile3D, {{nChAxis}, {etaBinAxis}, {spBinAxis}}); - - for (const auto& suf : pidSuffix) { - histos.add("hEtaPhiReco" + suf, ";vz;sign;pt;eta;phi", kTHnSparseF, {{vzAxis}, {chgAxis}, {pTAxis}, {etaAxis}, {phiAxis}}); - histos.add("hEtaPhiRecoEffWtd" + suf, ";vz;sign;pt;eta;phi", kTHnSparseF, {{vzAxis}, {chgAxis}, {pTAxis}, {etaAxis}, {phiAxis}}); - histos.add("hEtaPhiRecoWtd" + suf, ";vz;sign;pt;eta;phi", kTHnSparseF, {{vzAxis}, {chgAxis}, {pTAxis}, {etaAxis}, {phiAxis}}); - - histos.add(Form("MCGen/Prof_C2Sub2D_Cent_etaA_etaC%s", suf.c_str()), ";cent;etaA;etaC", kTProfile3D, {{centAxis1Per}, {etaAxis}, {etaAxis}}); - histos.add(Form("MCGen/Prof_GapSum2D%s", suf.c_str()), ";cent;gap;sum", kTProfile3D, {{centAxis1Per}, {gapAxis}, {sumAxis}}); - histos.add(Form("MCGen/Prof_Cov2D_Cent_etaA_etaC%s", suf.c_str()), ";cent;etaA;etaC", kTProfile3D, {{centAxis1Per}, {etaAxis}, {etaAxis}}); - histos.add(Form("MCGen/Prof_CovFT0A2D_Cent_etaA_etaC%s", suf.c_str()), ";cent;etaA;etaC", kTProfile3D, {{centAxis1Per}, {etaAxis}, {etaAxis}}); - histos.add(Form("MCGen/Prof_CovFT0C2D_Cent_etaA_etaC%s", suf.c_str()), ";cent;etaA;etaC", kTProfile3D, {{centAxis1Per}, {etaAxis}, {etaAxis}}); - } + histos.add("MCGen/Prof_Cent_NEta_Nchrec", ";cent;eta", kTProfile2D, {{centAxis1Per}, {etaBinAxis}}); + histos.add("MCGen/Prof_Mult_NEta_Nchrec", ";mult;eta", kTProfile2D, {{nChAxis}, {etaBinAxis}}); + histos.add("MCGen/Prof_Cent_NEta_MeanpT", ";cent;eta", kTProfile2D, {{centAxis1Per}, {etaBinAxis}}); + histos.add("MCGen/Prof_Mult_NEta_MeanpT", ";mult;eta", kTProfile2D, {{nChAxis}, {etaBinAxis}}); + + histos.add("MCGen/Prof_MeanpT_Cent_etabin", ";cent;eta", kTProfile2D, {{centAxis1Per}, {etaBinAxis}}); + histos.add("MCGen/Prof_C2_Cent_etabin", ";cent;eta", kTProfile2D, {{centAxis1Per}, {etaBinAxis}}); + histos.add("MCGen/Prof_C2Sub_Cent_etabin", ";cent;eta", kTProfile2D, {{centAxis1Per}, {etaBinAxis}}); + histos.add("MCGen/Prof_Cov_Cent_etabin", ";cent;eta", kTProfile2D, {{centAxis1Per}, {etaBinAxis}}); + histos.add("MCGen/Prof_CovFT0A_Cent_etabin", ";cent;eta", kTProfile2D, {{centAxis1Per}, {etaBinAxis}}); + histos.add("MCGen/Prof_CovFT0C_Cent_etabin", ";cent;eta", kTProfile2D, {{centAxis1Per}, {etaBinAxis}}); + + histos.add("MCGen/Prof_MeanpT_Mult_etabin", ";mult;eta", kTProfile2D, {{nChAxis}, {etaBinAxis}}); + histos.add("MCGen/Prof_C2_Mult_etabin", ";mult;eta", kTProfile2D, {{nChAxis}, {etaBinAxis}}); + histos.add("MCGen/Prof_C2Sub_Mult_etabin", ";mult;eta", kTProfile2D, {{nChAxis}, {etaBinAxis}}); + histos.add("MCGen/Prof_Cov_Mult_etabin", ";mult;eta", kTProfile2D, {{nChAxis}, {etaBinAxis}}); + histos.add("MCGen/Prof_CovFT0A_Mult_etabin", ";mult;eta", kTProfile2D, {{nChAxis}, {etaBinAxis}}); + histos.add("MCGen/Prof_CovFT0C_Mult_etabin", ";mult;eta", kTProfile2D, {{nChAxis}, {etaBinAxis}}); + + histos.add("MCGen/hEtaPhiReco", ";vz;sign;pt;eta;phi", kTHnSparseF, {{vzAxis}, {chgAxis}, {pTAxis}, {etaFlatAxis}, {phiAxis}}); + histos.add("MCGen/hEtaPhiRecoEffWtd", ";vz;sign;pt;eta;phi", kTHnSparseF, {{vzAxis}, {chgAxis}, {pTAxis}, {etaFlatAxis}, {phiAxis}}); + histos.add("MCGen/hEtaPhiRecoWtd", ";vz;sign;pt;eta;phi", kTHnSparseF, {{vzAxis}, {chgAxis}, {pTAxis}, {etaFlatAxis}, {phiAxis}}); + + histos.add("MCGen/Prof_C2Sub2D_Cent_etaA_etaC", ";cent;etaA;etaC", kTProfile3D, {{centAxis1Per}, {etaAxis}, {etaAxis}}); + histos.add("MCGen/Prof_GapSum2D", ";cent;gap;sum", kTProfile3D, {{centAxis1Per}, {gapAxis}, {sumAxis}}); + histos.add("MCGen/Prof_Cov2D_Cent_etaA_etaC", ";cent;etaA;etaC", kTProfile3D, {{centAxis1Per}, {etaAxis}, {etaAxis}}); + histos.add("MCGen/Prof_CovFT0A2D_Cent_etaA_etaC", ";cent;etaA;etaC", kTProfile3D, {{centAxis1Per}, {etaAxis}, {etaAxis}}); + histos.add("MCGen/Prof_CovFT0C2D_Cent_etaA_etaC", ";cent;etaA;etaC", kTProfile3D, {{centAxis1Per}, {etaAxis}, {etaAxis}}); } void declareDataGetFlatHists() { - for (const auto& suf : pidSuffix) { - histos.add("hEtaPhiReco" + suf, ";vz;sign;pt;eta;phi", kTHnSparseF, {{vzAxis}, {chgAxis}, {pTAxis}, {etaAxis}, {phiAxis}}); - histos.add("hEtaPhiRecoEffWtd" + suf, ";vz;sign;pt;eta;phi", kTHnSparseF, {{vzAxis}, {chgAxis}, {pTAxis}, {etaAxis}, {phiAxis}}); - histos.add("hEtaPhiRecoWtd" + suf, ";vz;sign;pt;eta;phi", kTHnSparseF, {{vzAxis}, {chgAxis}, {pTAxis}, {etaAxis}, {phiAxis}}); - } + histos.add("hEtaPhiReco", ";vz;sign;pt;eta;phi", kTHnSparseF, {{vzAxis}, {chgAxis}, {pTAxis}, {etaFlatAxis}, {phiAxis}}); + histos.add("hEtaPhiRecoEffWtd", ";vz;sign;pt;eta;phi", kTHnSparseF, {{vzAxis}, {chgAxis}, {pTAxis}, {etaFlatAxis}, {phiAxis}}); + histos.add("hEtaPhiRecoWtd", ";vz;sign;pt;eta;phi", kTHnSparseF, {{vzAxis}, {chgAxis}, {pTAxis}, {etaFlatAxis}, {phiAxis}}); histos.add("hnTrkPVZDC", ";N_{PV};ZDC_{A+C}", kTH2F, {{nChAxis2}, {200, 0, 3000}}); histos.add("hNchZDC", ";N_{trk};ZDC_{A+C}", kTH2F, {{nChAxis2}, {200, 0, 30000}}); } @@ -993,76 +766,104 @@ struct RadialFlowDecorr { { histos.add("pmeanFT0Amultpv", "N_{PV}; AmplitudeA", kTProfile, {nChAxis}); histos.add("pmeanFT0A_cent", "cent; AmplitudeA", kTProfile, {centAxis1Per}); - histos.add("pmeanFT0Cmultpv", "N_{PV}; AmplitudeA", kTProfile, {nChAxis}); - histos.add("pmeanFT0C_cent", "cent; AmplitudeA", kTProfile, {centAxis1Per}); + histos.add("pmeanFT0Cmultpv", "N_{PV}; AmplitudeC", kTProfile, {nChAxis}); + histos.add("pmeanFT0C_cent", "cent; AmplitudeC", kTProfile, {centAxis1Per}); histos.add("pmean_cent_id_eta_FT0", ";cent;channel id; #eta;amplitude", kTProfile3D, {{centAxis1Per}, {200, -0.5, 199.5}, {100, -5.0, 5.0}}); histos.add("h3_cent_id_eta_FT0", ";cent;channel id; #eta", kTH3F, {{centAxis1Per}, {200, -0.5, 199.5}, {100, -5.0, 5.0}}); - histos.add("Prof_Cent_Nsp_Nchrec", ";cent;Species;#LT N_{PV}#GT", kTProfile2D, {{centAxis1Per}, {spBinAxis}}); - histos.add("Prof_Mult_Nsp_Nchrec", ";N_{PV};Species;#LT N_{PV}#GT", kTProfile2D, {{nChAxis}, {spBinAxis}}); - histos.add("Prof_Cent_Nsp_MeanpT", ";cent;Species;#LT p_{T}#GT", kTProfile2D, {{centAxis1Per}, {spBinAxis}}); - histos.add("Prof_Mult_Nsp_MeanpT", ";N_{PV};Species;#LT p_{T}#GT", kTProfile2D, {{nChAxis}, {spBinAxis}}); + histos.add("Prof_Cent_Nchrec", ";cent;#LT N_{PV}#GT", kTProfile, {centAxis1Per}); + histos.add("Prof_Mult_Nchrec", ";N_{PV};#LT N_{PV}#GT", kTProfile, {nChAxis}); + histos.add("Prof_Cent_MeanpT", ";cent;#LT p_{T}#GT", kTProfile, {centAxis1Per}); + histos.add("Prof_Mult_MeanpT", ";N_{PV};#LT p_{T}#GT", kTProfile, {nChAxis}); - histos.add("pmean_nch_etabin_spbin", ";N_{PV};#eta-bin;Species", kTProfile3D, {{nChAxis}, {{etaBinAxis}}, {spBinAxis}}); - histos.add("pmeanMult_nch_etabin_spbin", ";N_{PV};#eta-bin;Species", kTProfile3D, {{nChAxis}, {{etaBinAxis}}, {spBinAxis}}); - histos.add("pmean_cent_etabin_spbin", ";Centrality (%) ;#eta-bin;Species", kTProfile3D, {{centAxis1Per}, {{etaBinAxis}}, {spBinAxis}}); - histos.add("pmeanMult_cent_etabin_spbin", ";Centrality (%) ;#eta-bin;Species", kTProfile3D, {{centAxis1Per}, {{etaBinAxis}}, {spBinAxis}}); + histos.add("pmean_nch_etabin", ";N_{PV};#eta-bin", kTProfile2D, {{nChAxis}, {etaBinAxis}}); + histos.add("pmeanMult_nch_etabin", ";N_{PV};#eta-bin", kTProfile2D, {{nChAxis}, {etaBinAxis}}); + histos.add("pmean_cent_etabin", ";Centrality (%);#eta-bin", kTProfile2D, {{centAxis1Per}, {etaBinAxis}}); + histos.add("pmeanMult_cent_etabin", ";Centrality (%);#eta-bin", kTProfile2D, {{centAxis1Per}, {etaBinAxis}}); - for (const auto& suf : pidSuffix) { - histos.add("hEtaPhiReco" + suf, ";vz;sign;pt;eta;phi", kTHnSparseF, {{vzAxis}, {chgAxis}, {pTAxis}, {etaAxis}, {phiAxis}}); - histos.add("hEtaPhiRecoEffWtd" + suf, ";vz;sign;pt;eta;phi", kTHnSparseF, {{vzAxis}, {chgAxis}, {pTAxis}, {etaAxis}, {phiAxis}}); - histos.add("hEtaPhiRecoWtd" + suf, ";vz;sign;pt;eta;phi", kTHnSparseF, {{vzAxis}, {chgAxis}, {pTAxis}, {etaAxis}, {phiAxis}}); + histos.add("hEtaPhiReco", ";vz;sign;pt;eta;phi", kTHnSparseF, {{vzAxis}, {chgAxis}, {pTAxis}, {etaFlatAxis}, {phiAxis}}); + histos.add("hEtaPhiRecoEffWtd", ";vz;sign;pt;eta;phi", kTHnSparseF, {{vzAxis}, {chgAxis}, {pTAxis}, {etaFlatAxis}, {phiAxis}}); + histos.add("hEtaPhiRecoWtd", ";vz;sign;pt;eta;phi", kTHnSparseF, {{vzAxis}, {chgAxis}, {pTAxis}, {etaFlatAxis}, {phiAxis}}); - histos.add("Prof2D_MeanpTSub" + suf, ";cent;#eta_{A} bin;#eta_{C} bin", kTProfile3D, {{centAxis1Per}, {{etaBinAxis}}, {{etaBinAxis}}}); - } + histos.add("Prof2D_MeanpTSub", ";cent;#eta_{A} bin;#eta_{C} bin", kTProfile3D, {{centAxis1Per}, {etaBinAxis}, {etaBinAxis}}); + + histos.add("pEffWeight_pt_eta_cent", ";p_{T} (GeV/c);#eta;cent;#LT eff#GT", kTProfile3D, {{KNbinsPtRes, KPtMin, KPtMax}, {KNbinsEtaFine, -KEtaFineMax, KEtaFineMax}, {centAxis1Per}}); + histos.add("pFakeWeight_pt_eta_cent", ";p_{T} (GeV/c);#eta;cent;#LT fake#GT", kTProfile3D, {{KNbinsPtRes, KPtMin, KPtMax}, {KNbinsEtaFine, -KEtaFineMax, KEtaFineMax}, {centAxis1Per}}); + histos.add("pFlatWeight_pt_eta_cent", ";p_{T} (GeV/c);#eta;cent;#LT w_{#phi}#GT", kTProfile3D, {{KNbinsPtRes, KPtMin, KPtMax}, {KNbinsEtaFine, -KEtaFineMax, KEtaFineMax}, {centAxis1Per}}); } void declareDataFlucHists() { - histos.add("Prof_MeanpT_Cent_etabin_spbin", ";cent;#eta-bin;Species", kTProfile3D, {{centAxis1Per}, {{etaBinAxis}}, {spBinAxis}}); - histos.add("Prof_MeanpT_Mult_etabin_spbin", ";N_{PV};#eta-bin;Species", kTProfile3D, {{nChAxis}, {{etaBinAxis}}, {spBinAxis}}); - histos.add("Prof_C2_Cent_etabin_spbin", ";cent;#eta-bin;Species", kTProfile3D, {{centAxis1Per}, {{etaBinAxis}}, {spBinAxis}}); - histos.add("Prof_C2_Mult_etabin_spbin", ";N_{PV};#eta-bin;Species", kTProfile3D, {{nChAxis}, {{etaBinAxis}}, {spBinAxis}}); - - histos.add("Prof_C2Sub_Cent_etabin_spbin", ";Centrality;#eta-bin;Species", kTProfile3D, {{centAxis1Per}, {{etaBinAxis}}, {spBinAxis}}); - histos.add("Prof_C2Sub_Mult_etabin_spbin", ";N_{PV};#eta-bin;Species", kTProfile3D, {{nChAxis}, {{etaBinAxis}}, {spBinAxis}}); - histos.add("Prof_Cov_Cent_etabin_spbin", ";Centrality;#eta-bin;Species", kTProfile3D, {{centAxis1Per}, {{etaBinAxis}}, {spBinAxis}}); - histos.add("Prof_Cov_Mult_etabin_spbin", ";N_{PV};#eta-bin;Species", kTProfile3D, {{nChAxis}, {{etaBinAxis}}, {spBinAxis}}); - - histos.add("Prof_CovFT0A_Cent_etabin_spbin", ";Centrality;#eta-bin;Species", kTProfile3D, {{centAxis1Per}, {{etaBinAxis}}, {spBinAxis}}); - histos.add("Prof_CovFT0A_Mult_etabin_spbin", ";N_{PV};#eta-bin;Species", kTProfile3D, {{nChAxis}, {{etaBinAxis}}, {spBinAxis}}); - histos.add("Prof_CovFT0C_Cent_etabin_spbin", ";Centrality;#eta-bin;Species", kTProfile3D, {{centAxis1Per}, {{etaBinAxis}}, {spBinAxis}}); - histos.add("Prof_CovFT0C_Mult_etabin_spbin", ";N_{PV};#eta-bin;Species", kTProfile3D, {{nChAxis}, {{etaBinAxis}}, {spBinAxis}}); - - for (const auto& suf : pidSuffix) { - histos.add("hEtaPhiReco" + suf, ";vz;sign;pt;eta;phi", kTHnSparseF, {{vzAxis}, {chgAxis}, {pTAxis}, {etaAxis}, {phiAxis}}); - histos.add("hEtaPhiRecoEffWtd" + suf, ";vz;sign;pt;eta;phi", kTHnSparseF, {{vzAxis}, {chgAxis}, {pTAxis}, {etaAxis}, {phiAxis}}); - histos.add("hEtaPhiRecoWtd" + suf, ";vz;sign;pt;eta;phi", kTHnSparseF, {{vzAxis}, {chgAxis}, {pTAxis}, {etaAxis}, {phiAxis}}); - histos.add("Prof_C2Sub2D_Cent_etaA_etaC" + suf, ";cent;#eta_{A};#eta_{C}", kTProfile3D, {{centAxis1Per}, {etaAxis}, {etaAxis}}); - histos.add("Prof_GapSum2D" + suf, ";cent;#Delta#eta (Gap);#Sigma#eta/2 (Sum)", kTProfile3D, {{centAxis1Per}, {gapAxis}, {sumAxis}}); - histos.add("Prof_Cov2D_Cent_etaA_etaC" + suf, ";cent;#eta_{A} bin;#eta_{C} bin", kTProfile3D, {{centAxis1Per}, {etaAxis}, {etaAxis}}); - histos.add("Prof_CovFT0A2D_Cent_etaA_etaC" + suf, ";cent;#eta_{A};#eta_{B}", kTProfile3D, {{centAxis1Per}, {etaAxis}, {etaAxis}}); - histos.add("Prof_CovFT0C2D_Cent_etaA_etaC" + suf, ";cent;#eta_{A};#eta_{B}", kTProfile3D, {{centAxis1Per}, {etaAxis}, {etaAxis}}); + histos.add("Prof_MeanpT_Cent_etabin", ";cent;#eta-bin", kTProfile2D, {{centAxis1Per}, {etaBinAxis}}); + histos.add("Prof_MeanpT_Mult_etabin", ";N_{PV};#eta-bin", kTProfile2D, {{nChAxis}, {etaBinAxis}}); + histos.add("Prof_C2_Cent_etabin", ";cent;#eta-bin", kTProfile2D, {{centAxis1Per}, {etaBinAxis}}); + histos.add("Prof_C2_Mult_etabin", ";N_{PV};#eta-bin", kTProfile2D, {{nChAxis}, {etaBinAxis}}); + histos.add("Prof_C2Sub_Cent_etabin", ";Centrality;#eta-bin", kTProfile2D, {{centAxis1Per}, {etaBinAxis}}); + histos.add("Prof_C2Sub_Mult_etabin", ";N_{PV};#eta-bin", kTProfile2D, {{nChAxis}, {etaBinAxis}}); + histos.add("Prof_Cov_Cent_etabin", ";Centrality;#eta-bin", kTProfile2D, {{centAxis1Per}, {etaBinAxis}}); + histos.add("Prof_Cov_Mult_etabin", ";N_{PV};#eta-bin", kTProfile2D, {{nChAxis}, {etaBinAxis}}); + histos.add("Prof_CovFT0A_Cent_etabin", ";Centrality;#eta-bin", kTProfile2D, {{centAxis1Per}, {etaBinAxis}}); + histos.add("Prof_CovFT0A_Mult_etabin", ";N_{PV};#eta-bin", kTProfile2D, {{nChAxis}, {etaBinAxis}}); + histos.add("Prof_CovFT0C_Cent_etabin", ";Centrality;#eta-bin", kTProfile2D, {{centAxis1Per}, {etaBinAxis}}); + histos.add("Prof_CovFT0C_Mult_etabin", ";N_{PV};#eta-bin", kTProfile2D, {{nChAxis}, {etaBinAxis}}); + + histos.add("hEtaPhiReco", ";vz;sign;pt;eta;phi", kTHnSparseF, {{vzAxis}, {chgAxis}, {pTAxis}, {etaFlatAxis}, {phiAxis}}); + histos.add("hEtaPhiRecoEffWtd", ";vz;sign;pt;eta;phi", kTHnSparseF, {{vzAxis}, {chgAxis}, {pTAxis}, {etaFlatAxis}, {phiAxis}}); + histos.add("hEtaPhiRecoWtd", ";vz;sign;pt;eta;phi", kTHnSparseF, {{vzAxis}, {chgAxis}, {pTAxis}, {etaFlatAxis}, {phiAxis}}); + + histos.add("Prof_C2Sub2D_Cent_etaA_etaC", ";cent;#eta_{A};#eta_{C}", kTProfile3D, {{centAxis1Per}, {etaAxis}, {etaAxis}}); + histos.add("Prof_GapSum2D", ";cent;#Delta#eta (Gap);#Sigma#eta/2 (Sum)", kTProfile3D, {{centAxis1Per}, {gapAxis}, {sumAxis}}); + histos.add("Prof_Cov2D_Cent_etaA_etaC", ";cent;#eta_{A};#eta_{C}", kTProfile3D, {{centAxis1Per}, {etaAxis}, {etaAxis}}); + histos.add("Prof_CovFT0A2D_Cent_etaA_etaC", ";cent;#eta_{A};#eta_{C}", kTProfile3D, {{centAxis1Per}, {etaAxis}, {etaAxis}}); + histos.add("Prof_CovFT0C2D_Cent_etaA_etaC", ";cent;#eta_{A};#eta_{C}", kTProfile3D, {{centAxis1Per}, {etaAxis}, {etaAxis}}); + } + + // 30 Poisson-bootstrap replicas of every final fluctuation observable + // (base data run only). Filled by pointer to bypass the compile-time HIST() + // macro, which cannot take a runtime sample index. + void declareBootstrapHists() + { + for (int s = 0; s < nBoot; ++s) { + bs.meanpTCent[s] = histos.add(Form("Bootstrap/Prof_MeanpT_Cent_etabin_sample%d", s), ";cent;#eta-bin", kTProfile2D, {{centAxis1Per}, {etaBinAxis}}); + bs.meanpTMult[s] = histos.add(Form("Bootstrap/Prof_MeanpT_Mult_etabin_sample%d", s), ";N_{PV};#eta-bin", kTProfile2D, {{nChAxis}, {etaBinAxis}}); + bs.c2Cent[s] = histos.add(Form("Bootstrap/Prof_C2_Cent_etabin_sample%d", s), ";cent;#eta-bin", kTProfile2D, {{centAxis1Per}, {etaBinAxis}}); + bs.c2Mult[s] = histos.add(Form("Bootstrap/Prof_C2_Mult_etabin_sample%d", s), ";N_{PV};#eta-bin", kTProfile2D, {{nChAxis}, {etaBinAxis}}); + bs.c2SubCent[s] = histos.add(Form("Bootstrap/Prof_C2Sub_Cent_etabin_sample%d", s), ";cent;#eta-bin", kTProfile2D, {{centAxis1Per}, {etaBinAxis}}); + bs.c2SubMult[s] = histos.add(Form("Bootstrap/Prof_C2Sub_Mult_etabin_sample%d", s), ";N_{PV};#eta-bin", kTProfile2D, {{nChAxis}, {etaBinAxis}}); + bs.covCent[s] = histos.add(Form("Bootstrap/Prof_Cov_Cent_etabin_sample%d", s), ";cent;#eta-bin", kTProfile2D, {{centAxis1Per}, {etaBinAxis}}); + bs.covMult[s] = histos.add(Form("Bootstrap/Prof_Cov_Mult_etabin_sample%d", s), ";N_{PV};#eta-bin", kTProfile2D, {{nChAxis}, {etaBinAxis}}); + bs.covFT0ACent[s] = histos.add(Form("Bootstrap/Prof_CovFT0A_Cent_etabin_sample%d", s), ";cent;#eta-bin", kTProfile2D, {{centAxis1Per}, {etaBinAxis}}); + bs.covFT0AMult[s] = histos.add(Form("Bootstrap/Prof_CovFT0A_Mult_etabin_sample%d", s), ";N_{PV};#eta-bin", kTProfile2D, {{nChAxis}, {etaBinAxis}}); + bs.covFT0CCent[s] = histos.add(Form("Bootstrap/Prof_CovFT0C_Cent_etabin_sample%d", s), ";cent;#eta-bin", kTProfile2D, {{centAxis1Per}, {etaBinAxis}}); + bs.covFT0CMult[s] = histos.add(Form("Bootstrap/Prof_CovFT0C_Mult_etabin_sample%d", s), ";N_{PV};#eta-bin", kTProfile2D, {{nChAxis}, {etaBinAxis}}); + bs.c2Sub2D[s] = histos.add(Form("Bootstrap/Prof_C2Sub2D_Cent_etaA_etaC_sample%d", s), ";cent;#eta_{A};#eta_{C}", kTProfile3D, {{centAxis1Per}, {etaAxis}, {etaAxis}}); + bs.gapSum2D[s] = histos.add(Form("Bootstrap/Prof_GapSum2D_sample%d", s), ";cent;gap;sum", kTProfile3D, {{centAxis1Per}, {gapAxis}, {sumAxis}}); + bs.cov2D[s] = histos.add(Form("Bootstrap/Prof_Cov2D_Cent_etaA_etaC_sample%d", s), ";cent;#eta_{A};#eta_{C}", kTProfile3D, {{centAxis1Per}, {etaAxis}, {etaAxis}}); + bs.covFT0A2D[s] = histos.add(Form("Bootstrap/Prof_CovFT0A2D_Cent_etaA_etaC_sample%d", s), ";cent;#eta_{A};#eta_{C}", kTProfile3D, {{centAxis1Per}, {etaAxis}, {etaAxis}}); + bs.covFT0C2D[s] = histos.add(Form("Bootstrap/Prof_CovFT0C2D_Cent_etaA_etaC_sample%d", s), ";cent;#eta_{A};#eta_{C}", kTProfile3D, {{centAxis1Per}, {etaAxis}, {etaAxis}}); } } + // =========================================================================== + // CCDB helpers + // =========================================================================== THnSparseF* buildWeightMapFromRaw(THnSparseF* hRaw, const char* mapName) { if (!hRaw) { LOGF(error, "Raw eta-phi map for '%s' is null; no flattening will be applied.", mapName); return nullptr; } - auto hWMap = reinterpret_cast(hRaw->Clone(mapName)); + auto hWMap = dynamic_cast(hRaw->Clone(mapName)); hWMap->SetTitle(Form("Flattening Weight Map %s (w_{#phi} = / N_{#phi})", mapName)); hWMap->Reset(); - auto axV = hRaw->GetAxis(0); // Vz - auto axChg = hRaw->GetAxis(1); // Charge - auto axPt = hRaw->GetAxis(2); // Pt - auto axE = hRaw->GetAxis(3); // Eta - auto axP = hRaw->GetAxis(4); // Phi + auto axV = hRaw->GetAxis(0); + auto axChg = hRaw->GetAxis(1); + auto axPt = hRaw->GetAxis(2); + auto axE = hRaw->GetAxis(3); + auto axP = hRaw->GetAxis(4); - int bins[5]; + std::array bins{}; for (int iv = 1; iv <= axV->GetNbins(); ++iv) { bins[0] = iv; for (int ichg = 1; ichg <= axChg->GetNbins(); ++ichg) { @@ -1075,47 +876,50 @@ struct RadialFlowDecorr { int nphi = axP->GetNbins(); for (int ip = 1; ip <= nphi; ++ip) { bins[4] = ip; - sum += hRaw->GetBinContent(bins); + sum += hRaw->GetBinContent(bins.data()); } const double avg = (nphi > 0 ? sum / nphi : 0.0); for (int ip = 1; ip <= nphi; ++ip) { bins[4] = ip; - const double raw = hRaw->GetBinContent(bins); + const double raw = hRaw->GetBinContent(bins.data()); const double w = (avg > 0.0 && raw > 0.0) ? (avg / raw) : 1.0; - hWMap->SetBinContent(bins, w); + hWMap->SetBinContent(bins.data(), w); } } } } } - LOGF(info, "Flattening weight map '%s' built.", mapName); return hWMap; } - inline void loadTProfile3D(TDirectory* dir, const char* name, TProfile3D*& target) + template + void loadProfileFromList(TList* src, const char* name, TP*& target) { - if (!dir) { - LOGF(error, "loadTProfile3D: directory is null for object %s", name); + if (!src) { return; } - auto* obj = dir->Get(name); + auto* obj = src->FindObject(name); if (!obj) { - LOGF(error, "loadTProfile3D: object '%s' not found in directory %s", name, dir->GetName()); + LOGF(error, "Profile %s missing in CCDB TList", name); return; } - auto* prof = dynamic_cast(obj); - if (!prof) { - LOGF(error, "loadTProfile3D: object '%s' is not a TProfile3D (it is %s)", name, obj->ClassName()); + auto* tp = dynamic_cast(obj); + if (!tp) { + LOGF(error, "%s is not the expected profile type (it is %s)", name, obj->ClassName()); return; } - target = reinterpret_cast(prof->Clone(Form("%s_clone", name))); + target = dynamic_cast(tp->Clone()); target->SetDirectory(nullptr); - LOGF(info, "Loaded TProfile3D '%s' with entries = %.0f", name, target->GetEntries()); + LOGF(info, "Loaded %s from list", name); } + // =========================================================================== + // init + // =========================================================================== void init(InitContext&) { + // Nch axes by system if (cfgSys == kPbPb) { nChAxis = {cfgNchPbMax / 2, KBinOffset, cfgNchPbMax + KBinOffset, "Nch", "PV-contributor track multiplicity"}; nChAxis2 = {cfgNchPbMax / 4, KBinOffset, cfgNchPbMax + KBinOffset, "Nch", "PV-contributor track multiplicity"}; @@ -1124,6 +928,44 @@ struct RadialFlowDecorr { nChAxis2 = {cfgNchOMax, KBinOffset, cfgNchOMax + KBinOffset, "Nch", "PV-contributor track multiplicity"}; } + // ---- observable eta binning: 0.1 for the base data run, 0.2 otherwise ---- + const bool isDataRun = (cfgRunGetDataFlat || cfgRunDataMean || cfgRunDataFluc); + const bool useFineBinning = isDataRun && (cfgSystType == kSystBase); + { + const float lo = -cfgCutEta; + const float hi = cfgCutEta; + const float width = useFineBinning ? 0.1f : 0.2f; + const int nbins = static_cast(std::lround((hi - lo) / width)); + etaLw.clear(); + etaUp.clear(); + etaLw.push_back(lo); // index 0: full-range reference bin + etaUp.push_back(hi); + for (int i = 0; i < nbins; ++i) { + etaLw.push_back(lo + i * width); + etaUp.push_back(lo + (i + 1) * width); + } + nEta = nbins + 1; + if (nEta > KNEtaMax) { + LOGF(fatal, "nEta=%d exceeds KNEtaMax=%d", nEta, KNEtaMax); + } + + std::vector obsEdges; + obsEdges.reserve(nEta); + obsEdges.push_back(etaLw[1]); + for (int i = 1; i < nEta; ++i) { + obsEdges.push_back(etaUp[i]); + } + etaAxis = AxisSpec{obsEdges, "#eta"}; + etaBinAxis = AxisSpec{nEta + 1, -0.5, static_cast(nEta) + 0.5, "#eta bin Number"}; + LOGF(info, "Observable eta binning: %d bins of width %.2f (+ reference), nEta=%d", nbins, width, nEta); + } + + // bootstrap active only for the base data fluctuation pass + doBoot = cfgRunDataFluc && (cfgSystType == kSystBase) && (cfgNBootstrap > 0); + nBoot = std::min(cfgNBootstrap, KMaxBoot); + rng.SetSeed(cfgBootstrapSeed); + + // ---- CCDB ---- ccdb->setURL(cfgCCDBurl.value); ccdb->setCaching(true); ccdb->setLocalObjectValidityChecking(); @@ -1133,7 +975,7 @@ struct RadialFlowDecorr { loadAlignParam(now); ft0Det.calculateChannelCenter(); - std::string sysDir = ""; + std::string sysDir; switch (cfgSys) { case kPbPb: sysDir = "PbPbTest"; @@ -1150,20 +992,14 @@ struct RadialFlowDecorr { default: LOGF(fatal, "Invalid cfgSys value: %d", cfgSys.value); } - std::string pathNsig = cfgCCDBUserPath.value + "/" + sysDir + "/Job0_nSigMaps"; std::string pathEff = cfgCCDBUserPath.value + "/" + sysDir + "/Job1_EffMaps"; std::string pathMCFlat = cfgCCDBUserPath.value + "/" + sysDir + "/Job1_MCFlatMaps"; std::string pathMCMean = cfgCCDBUserPath.value + "/" + sysDir + "/Job2_MCMean"; - - std::string pathDataNsig = cfgCCDBUserPath.value + "/" + sysDir + "/Job0_DatanSigMaps"; std::string pathDataFlat = cfgCCDBUserPath.value + "/" + sysDir + "/Job1_DataFlatMaps"; std::string pathDataMean = cfgCCDBUserPath.value + "/" + sysDir + "/Job2_DataMean"; + // ---- declarations ---- declareCommonQA(); - if (cfgRunMCGetNSig || cfgRunGetEff || cfgRunDataGetNSig || cfgRunGetDataFlat) { - declarenSigHists(); - } - if (cfgRunMCMean || cfgRunMCFluc || cfgRunGetEff) { declareMCCommonHists(); } @@ -1182,9 +1018,6 @@ struct RadialFlowDecorr { histos.addClone("MCGen/", "MCReco/"); histos.addClone("MCGen/", "MCRecoEffCorr/"); } - if (cfgRunDataGetNSig) { - declareDataGetFlatHists(); - } if (cfgRunGetDataFlat) { declareDataGetFlatHists(); } @@ -1193,269 +1026,160 @@ struct RadialFlowDecorr { } if (cfgRunDataFluc) { declareDataFlucHists(); + if (doBoot) { + declareBootstrapHists(); + } } - if (!cfgRunGetEff && !cfgRunMCGetNSig && (cfgEff)) { - TList* lst = ccdb->getForTimeStamp(pathEff, now); - + // ---- efficiency/fake maps (inclusive) ---- + if (!cfgRunGetEff && cfgEff) { + auto* lst = ccdb->getForTimeStamp(pathEff, now); if (!lst) { LOGF(fatal, "Efficiency maps required but CCDB list is null at %s!", pathEff.c_str()); + return; } - LOGF(info, "Loading Eff/Fake maps from TList for all species..."); - - auto loadEffFakeForPID = [&](PIDIdx pidType) { - std::string suffix = pidSuffix[pidType]; - std::string hEffNumName = "h3_RecoMatchedToPrimary" + suffix; - std::string hEffDenName = "h3_AllPrimary" + suffix; - std::string hFakeNumSecName = "h3_RecoUnMatchedToPrimary_Secondary" + suffix; - std::string hFakeNumFakName = "h3_RecoUnMatchedToPrimary_Fake" + suffix; - std::string hFakeNumFakName2 = "h3_RecoMatchedToPrimary_MisID" + suffix; - std::string hFakeDenName = "h3_AllReco" + suffix; + auto* hNum = dynamic_cast(lst->FindObject("h3_RecoMatchedToPrimary")); + auto* hDen = dynamic_cast(lst->FindObject("h3_AllPrimary")); + if (hNum && hDen) { + state.hEff = dynamic_cast(hNum->Clone("hEff")); + state.hEff->SetDirectory(nullptr); + state.hEff->Divide(hDen); + } else { + LOGF(error, "Missing CCDB objects for efficiency (h3_RecoMatchedToPrimary / h3_AllPrimary)."); + } - auto* hNum = reinterpret_cast(lst->FindObject(hEffNumName.c_str())); - auto* hDen = reinterpret_cast(lst->FindObject(hEffDenName.c_str())); + auto* hNumS = dynamic_cast(lst->FindObject("h3_RecoUnMatchedToPrimary_Secondary")); + auto* hNumF = dynamic_cast(lst->FindObject("h3_RecoUnMatchedToPrimary_Fake")); + auto* hDenF = dynamic_cast(lst->FindObject("h3_AllReco")); + if (hNumS && hNumF && hDenF) { + state.hFake = dynamic_cast(hNumS->Clone("hFake")); + state.hFake->Add(hNumF); + state.hFake->SetDirectory(nullptr); + state.hFake->Divide(hDenF); + } else { + LOGF(error, "Missing CCDB objects for fakes."); + } + } - if (hNum && hDen) { - state.hEff[pidType] = reinterpret_cast(hNum->Clone(Form("hEff%s", suffix.c_str()))); - state.hEff[pidType]->SetDirectory(nullptr); - state.hEff[pidType]->Divide(hDen); + // ---- flattening maps (inclusive) ---- + if (!cfgRunGetEff && cfgFlat) { + if (cfgRunDataMean || cfgRunDataFluc) { + auto* lstDataFlat = ccdb->getForTimeStamp(pathDataFlat, now); + if (lstDataFlat) { + std::string hName = cfgEff ? "hEtaPhiRecoWtd" : "hEtaPhiReco"; + auto* hRaw = dynamic_cast(lstDataFlat->FindObject(hName.c_str())); + if (hRaw) { + state.hFlatWeight = buildWeightMapFromRaw(hRaw, "hFlatWeight"); + } else { + LOGF(error, "Data flattening map '%s' not found.", hName.c_str()); + } } else { - LOGF(error, "Missing CCDB objects for efficiency. Checked: %s, %s", hEffNumName.c_str(), hEffDenName.c_str()); + LOGF(error, "Could not retrieve Data Flattening TList from: %s", pathDataFlat.c_str()); } - - auto* hNumS = reinterpret_cast(lst->FindObject(hFakeNumSecName.c_str())); - auto* hNumF = reinterpret_cast(lst->FindObject(hFakeNumFakName.c_str())); - auto* hNumF2 = reinterpret_cast(lst->FindObject(hFakeNumFakName2.c_str())); - auto* hDenF = reinterpret_cast(lst->FindObject(hFakeDenName.c_str())); - - if (hNumS && hNumF && hDenF) { - state.hFake[pidType] = reinterpret_cast(hNumS->Clone(Form("hFake%s", suffix.c_str()))); - state.hFake[pidType]->Add(hNumF); - if (pidType != kInclusiveIdx && hNumF2) { - state.hFake[pidType]->Add(hNumF2); + } + if (cfgRunMCMean || cfgRunMCFluc) { + auto* lstMCFlat = ccdb->getForTimeStamp(pathMCFlat, now); + if (lstMCFlat) { + std::string hName = cfgEff ? "MCReco/hEtaPhiRecoWtd" : "MCReco/hEtaPhiReco"; + auto* hRaw = dynamic_cast(lstMCFlat->FindObject(hName.c_str())); + if (hRaw) { + state.hFlatWeight = buildWeightMapFromRaw(hRaw, "hFlatWeight"); + } else { + LOGF(warning, "MC flattening source '%s' not found.", hName.c_str()); } - state.hFake[pidType]->SetDirectory(nullptr); - state.hFake[pidType]->Divide(hDenF); } else { - LOGF(error, "Missing CCDB object(s) for fakes for %s in list.", suffix.c_str()); + LOGF(error, "Could not retrieve MC Flattening TList from: %s", pathMCFlat.c_str()); } - }; - - for (int i = 0; i < KNsp; ++i) { - loadEffFakeForPID(static_cast(i)); } } - bool requiresMCMap = (cfgRunGetEff || cfgRunGetMCFlat || cfgRunMCMean || cfgRunMCFluc); - bool requiresDataMap = (cfgRunGetDataFlat || cfgRunDataMean || cfgRunDataFluc); - - if (requiresMCMap || requiresDataMap) { - std::string currentPath = requiresMCMap ? pathNsig : pathDataNsig; - TList* pidList = ccdb->getForTimeStamp(currentPath, now); - - if (!pidList) { - LOGF(warn, "nSigma maps required but CCDB list is null at %s! Using raw values.", currentPath.c_str()); - } else { - if (!pidMeanSigmaMap) { - pidMeanSigmaMap = new PIDMeanSigmaMap(); - } - LOGF(info, "Performing 2D Gaussian fits on PID maps from CCDB..."); - auto loadPIDMeans = [&](PIDIdx pidType) { - std::string suffix = pidSuffix[pidType]; - std::string hName = "h3DnsigmaTpcVsTofBefCut_Cent" + suffix; - auto* h3 = reinterpret_cast(pidList->FindObject(hName.c_str())); - if (!h3) { - LOGF(warn, " [!] PID Hist %s not found in CCDB list.", hName.c_str()); - return; - } - int nCentBins = std::min(h3->GetXaxis()->GetNbins(), PIDMeanSigmaMap::MaxCentBins - 1); - LOGF(info, " -> Species: %s (Bins: %d)", hName.c_str(), nCentBins); - for (int iCent = 1; iCent <= nCentBins; ++iCent) { - h3->GetXaxis()->SetRange(iCent, iCent); - // Projecting: Z(TPC) vs Y(TOF). Result: X_axis=TOF, Y_axis=TPC - std::unique_ptr h2(reinterpret_cast(h3->Project3D("zy"))); - if (h2) { - int binX, binY, binZ; - h2->GetMaximumBin(binX, binY, binZ); - double guessMeanTOF = h2->GetXaxis()->GetBinCenter(binX); - double guessMeanTPC = h2->GetYaxis()->GetBinCenter(binY); - TF2 f2("f2", "[0]*TMath::Gaus(x,[1],[2])*TMath::Gaus(y,[3],[4])", -3, 3, -3, 3); - f2.SetParameters(h2->GetMaximum(), guessMeanTOF, 1.0, guessMeanTPC, 1.0); - h2->Fit(&f2, "QRN"); // Q=Quiet, R=Range, N=NoDraw - pidMeanSigmaMap->meanTOF[pidType][iCent - 1] = f2.GetParameter(1); - pidMeanSigmaMap->meanTPC[pidType][iCent - 1] = f2.GetParameter(3); - pidMeanSigmaMap->sigmaTOF[pidType][iCent - 1] = std::abs(f2.GetParameter(2)); - pidMeanSigmaMap->sigmaTPC[pidType][iCent - 1] = std::abs(f2.GetParameter(4)); - if (iCent % KConstTen == 0) { - LOGF(info, " Derived: For Species: %s (Bins: %d), Mean TOF = %.3f, Mean TPC = %.3f, Sigma TOF = %.3f, Sigma TPC = %.3f", hName.c_str(), iCent - 1, f2.GetParameter(1), f2.GetParameter(3), f2.GetParameter(2), f2.GetParameter(4)); - } - } - } - }; - for (int i = 1; i < KNsp; ++i) { - loadPIDMeans(static_cast(i)); - } - + // ---- sigma-pileup limits ---- + // These were produced by the (now removed) nSigma pass. They are re-sourced + // from the flattening-map file, which also stores Hist2D_globalTracks_cent. + // Loaded only when the sigma-pileup cut is actually requested. + if (cfgApplySigPupCut) { + const bool mcSide = (cfgRunGetEff || cfgRunGetMCFlat || cfgRunMCMean || cfgRunMCFluc); + std::string limPath = mcSide ? pathMCFlat : pathDataFlat; + auto* limList = ccdb->getForTimeStamp(limPath, now); + if (limList) { auto loadLimits = [&](const char* name, std::vector>& limits, float& xMin, float& xMax) { - auto* h2 = reinterpret_cast(pidList->FindObject(name)); - if (!h2) + auto* h2 = dynamic_cast(limList->FindObject(name)); + if (!h2) { return; - + } std::unique_ptr prof(h2->ProfileX("ptmp", 1, -1, "S")); - int nBins = prof->GetNbinsX(); xMin = prof->GetXaxis()->GetXmin(); xMax = prof->GetXaxis()->GetXmax(); - limits.assign(nBins + 2, {-99999.f, 999999.f}); - for (int i = 1; i <= nBins; ++i) { float mean = prof->GetBinContent(i); float rms = prof->GetBinError(i); - limits[i] = {mean - cfgPupnSig * rms, mean + cfgPupnSig * rms}; } }; loadLimits("Hist2D_globalTracks_cent", state.mLimitsNchCent, state.mMinXNchCent, state.mMaxXNchCent); - } - } - if (!cfgRunGetEff && (cfgFlat)) { - if (cfgRunDataMean || cfgRunDataFluc) { - LOGF(info, "Data Run: Loading flattening maps from %s", pathDataFlat.c_str()); - TList* lstDataFlat = ccdb->getForTimeStamp(pathDataFlat, now); - - if (lstDataFlat) { - for (int i = 0; i < KNsp; ++i) { - std::string suffix = pidSuffix[i]; - std::string hName; - - if (cfgEff && cfgFlat) { - hName = "hEtaPhiRecoWtd" + suffix; - } else if (cfgEff) { - hName = "hEtaPhiRecoEffWtd" + suffix; - } else { - hName = "hEtaPhiReco" + suffix; - } - auto* hRaw = reinterpret_cast(lstDataFlat->FindObject(hName.c_str())); - - if (hRaw) { - state.hFlatWeight[i] = buildWeightMapFromRaw(hRaw, Form("hFlatWeight%s", suffix.c_str())); - } else { - LOGF(error, "Data flattening map '%s' not found.", hName.c_str()); - } - } - } else { - LOGF(error, "Could not retrieve Data Flattening TList from: %s", pathDataFlat.c_str()); - } - } - - if (cfgRunMCMean || cfgRunMCFluc) { - LOGF(info, "MC Run: Loading flattening maps from %s", pathMCFlat.c_str()); - TList* lstMCFlat = ccdb->getForTimeStamp(pathMCFlat, now); - - if (lstMCFlat) { - auto loadFlatForPID = [&](PIDIdx pidType) { - std::string suffix = pidSuffix[pidType]; - std::string hFlatSrcName; - if (cfgEff && cfgFlat) { - hFlatSrcName = "MCReco/hEtaPhiRecoWtd" + suffix; - } else if (cfgEff) { - hFlatSrcName = "MCReco/hEtaPhiRecoEffWtd" + suffix; - } else { - hFlatSrcName = "MCReco/hEtaPhiReco" + suffix; - } - - auto* hRaw = reinterpret_cast(lstMCFlat->FindObject(hFlatSrcName.c_str())); - - if (hRaw) { - state.hFlatWeight[pidType] = buildWeightMapFromRaw(hRaw, Form("hFlatWeight%s", suffix.c_str())); - } else { - LOGF(warning, "MC flattening source '%s' not found in list.", hFlatSrcName.c_str()); - } - }; - - for (int i = 0; i < KNsp; ++i) { - loadFlatForPID(static_cast(i)); - } - } else { - LOGF(error, "Could not retrieve MC Flattening TList from: %s", pathMCFlat.c_str()); - } - } - } - - auto loadTProfile3DFromList = [&](TList* sourceList, const char* objName, TProfile3D*& target) { - if (!sourceList) - return; - auto* tp = reinterpret_cast(sourceList->FindObject(objName)); - if (tp) { - target = reinterpret_cast(tp->Clone()); - target->SetDirectory(nullptr); - LOGF(info, "Loaded %s from list", objName); - } else { - LOGF(error, "Histogram %s missing in CCDB TList", objName); - } - }; - - auto loadTProfileFromList = [&](TList* sourceList, const char* objName, TProfile*& target) { - if (!sourceList) - return; - auto* tp = reinterpret_cast(sourceList->FindObject(objName)); - if (tp) { - target = reinterpret_cast(tp->Clone()); - target->SetDirectory(nullptr); - LOGF(info, "Loaded %s from list", objName); } else { - LOGF(error, "Histogram %s missing in CCDB TList", objName); + LOGF(warning, "sigma-pileup limits source list missing at %s; sigma-pileup cut effectively disabled.", limPath.c_str()); } - }; + } + // ---- MC mean profiles for MC fluc ---- if (cfgRunMCFluc) { LOGF(info, "Loading MC Mean profiles from CCDB path: %s", pathMCMean.c_str()); - TList* lstMCMean = ccdb->getForTimeStamp(pathMCMean, now); - + auto* lstMCMean = ccdb->getForTimeStamp(pathMCMean, now); if (lstMCMean) { - loadTProfileFromList(lstMCMean, "pmeanFT0Amultpv", state.pmeanFT0AmultpvStep2); - loadTProfileFromList(lstMCMean, "pmeanFT0Cmultpv", state.pmeanFT0CmultpvStep2); - - loadTProfile3DFromList(lstMCMean, "pmeanTru_nch_etabin_spbin", state.pmeanTruNchEtabinSpbinStep2); - loadTProfile3DFromList(lstMCMean, "pmeanReco_nch_etabin_spbin", state.pmeanRecoNchEtabinSpbinStep2); - loadTProfile3DFromList(lstMCMean, "pmeanRecoEffcorr_nch_etabin_spbin", state.pmeanRecoEffcorrNchEtabinSpbinStep2); - - loadTProfile3DFromList(lstMCMean, "pmeanMultTru_nch_etabin_spbin", state.pmeanMultTruNchEtabinSpbinStep2); - loadTProfile3DFromList(lstMCMean, "pmeanMultReco_nch_etabin_spbin", state.pmeanMultRecoNchEtabinSpbinStep2); - loadTProfile3DFromList(lstMCMean, "pmeanMultRecoEffcorr_nch_etabin_spbin", state.pmeanMultRecoEffcorrNchEtabinSpbinStep2); + loadProfileFromList(lstMCMean, "pmeanFT0Amultpv", state.pmeanFT0AmultpvStep2); + loadProfileFromList(lstMCMean, "pmeanFT0Cmultpv", state.pmeanFT0CmultpvStep2); + loadProfileFromList(lstMCMean, "pmeanTru_nch_etabin", state.pmeanTruNchEtabinStep2); + loadProfileFromList(lstMCMean, "pmeanReco_nch_etabin", state.pmeanRecoNchEtabinStep2); + loadProfileFromList(lstMCMean, "pmeanRecoEffcorr_nch_etabin", state.pmeanRecoEffcorrNchEtabinStep2); + loadProfileFromList(lstMCMean, "pmeanMultTru_nch_etabin", state.pmeanMultTruNchEtabinStep2); + loadProfileFromList(lstMCMean, "pmeanMultReco_nch_etabin", state.pmeanMultRecoNchEtabinStep2); + loadProfileFromList(lstMCMean, "pmeanMultRecoEffcorr_nch_etabin", state.pmeanMultRecoEffcorrNchEtabinStep2); } else { LOGF(error, "Could not retrieve TList for MC Mean from: %s", pathMCMean.c_str()); } } + // ---- Data mean profiles for data fluc (one per systematic variation) ---- if (cfgRunDataFluc) { - LOGF(info, "Loading Data Mean profiles from CCDB path: %s", pathDataMean.c_str()); - TList* lstDataMean = ccdb->getForTimeStamp(pathDataMean, now); - + int st = std::clamp(cfgSystType, 0, kNSystType - 1); + std::string meanPath = pathDataMean + systSuffix[st]; + LOGF(info, "Loading Data Mean profiles for systematic '%s' from: %s", + (st == kSystBase ? "Base" : systSuffix[st].c_str()), meanPath.c_str()); + auto* lstDataMean = ccdb->getForTimeStamp(meanPath, now); if (lstDataMean) { - loadTProfileFromList(lstDataMean, "pmeanFT0Amultpv", state.pmeanFT0AmultpvStep2); - loadTProfileFromList(lstDataMean, "pmeanFT0Cmultpv", state.pmeanFT0CmultpvStep2); - - loadTProfile3DFromList(lstDataMean, "pmean_nch_etabin_spbin", state.pmeanNchEtabinSpbinStep2); - loadTProfile3DFromList(lstDataMean, "pmeanMult_nch_etabin_spbin", state.pmeanMultNchEtabinSpbinStep2); + loadProfileFromList(lstDataMean, "pmeanFT0Amultpv", state.pmeanFT0AmultpvStep2); + loadProfileFromList(lstDataMean, "pmeanFT0Cmultpv", state.pmeanFT0CmultpvStep2); + loadProfileFromList(lstDataMean, "pmean_nch_etabin", state.pmeanNchEtabinStep2); + loadProfileFromList(lstDataMean, "pmeanMult_nch_etabin", state.pmeanMultNchEtabinStep2); } else { - LOGF(error, "Could not retrieve TList for Data Mean from: %s", pathDataMean.c_str()); + LOGF(error, "Could not retrieve TList for Data Mean from: %s", meanPath.c_str()); } } LOGF(info, "CCDB initialization complete for RadialFlowDecorr."); } - void processMCGetMeanNsig(MyRun3MCCollisions::iterator const& mcCollision, FilteredTCs const& mcTracks) + // =========================================================================== + // MC: build efficiency / fake maps (inclusive) + // =========================================================================== + void processGetEffHists(MyRun3MCCollisions::iterator const& mcCollision, FilteredTCs const& mcTracks, aod::McParticles const& mcParticles) { histos.fill(HIST("hVtxZ"), mcCollision.posZ()); - if (!mcCollision.has_mcCollision() || !isEventSelected(mcCollision)) + if (!mcCollision.has_mcCollision() || !isEventSelected(mcCollision)) { return; + } float cent = getCentrality(mcCollision); - if (cent > KCentMax) + if (cent > KCentMax) { return; + } float multPV = mcCollision.multNTracksPV(); + float vz = mcCollision.posZ(); + if (!isPassAddPileup(multPV, mcTracks.size(), cent)) { + return; + } histos.fill(HIST("hVtxZ_after_sel"), mcCollision.posZ()); histos.fill(HIST("hCentrality"), cent); @@ -1464,682 +1188,248 @@ struct RadialFlowDecorr { histos.fill(HIST("Hist2D_globalTracks_cent"), cent, mcTracks.size()); histos.fill(HIST("Hist2D_PVTracks_cent"), cent, multPV); - for (const auto& track : mcTracks) { - if (track.collisionId() != mcCollision.index()) + for (const auto& particle : mcParticles) { + if (particle.mcCollisionId() != mcCollision.mcCollisionId()) { + continue; + } + if (!isParticleSelected(particle) || !particle.isPhysicalPrimary()) { continue; + } + histos.fill(HIST("h3_AllPrimary"), multPV, particle.pt(), particle.eta()); + } - if (!isTrackSelected(track)) + for (const auto& track : mcTracks) { + if (track.collisionId() != mcCollision.index()) { + continue; + } + if (!isTrackSelected(track)) { continue; - fillNSigmaBefCut(track, cent); + } + + float pt = track.pt(), eta = track.eta(); + histos.fill(HIST("hPt"), pt); + histos.fill(HIST("hEta"), eta); + histos.fill(HIST("h3_AllReco"), multPV, pt, eta); + + if (track.has_mcParticle()) { + auto mcP = track.mcParticle(); + if (mcP.isPhysicalPrimary()) { + histos.fill(HIST("ptResolution"), mcP.pt(), (pt - mcP.pt()) / mcP.pt()); + histos.fill(HIST("etaResolution"), mcP.eta(), eta - mcP.eta()); + histos.fill(HIST("etaTruthReco"), mcP.eta(), eta); + histos.fill(HIST("vzResolution"), mcP.vz(), (vz - mcP.vz()) / mcP.vz()); + histos.fill(HIST("TruthTracKVz"), mcP.vz(), vz); + histos.fill(HIST("h3_RecoMatchedToPrimary"), multPV, mcP.pt(), mcP.eta()); + } else { + histos.fill(HIST("h3_RecoUnMatchedToPrimary_Secondary"), multPV, pt, eta); + } + } else { + histos.fill(HIST("h3_RecoUnMatchedToPrimary_Fake"), multPV, pt, eta); + } } } - PROCESS_SWITCH(RadialFlowDecorr, processMCGetMeanNsig, "process MC to calculate Mean values of nSig Plots", cfgRunMCGetNSig); + PROCESS_SWITCH(RadialFlowDecorr, processGetEffHists, "process MC to calculate EffWeights", cfgRunGetEff); - void processGetEffHists(MyRun3MCCollisions::iterator const& mcCollision, FilteredTCs const& mcTracks, aod::McParticles const& mcParticles) + // =========================================================================== + // MC: build flattening maps (inclusive) + // =========================================================================== + void processMCFlat(MyRun3MCCollisions::iterator const& mcCollision, FilteredTCs const& mcTracks) { histos.fill(HIST("hVtxZ"), mcCollision.posZ()); - if (!mcCollision.has_mcCollision() || !isEventSelected(mcCollision)) + if (!mcCollision.has_mcCollision() || !isEventSelected(mcCollision)) { return; + } float cent = getCentrality(mcCollision); - if (cent > KCentMax) + if (cent > KCentMax) { return; + } float multPV = mcCollision.multNTracksPV(); float vz = mcCollision.posZ(); - if (!isPassAddPileup(multPV, mcTracks.size(), cent)) + if (!isPassAddPileup(multPV, mcTracks.size(), cent)) { return; + } + histos.fill(HIST("hVtxZ_after_sel"), mcCollision.posZ()); histos.fill(HIST("hCentrality"), cent); - histos.fill(HIST("Hist2D_globalTracks_PVTracks"), multPV, mcTracks.size()); histos.fill(HIST("Hist2D_cent_nch"), mcTracks.size(), cent); histos.fill(HIST("Hist2D_globalTracks_cent"), cent, mcTracks.size()); histos.fill(HIST("Hist2D_PVTracks_cent"), cent, multPV); - for (const auto& particle : mcParticles) { - if (particle.mcCollisionId() != mcCollision.mcCollisionId()) + for (const auto& track : mcTracks) { + if (track.collisionId() != mcCollision.index()) { continue; - - if (!isParticleSelected(particle) || !particle.isPhysicalPrimary()) - continue; - - const int pdg = particle.pdgCode(); - const int absPdg = std::abs(pdg); - float pt = particle.pt(), eta = particle.eta(); - - bool isSpecies[KNsp] = { - true, // kInclusiveIdx - pdg == -KPiPlus, // kPiMinusIdx - pdg == KPiPlus, // kPiPlusIdx - absPdg == KPiPlus, // kPiAllIdx - pdg == -KKPlus, // kKaMinusIdx - pdg == KKPlus, // kKaPlusIdx - absPdg == KKPlus, // kKaAllIdx - pdg == -KProton, // kAntiPrIdx - pdg == KProton, // kPrIdx - absPdg == KProton // kPrAllIdx - }; - - histos.fill(HIST("h3_AllPrimary"), multPV, pt, eta); - if (isSpecies[kPiMinusIdx]) - histos.fill(HIST("h3_AllPrimary_PiMinus"), multPV, pt, eta); - else if (isSpecies[kPiPlusIdx]) - histos.fill(HIST("h3_AllPrimary_PiPlus"), multPV, pt, eta); - if (isSpecies[kPiAllIdx]) - histos.fill(HIST("h3_AllPrimary_PiAll"), multPV, pt, eta); - - if (isSpecies[kKaMinusIdx]) - histos.fill(HIST("h3_AllPrimary_KaMinus"), multPV, pt, eta); - else if (isSpecies[kKaPlusIdx]) - histos.fill(HIST("h3_AllPrimary_KaPlus"), multPV, pt, eta); - if (isSpecies[kKaAllIdx]) - histos.fill(HIST("h3_AllPrimary_KaAll"), multPV, pt, eta); - - if (isSpecies[kAntiPrIdx]) - histos.fill(HIST("h3_AllPrimary_AntiPr"), multPV, pt, eta); - else if (isSpecies[kPrIdx]) - histos.fill(HIST("h3_AllPrimary_Pr"), multPV, pt, eta); - if (isSpecies[kPrAllIdx]) - histos.fill(HIST("h3_AllPrimary_PrAll"), multPV, pt, eta); - } - - for (const auto& track : mcTracks) { - if (track.collisionId() != mcCollision.index()) - continue; - - if (!isTrackSelected(track)) - continue; - - float pt = track.pt(), eta = track.eta(); - histos.fill(HIST("hPt"), pt); - histos.fill(HIST("hEta"), eta); - auto sign = track.sign(); - fillNSigmaBefCut(track, cent); - - int id = identifyTrack(track, cent); - bool isPi = (id == KPidPionOne); - bool isKa = (id == KPidKaonTwo); - bool isPr = (id == KPidProtonThree); - bool isSpecies[KNsp] = { - true, - isPi && sign < 0, isPi && sign > 0, isPi, - isKa && sign < 0, isKa && sign > 0, isKa, - isPr && sign < 0, isPr && sign > 0, isPr}; - - fillNSigmaAftCut(track, cent, isSpecies); - - for (int isp = 0; isp < KNsp; ++isp) { - if (!isSpecies[isp]) - continue; - - if (isp == kInclusiveIdx) { - histos.fill(HIST("h3_AllReco"), multPV, pt, eta); - if (track.has_mcParticle()) { - auto mcP = track.mcParticle(); - if (mcP.isPhysicalPrimary()) { - histos.fill(HIST("ptResolution"), mcP.pt(), (pt - mcP.pt()) / mcP.pt()); - histos.fill(HIST("etaResolution"), mcP.eta(), eta - mcP.eta()); - histos.fill(HIST("etaTruthReco"), mcP.eta(), eta); - histos.fill(HIST("vzResolution"), mcP.vz(), (vz - mcP.vz()) / mcP.vz()); - histos.fill(HIST("TruthTracKVz"), mcP.vz(), vz); - histos.fill(HIST("h3_RecoMatchedToPrimary"), multPV, mcP.pt(), mcP.eta()); - } else { - histos.fill(HIST("h3_RecoUnMatchedToPrimary_Secondary"), multPV, pt, eta); - } - } else { - histos.fill(HIST("h3_RecoUnMatchedToPrimary_Fake"), multPV, pt, eta); - } - } else if (isp == kPiMinusIdx) { - histos.fill(HIST("h3_AllReco_PiMinus"), multPV, pt, eta); - if (track.has_mcParticle()) { - auto mcP = track.mcParticle(); - if (mcP.isPhysicalPrimary()) { - if (mcP.pdgCode() == -KPiPlus) { - histos.fill(HIST("h3_RecoMatchedToPrimary_PiMinus"), multPV, mcP.pt(), mcP.eta()); - } else { // Misidentified - histos.fill(HIST("h3_RecoMatchedToPrimary_MisID_PiMinus"), multPV, pt, eta); - } - } else { - histos.fill(HIST("h3_RecoUnMatchedToPrimary_Secondary_PiMinus"), multPV, pt, eta); - } - } else { // No MC - histos.fill(HIST("h3_RecoUnMatchedToPrimary_Fake_PiMinus"), multPV, pt, eta); - } - } else if (isp == kPiPlusIdx) { - histos.fill(HIST("h3_AllReco_PiPlus"), multPV, pt, eta); - if (track.has_mcParticle()) { - auto mcP = track.mcParticle(); - if (mcP.isPhysicalPrimary()) { - if (mcP.pdgCode() == KPiPlus) { - histos.fill(HIST("h3_RecoMatchedToPrimary_PiPlus"), multPV, mcP.pt(), mcP.eta()); - } else { // Misidentified - histos.fill(HIST("h3_RecoMatchedToPrimary_MisID_PiPlus"), multPV, pt, eta); - } - } else { - histos.fill(HIST("h3_RecoUnMatchedToPrimary_Secondary_PiPlus"), multPV, pt, eta); - } - } else { - histos.fill(HIST("h3_RecoUnMatchedToPrimary_Fake_PiPlus"), multPV, pt, eta); - } - } else if (isp == kPiAllIdx) { - histos.fill(HIST("h3_AllReco_PiAll"), multPV, pt, eta); - if (track.has_mcParticle()) { - auto mcP = track.mcParticle(); - if (mcP.isPhysicalPrimary()) { - if (std::abs(mcP.pdgCode()) == KPiPlus) { - histos.fill(HIST("h3_RecoMatchedToPrimary_PiAll"), multPV, mcP.pt(), mcP.eta()); - } else { // Misidentified - histos.fill(HIST("h3_RecoMatchedToPrimary_MisID_PiAll"), multPV, pt, eta); - } - } else { - histos.fill(HIST("h3_RecoUnMatchedToPrimary_Secondary_PiAll"), multPV, pt, eta); - } - } else { - histos.fill(HIST("h3_RecoUnMatchedToPrimary_Fake_PiAll"), multPV, pt, eta); - } - } else if (isp == kKaMinusIdx) { - histos.fill(HIST("h3_AllReco_KaMinus"), multPV, pt, eta); - if (track.has_mcParticle()) { - auto mcP = track.mcParticle(); - if (mcP.isPhysicalPrimary()) { - if (mcP.pdgCode() == -KKPlus) { - histos.fill(HIST("h3_RecoMatchedToPrimary_KaMinus"), multPV, mcP.pt(), mcP.eta()); - } else { // Misidentified - histos.fill(HIST("h3_RecoMatchedToPrimary_MisID_KaMinus"), multPV, pt, eta); - } - } else { - histos.fill(HIST("h3_RecoUnMatchedToPrimary_Secondary_KaMinus"), multPV, pt, eta); - } - } else { - histos.fill(HIST("h3_RecoUnMatchedToPrimary_Fake_KaMinus"), multPV, pt, eta); - } - } else if (isp == kKaPlusIdx) { - histos.fill(HIST("h3_AllReco_KaPlus"), multPV, pt, eta); - if (track.has_mcParticle()) { - auto mcP = track.mcParticle(); - if (mcP.isPhysicalPrimary()) { - if (mcP.pdgCode() == KKPlus) { - histos.fill(HIST("h3_RecoMatchedToPrimary_KaPlus"), multPV, mcP.pt(), mcP.eta()); - } else { // Misidentified - histos.fill(HIST("h3_RecoMatchedToPrimary_MisID_KaPlus"), multPV, pt, eta); - } - } else { - histos.fill(HIST("h3_RecoUnMatchedToPrimary_Secondary_KaPlus"), multPV, pt, eta); - } - } else { - histos.fill(HIST("h3_RecoUnMatchedToPrimary_Fake_KaPlus"), multPV, pt, eta); - } - } else if (isp == kKaAllIdx) { - histos.fill(HIST("h3_AllReco_KaAll"), multPV, pt, eta); - if (track.has_mcParticle()) { - auto mcP = track.mcParticle(); - if (mcP.isPhysicalPrimary()) { - if (std::abs(mcP.pdgCode()) == KKPlus) { - histos.fill(HIST("h3_RecoMatchedToPrimary_KaAll"), multPV, mcP.pt(), mcP.eta()); - } else { // Misidentified - histos.fill(HIST("h3_RecoMatchedToPrimary_MisID_KaAll"), multPV, pt, eta); - } - } else { - histos.fill(HIST("h3_RecoUnMatchedToPrimary_Secondary_KaAll"), multPV, pt, eta); - } - } else { - histos.fill(HIST("h3_RecoUnMatchedToPrimary_Fake_KaAll"), multPV, pt, eta); - } - } else if (isp == kAntiPrIdx) { - histos.fill(HIST("h3_AllReco_AntiPr"), multPV, pt, eta); - if (track.has_mcParticle()) { - auto mcP = track.mcParticle(); - if (mcP.isPhysicalPrimary()) { - if (mcP.pdgCode() == -KProton) { - histos.fill(HIST("h3_RecoMatchedToPrimary_AntiPr"), multPV, mcP.pt(), mcP.eta()); - } else { // Misidentified - histos.fill(HIST("h3_RecoMatchedToPrimary_MisID_AntiPr"), multPV, pt, eta); - } - } else { - histos.fill(HIST("h3_RecoUnMatchedToPrimary_Secondary_AntiPr"), multPV, pt, eta); - } - } else { - histos.fill(HIST("h3_RecoUnMatchedToPrimary_Fake_AntiPr"), multPV, pt, eta); - } - } else if (isp == kPrIdx) { - histos.fill(HIST("h3_AllReco_Pr"), multPV, pt, eta); - if (track.has_mcParticle()) { - auto mcP = track.mcParticle(); - if (mcP.isPhysicalPrimary()) { - if (mcP.pdgCode() == KProton) { - histos.fill(HIST("h3_RecoMatchedToPrimary_Pr"), multPV, mcP.pt(), mcP.eta()); - } else { // Misidentified - histos.fill(HIST("h3_RecoMatchedToPrimary_MisID_Pr"), multPV, pt, eta); - } - } else { - histos.fill(HIST("h3_RecoUnMatchedToPrimary_Secondary_Pr"), multPV, pt, eta); - } - } else { - histos.fill(HIST("h3_RecoUnMatchedToPrimary_Fake_Pr"), multPV, pt, eta); - } - } else if (isp == kPrAllIdx) { - histos.fill(HIST("h3_AllReco_PrAll"), multPV, pt, eta); - if (track.has_mcParticle()) { - auto mcP = track.mcParticle(); - if (mcP.isPhysicalPrimary()) { - if (std::abs(mcP.pdgCode()) == KProton) { - histos.fill(HIST("h3_RecoMatchedToPrimary_PrAll"), multPV, mcP.pt(), mcP.eta()); - } else { // Misidentified - histos.fill(HIST("h3_RecoMatchedToPrimary_MisID_PrAll"), multPV, pt, eta); - } - } else { - histos.fill(HIST("h3_RecoUnMatchedToPrimary_Secondary_PrAll"), multPV, pt, eta); - } - } else { - histos.fill(HIST("h3_RecoUnMatchedToPrimary_Fake_PrAll"), multPV, pt, eta); - } - } } - } - } - PROCESS_SWITCH(RadialFlowDecorr, processGetEffHists, "process MC to calculate EffWeights", cfgRunGetEff); - - void processMCFlat(MyRun3MCCollisions::iterator const& mcCollision, FilteredTCs const& mcTracks) - { - histos.fill(HIST("hVtxZ"), mcCollision.posZ()); - if (!mcCollision.has_mcCollision() || !isEventSelected(mcCollision)) - return; - - float cent = getCentrality(mcCollision); - if (cent > KCentMax) - return; - - float multPV = mcCollision.multNTracksPV(); - float vz = mcCollision.posZ(); - - if (!isPassAddPileup(multPV, mcTracks.size(), cent)) - return; - histos.fill(HIST("hVtxZ_after_sel"), mcCollision.posZ()); - histos.fill(HIST("hCentrality"), cent); - - histos.fill(HIST("Hist2D_globalTracks_PVTracks"), mcCollision.multNTracksPV(), mcTracks.size()); - histos.fill(HIST("Hist2D_cent_nch"), mcTracks.size(), cent); - histos.fill(HIST("Hist2D_globalTracks_cent"), cent, mcTracks.size()); - histos.fill(HIST("Hist2D_PVTracks_cent"), cent, multPV); - for (const auto& track : mcTracks) { - if (track.collisionId() != mcCollision.index()) - continue; - - if (!isTrackSelected(track)) + if (!isTrackSelected(track)) { continue; + } float pt = track.pt(), eta = track.eta(), phi = track.phi(); auto sign = track.sign(); histos.fill(HIST("hPt"), pt); histos.fill(HIST("hEta"), eta); histos.fill(HIST("hPhi"), phi); - int id = identifyTrack(track, cent); - bool isPi = (id == KPidPionOne); - bool isKa = (id == KPidKaonTwo); - bool isPr = (id == KPidProtonThree); - bool isSpecies[KNsp] = { - true, - isPi && sign < 0, isPi && sign > 0, isPi, - isKa && sign < 0, isKa && sign > 0, isKa, - isPr && sign < 0, isPr && sign > 0, isPr}; - - for (int isp = 0; isp < KNsp; ++isp) { - if (!isSpecies[isp]) - continue; - float eff = getEfficiency(multPV, pt, eta, static_cast(isp), 0, cfgEff); - float fake = getEfficiency(multPV, pt, eta, static_cast(isp), 1, cfgEff); - float w = (eff > KFloatEpsilon) ? (1.0f - fake) / eff : 0.0f; - - if (std::isfinite(w) && w > 0.f) { - if (isp == kInclusiveIdx) { - histos.fill(HIST("MCReco/hEtaPhiRecoEffWtd"), vz, sign, pt, eta, phi, w); - histos.fill(HIST("MCReco/hEtaPhiReco"), vz, sign, pt, eta, phi, 1.0); - histos.fill(HIST("MCReco/hEtaPhiRecoWtd"), vz, sign, pt, eta, phi, w); - } else if (isp == kPiMinusIdx) { - histos.fill(HIST("MCReco/hEtaPhiRecoEffWtd_PiMinus"), vz, sign, pt, eta, phi, w); - histos.fill(HIST("MCReco/hEtaPhiReco_PiMinus"), vz, sign, pt, eta, phi, 1.0); - histos.fill(HIST("MCReco/hEtaPhiRecoWtd_PiMinus"), vz, sign, pt, eta, phi, w); - } else if (isp == kPiPlusIdx) { - histos.fill(HIST("MCReco/hEtaPhiRecoEffWtd_PiPlus"), vz, sign, pt, eta, phi, w); - histos.fill(HIST("MCReco/hEtaPhiReco_PiPlus"), vz, sign, pt, eta, phi, 1.0); - histos.fill(HIST("MCReco/hEtaPhiRecoWtd_PiPlus"), vz, sign, pt, eta, phi, w); - } else if (isp == kPiAllIdx) { - histos.fill(HIST("MCReco/hEtaPhiRecoEffWtd_PiAll"), vz, sign, pt, eta, phi, w); - histos.fill(HIST("MCReco/hEtaPhiReco_PiAll"), vz, sign, pt, eta, phi, 1.0); - histos.fill(HIST("MCReco/hEtaPhiRecoWtd_PiAll"), vz, sign, pt, eta, phi, w); - } else if (isp == kKaMinusIdx) { - histos.fill(HIST("MCReco/hEtaPhiRecoEffWtd_KaMinus"), vz, sign, pt, eta, phi, w); - histos.fill(HIST("MCReco/hEtaPhiReco_KaMinus"), vz, sign, pt, eta, phi, 1.0); - histos.fill(HIST("MCReco/hEtaPhiRecoWtd_KaMinus"), vz, sign, pt, eta, phi, w); - } else if (isp == kKaPlusIdx) { - histos.fill(HIST("MCReco/hEtaPhiRecoEffWtd_KaPlus"), vz, sign, pt, eta, phi, w); - histos.fill(HIST("MCReco/hEtaPhiReco_KaPlus"), vz, sign, pt, eta, phi, 1.0); - histos.fill(HIST("MCReco/hEtaPhiRecoWtd_KaPlus"), vz, sign, pt, eta, phi, w); - } else if (isp == kKaAllIdx) { - histos.fill(HIST("MCReco/hEtaPhiRecoEffWtd_KaAll"), vz, sign, pt, eta, phi, w); - histos.fill(HIST("MCReco/hEtaPhiReco_KaAll"), vz, sign, pt, eta, phi, 1.0); - histos.fill(HIST("MCReco/hEtaPhiRecoWtd_KaAll"), vz, sign, pt, eta, phi, w); - } else if (isp == kAntiPrIdx) { - histos.fill(HIST("MCReco/hEtaPhiRecoEffWtd_AntiPr"), vz, sign, pt, eta, phi, w); - histos.fill(HIST("MCReco/hEtaPhiReco_AntiPr"), vz, sign, pt, eta, phi, 1.0); - histos.fill(HIST("MCReco/hEtaPhiRecoWtd_AntiPr"), vz, sign, pt, eta, phi, w); - } else if (isp == kPrIdx) { - histos.fill(HIST("MCReco/hEtaPhiRecoEffWtd_Pr"), vz, sign, pt, eta, phi, w); - histos.fill(HIST("MCReco/hEtaPhiReco_Pr"), vz, sign, pt, eta, phi, 1.0); - histos.fill(HIST("MCReco/hEtaPhiRecoWtd_Pr"), vz, sign, pt, eta, phi, w); - } else if (isp == kPrAllIdx) { - histos.fill(HIST("MCReco/hEtaPhiRecoEffWtd_PrAll"), vz, sign, pt, eta, phi, w); - histos.fill(HIST("MCReco/hEtaPhiReco_PrAll"), vz, sign, pt, eta, phi, 1.0); - histos.fill(HIST("MCReco/hEtaPhiRecoWtd_PrAll"), vz, sign, pt, eta, phi, w); - } - } + float eff = getEfficiency(multPV, pt, eta, 0, cfgEff); + float fake = getEfficiency(multPV, pt, eta, 1, cfgEff); + float w = (eff > KFloatEpsilon) ? (1.0f - fake) / eff : 0.0f; + if (std::isfinite(w) && w > 0.f) { + histos.fill(HIST("MCReco/hEtaPhiRecoEffWtd"), vz, sign, pt, eta, phi, w); + histos.fill(HIST("MCReco/hEtaPhiReco"), vz, sign, pt, eta, phi, 1.0); + histos.fill(HIST("MCReco/hEtaPhiRecoWtd"), vz, sign, pt, eta, phi, w); } } } PROCESS_SWITCH(RadialFlowDecorr, processMCFlat, "process MC to calculate FlatWeights", cfgRunGetMCFlat); + // =========================================================================== + // MC: mean pT (truth / reco / reco-eff-corrected) + // =========================================================================== void processMCMean(MyRun3MCCollisions::iterator const& mcCollision, FilteredTCs const& mcTracks, aod::FT0s const&, aod::McParticles const& mcParticles) { - double sumWiTruth[KNsp][KNEta]{}, sumWiptiTruth[KNsp][KNEta]{}; - double sumWiReco[KNsp][KNEta]{}, sumWiptiReco[KNsp][KNEta]{}; - double sumWiRecoEffCorr[KNsp][KNEta]{}, sumWiptiRecoEffCorr[KNsp][KNEta]{}; + std::array sumWiTruth{}, sumWiptiTruth{}; + std::array sumWiReco{}, sumWiptiReco{}; + std::array sumWiRecoEffCorr{}, sumWiptiRecoEffCorr{}; + histos.fill(HIST("hVtxZ"), mcCollision.posZ()); - if (!mcCollision.has_mcCollision() || !isEventSelected(mcCollision)) + if (!mcCollision.has_mcCollision() || !isEventSelected(mcCollision)) { return; + } float cent = getCentrality(mcCollision); - if (cent > KCentMax) + if (cent > KCentMax) { return; + } float multPV = mcCollision.multNTracksPV(); float vz = mcCollision.posZ(); - if (!isPassAddPileup(multPV, mcTracks.size(), cent)) + if (!isPassAddPileup(multPV, mcTracks.size(), cent)) { return; + } histos.fill(HIST("hVtxZ_after_sel"), mcCollision.posZ()); histos.fill(HIST("hCentrality"), cent); - - histos.fill(HIST("Hist2D_globalTracks_PVTracks"), mcCollision.multNTracksPV(), mcTracks.size()); + histos.fill(HIST("Hist2D_globalTracks_PVTracks"), multPV, mcTracks.size()); histos.fill(HIST("Hist2D_cent_nch"), mcTracks.size(), cent); histos.fill(HIST("Hist2D_globalTracks_cent"), cent, mcTracks.size()); histos.fill(HIST("Hist2D_PVTracks_cent"), cent, multPV); - memset(sumWiTruth, 0, sizeof(sumWiTruth)); - memset(sumWiptiTruth, 0, sizeof(sumWiptiTruth)); - memset(sumWiReco, 0, sizeof(sumWiReco)); - memset(sumWiptiReco, 0, sizeof(sumWiptiReco)); - memset(sumWiRecoEffCorr, 0, sizeof(sumWiRecoEffCorr)); - memset(sumWiptiRecoEffCorr, 0, sizeof(sumWiptiRecoEffCorr)); - + // --- truth --- for (const auto& particle : mcParticles) { - if (particle.mcCollisionId() != mcCollision.mcCollisionId()) + if (particle.mcCollisionId() != mcCollision.mcCollisionId()) { continue; - - if (!isParticleSelected(particle) || !particle.isPhysicalPrimary()) + } + if (!isParticleSelected(particle) || !particle.isPhysicalPrimary()) { continue; + } float pt = particle.pt(), eta = particle.eta(); - if (pt <= cfgPtMin || pt > cfgPtMax) + if (pt <= cfgPtMin || pt > cfgPtMax) { continue; - int pdgCode = particle.pdgCode(); - int absPdg = std::abs(pdgCode); - - bool isSpecies[KNsp] = { - true, // kInclusiveIdx - pdgCode == -KPiPlus, // kPiMinusIdx - pdgCode == KPiPlus, // kPiPlusIdx - absPdg == KPiPlus, // kPiAllIdx - pdgCode == -KKPlus, // kKaMinusIdx - pdgCode == KKPlus, // kKaPlusIdx - absPdg == KKPlus, // kKaAllIdx - pdgCode == -KProton, // kAntiPrIdx - pdgCode == KProton, // kPrIdx - absPdg == KProton // kPrAllIdx - }; - - for (int ieta = 0; ieta < KNEta; ++ieta) { - if (eta <= etaLw[ieta] || eta > etaUp[ieta]) + } + for (int ieta = 0; ieta < nEta; ++ieta) { + if (eta <= etaLw[ieta] || eta > etaUp[ieta]) { continue; - - for (int isp = 0; isp < KNsp; ++isp) { - if (isSpecies[isp]) { - sumWiTruth[isp][ieta]++; - sumWiptiTruth[isp][ieta] += pt; - } } + sumWiTruth[ieta]++; + sumWiptiTruth[ieta] += pt; } } - for (int isp = 0; isp < KNsp; ++isp) { - histos.fill(HIST("MCGen/Prof_Cent_Nsp_Nchrec"), cent, isp, sumWiTruth[isp][0]); - histos.fill(HIST("MCGen/Prof_Mult_Nsp_Nchrec"), multPV, isp, sumWiTruth[isp][0]); - if (sumWiTruth[isp][0] > 1.0f) { - histos.fill(HIST("MCGen/Prof_Cent_Nsp_MeanpT"), cent, isp, sumWiptiTruth[isp][0] / sumWiTruth[isp][0]); - histos.fill(HIST("MCGen/Prof_Mult_Nsp_MeanpT"), multPV, isp, sumWiptiTruth[isp][0] / sumWiTruth[isp][0]); - } + histos.fill(HIST("MCGen/Prof_Cent_Nchrec"), cent, sumWiTruth[0]); + histos.fill(HIST("MCGen/Prof_Mult_Nchrec"), multPV, sumWiTruth[0]); + if (sumWiTruth[0] > 1.0f) { + histos.fill(HIST("MCGen/Prof_Cent_MeanpT"), cent, sumWiptiTruth[0] / sumWiTruth[0]); + histos.fill(HIST("MCGen/Prof_Mult_MeanpT"), multPV, sumWiptiTruth[0] / sumWiTruth[0]); } + // --- reco --- for (const auto& track : mcTracks) { - if (track.collisionId() != mcCollision.index()) + if (track.collisionId() != mcCollision.index()) { continue; - - if (!isTrackSelected(track)) + } + if (!isTrackSelected(track)) { continue; + } float pt = track.pt(), eta = track.eta(), phi = track.phi(); - if (pt <= cfgPtMin || pt > cfgPtMax) + if (pt <= cfgPtMin || pt > cfgPtMax) { continue; + } auto sign = track.sign(); histos.fill(HIST("hPt"), pt); histos.fill(HIST("hEta"), eta); histos.fill(HIST("hPhi"), phi); - int id = identifyTrack(track, cent); - bool isPi = (id == KPidPionOne); - bool isKa = (id == KPidKaonTwo); - bool isPr = (id == KPidProtonThree); - bool isSpecies[KNsp] = { - true, - isPi && sign < 0, isPi && sign > 0, isPi, - isKa && sign < 0, isKa && sign > 0, isKa, - isPr && sign < 0, isPr && sign > 0, isPr}; - - for (int isp = 0; isp < KNsp; ++isp) { - if (!isSpecies[isp]) - continue; - float eff = getEfficiency(multPV, pt, eta, static_cast(isp), 0, cfgEff); - float fake = getEfficiency(multPV, pt, eta, static_cast(isp), 1, cfgEff); - float flatW = getFlatteningWeight(vz, sign, pt, eta, phi, static_cast(isp), cfgFlat); - float w = flatW * (1.0 - fake) / eff; - if (!std::isfinite(w) || w <= 0.f || eff <= KFloatEpsilon) - continue; - - for (int ieta = 0; ieta < KNEta; ++ieta) { - if (eta <= etaLw[ieta] || eta > etaUp[ieta]) - continue; - sumWiReco[isp][ieta]++; - sumWiptiReco[isp][ieta] += pt; - sumWiRecoEffCorr[isp][ieta] += w; - sumWiptiRecoEffCorr[isp][ieta] += w * pt; - } - - if (isp == kInclusiveIdx) { - histos.fill(HIST("Eff_cent"), cent, eff); - histos.fill(HIST("Fake_cent"), cent, fake); - histos.fill(HIST("wgt_cent"), cent, w); - - histos.fill(HIST("Eff_Ntrk"), multPV, eff); - histos.fill(HIST("Fake_Ntrk"), multPV, fake); - histos.fill(HIST("wgt_Ntrk"), multPV, w); - histos.fill(HIST("Eff_pT"), pt, eff); - histos.fill(HIST("Fake_pT"), pt, fake); - histos.fill(HIST("wgt_pT"), pt, w); + float eff = getEfficiency(multPV, pt, eta, 0, cfgEff); + float fake = getEfficiency(multPV, pt, eta, 1, cfgEff); + float flatW = getFlatteningWeight(vz, sign, pt, eta, phi, cfgFlat); + float w = flatW * (1.0 - fake) / eff; + if (!std::isfinite(w) || w <= 0.f || eff <= KFloatEpsilon) { + continue; + } - histos.fill(HIST("Eff_eta"), eta, eff); - histos.fill(HIST("Fake_eta"), eta, fake); - histos.fill(HIST("wgt_eta"), eta, w); - } - if (isp == kInclusiveIdx) { - histos.fill(HIST("hEtaPhiReco"), vz, sign, pt, eta, phi); - histos.fill(HIST("hEtaPhiRecoWtd"), vz, sign, pt, eta, phi, w); - histos.fill(HIST("hEtaPhiRecoEffWtd"), vz, sign, pt, eta, phi, (1.0 - fake) / eff); - } else if (isp == kPiMinusIdx) { - histos.fill(HIST("hEtaPhiReco_PiMinus"), vz, sign, pt, eta, phi); - histos.fill(HIST("hEtaPhiRecoWtd_PiMinus"), vz, sign, pt, eta, phi, w); - histos.fill(HIST("hEtaPhiRecoEffWtd_PiMinus"), vz, sign, pt, eta, phi, (1.0 - fake) / eff); - } else if (isp == kPiPlusIdx) { - histos.fill(HIST("hEtaPhiReco_PiPlus"), vz, sign, pt, eta, phi); - histos.fill(HIST("hEtaPhiRecoWtd_PiPlus"), vz, sign, pt, eta, phi, w); - histos.fill(HIST("hEtaPhiRecoEffWtd_PiPlus"), vz, sign, pt, eta, phi, (1.0 - fake) / eff); - } else if (isp == kPiAllIdx) { - histos.fill(HIST("hEtaPhiReco_PiAll"), vz, sign, pt, eta, phi); - histos.fill(HIST("hEtaPhiRecoWtd_PiAll"), vz, sign, pt, eta, phi, w); - histos.fill(HIST("hEtaPhiRecoEffWtd_PiAll"), vz, sign, pt, eta, phi, (1.0 - fake) / eff); - } else if (isp == kKaMinusIdx) { - histos.fill(HIST("hEtaPhiReco_KaMinus"), vz, sign, pt, eta, phi); - histos.fill(HIST("hEtaPhiRecoWtd_KaMinus"), vz, sign, pt, eta, phi, w); - histos.fill(HIST("hEtaPhiRecoEffWtd_KaMinus"), vz, sign, pt, eta, phi, (1.0 - fake) / eff); - } else if (isp == kKaPlusIdx) { - histos.fill(HIST("hEtaPhiReco_KaPlus"), vz, sign, pt, eta, phi); - histos.fill(HIST("hEtaPhiRecoWtd_KaPlus"), vz, sign, pt, eta, phi, w); - histos.fill(HIST("hEtaPhiRecoEffWtd_KaPlus"), vz, sign, pt, eta, phi, (1.0 - fake) / eff); - } else if (isp == kKaAllIdx) { - histos.fill(HIST("hEtaPhiReco_KaAll"), vz, sign, pt, eta, phi); - histos.fill(HIST("hEtaPhiRecoWtd_KaAll"), vz, sign, pt, eta, phi, w); - histos.fill(HIST("hEtaPhiRecoEffWtd_KaAll"), vz, sign, pt, eta, phi, (1.0 - fake) / eff); - } else if (isp == kPrIdx) { - histos.fill(HIST("hEtaPhiReco_Pr"), vz, sign, pt, eta, phi); - histos.fill(HIST("hEtaPhiRecoWtd_Pr"), vz, sign, pt, eta, phi, w); - histos.fill(HIST("hEtaPhiRecoEffWtd_Pr"), vz, sign, pt, eta, phi, (1.0 - fake) / eff); - } else if (isp == kAntiPrIdx) { - histos.fill(HIST("hEtaPhiReco_AntiPr"), vz, sign, pt, eta, phi); - histos.fill(HIST("hEtaPhiRecoWtd_AntiPr"), vz, sign, pt, eta, phi, w); - histos.fill(HIST("hEtaPhiRecoEffWtd_AntiPr"), vz, sign, pt, eta, phi, (1.0 - fake) / eff); - } else if (isp == kPrAllIdx) { - histos.fill(HIST("hEtaPhiReco_PrAll"), vz, sign, pt, eta, phi); - histos.fill(HIST("hEtaPhiRecoWtd_PrAll"), vz, sign, pt, eta, phi, w); - histos.fill(HIST("hEtaPhiRecoEffWtd_PrAll"), vz, sign, pt, eta, phi, (1.0 - fake) / eff); + for (int ieta = 0; ieta < nEta; ++ieta) { + if (eta <= etaLw[ieta] || eta > etaUp[ieta]) { + continue; } + sumWiReco[ieta]++; + sumWiptiReco[ieta] += pt; + sumWiRecoEffCorr[ieta] += w; + sumWiptiRecoEffCorr[ieta] += w * pt; } - } - if (!hasMinTracksInAllEtaBins(sumWiTruth) || !hasMinTracksInAllEtaBins(sumWiReco)) + histos.fill(HIST("Eff_cent"), cent, eff); + histos.fill(HIST("Fake_cent"), cent, fake); + histos.fill(HIST("wgt_cent"), cent, w); + histos.fill(HIST("Eff_Ntrk"), multPV, eff); + histos.fill(HIST("Fake_Ntrk"), multPV, fake); + histos.fill(HIST("wgt_Ntrk"), multPV, w); + histos.fill(HIST("Eff_pT"), pt, eff); + histos.fill(HIST("Fake_pT"), pt, fake); + histos.fill(HIST("wgt_pT"), pt, w); + histos.fill(HIST("Eff_eta"), eta, eff); + histos.fill(HIST("Fake_eta"), eta, fake); + histos.fill(HIST("wgt_eta"), eta, w); + + histos.fill(HIST("hEtaPhiReco"), vz, sign, pt, eta, phi); + histos.fill(HIST("hEtaPhiRecoWtd"), vz, sign, pt, eta, phi, w); + histos.fill(HIST("hEtaPhiRecoEffWtd"), vz, sign, pt, eta, phi, (1.0 - fake) / eff); + } + + if (!hasMinTracksInAllEtaBins(sumWiTruth) || !hasMinTracksInAllEtaBins(sumWiReco)) { return; + } - for (int isp = 0; isp < KNsp; ++isp) { - histos.fill(HIST("MCReco/Prof_Cent_Nsp_Nchrec"), cent, isp, sumWiReco[isp][0]); - histos.fill(HIST("MCReco/Prof_Mult_Nsp_Nchrec"), multPV, isp, sumWiReco[isp][0]); - - histos.fill(HIST("MCRecoEffCorr/Prof_Cent_Nsp_Nchrec"), cent, isp, sumWiRecoEffCorr[isp][0]); - histos.fill(HIST("MCRecoEffCorr/Prof_Mult_Nsp_Nchrec"), multPV, isp, sumWiRecoEffCorr[isp][0]); - - if (sumWiReco[isp][0] > 1.0f) { - histos.fill(HIST("MCReco/Prof_Cent_Nsp_MeanpT"), cent, isp, sumWiptiReco[isp][0] / sumWiReco[isp][0]); - histos.fill(HIST("MCReco/Prof_Mult_Nsp_MeanpT"), multPV, isp, sumWiptiReco[isp][0] / sumWiReco[isp][0]); - } - if (sumWiRecoEffCorr[isp][0] > 1.0f) { - histos.fill(HIST("MCRecoEffCorr/Prof_Cent_Nsp_MeanpT"), cent, isp, sumWiptiRecoEffCorr[isp][0] / sumWiRecoEffCorr[isp][0]); - histos.fill(HIST("MCRecoEffCorr/Prof_Mult_Nsp_MeanpT"), multPV, isp, sumWiptiRecoEffCorr[isp][0] / sumWiRecoEffCorr[isp][0]); - } - } - - for (int ietaA = 0; ietaA < KNEta; ++ietaA) { - for (int ietaC = 0; ietaC < KNEta; ++ietaC) { - for (int isp = 0; isp < KNsp; ++isp) { - float nTruAB = sumWiTruth[isp][ietaA] + sumWiTruth[isp][ietaC]; - float nRecoAB = sumWiReco[isp][ietaA] + sumWiReco[isp][ietaC]; - float nCorrAB = sumWiRecoEffCorr[isp][ietaA] + sumWiRecoEffCorr[isp][ietaC]; - - float mptsubTru = (sumWiptiTruth[isp][ietaA] + sumWiptiTruth[isp][ietaC]) / nTruAB; - float mptsubReco = (sumWiptiReco[isp][ietaA] + sumWiptiReco[isp][ietaC]) / nRecoAB; - float mptsubRecoEffCorr = (sumWiptiRecoEffCorr[isp][ietaA] + sumWiptiRecoEffCorr[isp][ietaC]) / nCorrAB; - - if (nTruAB > 0) { - if (isp == kInclusiveIdx) - histos.fill(HIST("Prof2D_MeanpTSub_Tru"), cent, ietaA, ietaC, mptsubTru); - else if (isp == kPiMinusIdx) - histos.fill(HIST("Prof2D_MeanpTSub_Tru_PiMinus"), cent, ietaA, ietaC, mptsubTru); - else if (isp == kPiPlusIdx) - histos.fill(HIST("Prof2D_MeanpTSub_Tru_PiPlus"), cent, ietaA, ietaC, mptsubTru); - else if (isp == kPiAllIdx) - histos.fill(HIST("Prof2D_MeanpTSub_Tru_PiAll"), cent, ietaA, ietaC, mptsubTru); - else if (isp == kKaMinusIdx) - histos.fill(HIST("Prof2D_MeanpTSub_Tru_KaMinus"), cent, ietaA, ietaC, mptsubTru); - else if (isp == kKaPlusIdx) - histos.fill(HIST("Prof2D_MeanpTSub_Tru_KaPlus"), cent, ietaA, ietaC, mptsubTru); - else if (isp == kKaAllIdx) - histos.fill(HIST("Prof2D_MeanpTSub_Tru_KaAll"), cent, ietaA, ietaC, mptsubTru); - else if (isp == kPrIdx) - histos.fill(HIST("Prof2D_MeanpTSub_Tru_Pr"), cent, ietaA, ietaC, mptsubTru); - else if (isp == kAntiPrIdx) - histos.fill(HIST("Prof2D_MeanpTSub_Tru_AntiPr"), cent, ietaA, ietaC, mptsubTru); - else if (isp == kPrAllIdx) - histos.fill(HIST("Prof2D_MeanpTSub_Tru_PrAll"), cent, ietaA, ietaC, mptsubTru); - } - - if (nRecoAB > 0) { - if (isp == kInclusiveIdx) - histos.fill(HIST("Prof2D_MeanpTSub_Reco"), cent, ietaA, ietaC, mptsubReco); - else if (isp == kPiMinusIdx) - histos.fill(HIST("Prof2D_MeanpTSub_Reco_PiMinus"), cent, ietaA, ietaC, mptsubReco); - else if (isp == kPiPlusIdx) - histos.fill(HIST("Prof2D_MeanpTSub_Reco_PiPlus"), cent, ietaA, ietaC, mptsubReco); - else if (isp == kPiAllIdx) - histos.fill(HIST("Prof2D_MeanpTSub_Reco_PiAll"), cent, ietaA, ietaC, mptsubReco); - else if (isp == kKaMinusIdx) - histos.fill(HIST("Prof2D_MeanpTSub_Reco_KaMinus"), cent, ietaA, ietaC, mptsubReco); - else if (isp == kKaPlusIdx) - histos.fill(HIST("Prof2D_MeanpTSub_Reco_KaPlus"), cent, ietaA, ietaC, mptsubReco); - else if (isp == kKaAllIdx) - histos.fill(HIST("Prof2D_MeanpTSub_Reco_KaAll"), cent, ietaA, ietaC, mptsubReco); - else if (isp == kPrIdx) - histos.fill(HIST("Prof2D_MeanpTSub_Reco_Pr"), cent, ietaA, ietaC, mptsubReco); - else if (isp == kAntiPrIdx) - histos.fill(HIST("Prof2D_MeanpTSub_Reco_AntiPr"), cent, ietaA, ietaC, mptsubReco); - else if (isp == kPrAllIdx) - histos.fill(HIST("Prof2D_MeanpTSub_Reco_PrAll"), cent, ietaA, ietaC, mptsubReco); - } - - if (nCorrAB > 0) { - if (isp == kInclusiveIdx) - histos.fill(HIST("Prof2D_MeanpTSub_RecoEffCorr"), cent, ietaA, ietaC, mptsubRecoEffCorr); - else if (isp == kPiMinusIdx) - histos.fill(HIST("Prof2D_MeanpTSub_RecoEffCorr_PiMinus"), cent, ietaA, ietaC, mptsubRecoEffCorr); - else if (isp == kPiPlusIdx) - histos.fill(HIST("Prof2D_MeanpTSub_RecoEffCorr_PiPlus"), cent, ietaA, ietaC, mptsubRecoEffCorr); - else if (isp == kPiAllIdx) - histos.fill(HIST("Prof2D_MeanpTSub_RecoEffCorr_PiAll"), cent, ietaA, ietaC, mptsubRecoEffCorr); - else if (isp == kKaMinusIdx) - histos.fill(HIST("Prof2D_MeanpTSub_RecoEffCorr_KaMinus"), cent, ietaA, ietaC, mptsubRecoEffCorr); - else if (isp == kKaPlusIdx) - histos.fill(HIST("Prof2D_MeanpTSub_RecoEffCorr_KaPlus"), cent, ietaA, ietaC, mptsubRecoEffCorr); - else if (isp == kKaAllIdx) - histos.fill(HIST("Prof2D_MeanpTSub_RecoEffCorr_KaAll"), cent, ietaA, ietaC, mptsubRecoEffCorr); - else if (isp == kPrIdx) - histos.fill(HIST("Prof2D_MeanpTSub_RecoEffCorr_Pr"), cent, ietaA, ietaC, mptsubRecoEffCorr); - else if (isp == kAntiPrIdx) - histos.fill(HIST("Prof2D_MeanpTSub_RecoEffCorr_AntiPr"), cent, ietaA, ietaC, mptsubRecoEffCorr); - else if (isp == kPrAllIdx) - histos.fill(HIST("Prof2D_MeanpTSub_RecoEffCorr_PrAll"), cent, ietaA, ietaC, mptsubRecoEffCorr); - } - } - } + // subevent mean-pT maps + for (int ietaA = 0; ietaA < nEta; ++ietaA) { + for (int ietaC = 0; ietaC < nEta; ++ietaC) { + float nTruAB = sumWiTruth[ietaA] + sumWiTruth[ietaC]; + float nRecoAB = sumWiReco[ietaA] + sumWiReco[ietaC]; + float nCorrAB = sumWiRecoEffCorr[ietaA] + sumWiRecoEffCorr[ietaC]; - for (int isp = 0; isp < KNsp; ++isp) { - if (sumWiTruth[isp][ietaA] > 0) { - float val = sumWiptiTruth[isp][ietaA] / sumWiTruth[isp][ietaA]; - histos.fill(HIST("pmeanTru_nch_etabin_spbin"), multPV, ietaA, isp, val); - histos.fill(HIST("pmeanMultTru_nch_etabin_spbin"), multPV, ietaA, isp, sumWiTruth[isp][ietaA]); + if (nTruAB > 0) { + histos.fill(HIST("Prof2D_MeanpTSub_Tru"), cent, ietaA, ietaC, (sumWiptiTruth[ietaA] + sumWiptiTruth[ietaC]) / nTruAB); } - if (sumWiReco[isp][ietaA] > 0) { - float val = sumWiptiReco[isp][ietaA] / sumWiReco[isp][ietaA]; - histos.fill(HIST("pmeanReco_nch_etabin_spbin"), multPV, ietaA, isp, val); - histos.fill(HIST("pmeanMultReco_nch_etabin_spbin"), multPV, ietaA, isp, sumWiReco[isp][ietaA]); + if (nRecoAB > 0) { + histos.fill(HIST("Prof2D_MeanpTSub_Reco"), cent, ietaA, ietaC, (sumWiptiReco[ietaA] + sumWiptiReco[ietaC]) / nRecoAB); } - if (sumWiRecoEffCorr[isp][ietaA] > 0) { - float val = sumWiptiRecoEffCorr[isp][ietaA] / sumWiRecoEffCorr[isp][ietaA]; - histos.fill(HIST("pmeanRecoEffcorr_nch_etabin_spbin"), multPV, ietaA, isp, val); - histos.fill(HIST("pmeanMultRecoEffcorr_nch_etabin_spbin"), multPV, ietaA, isp, sumWiRecoEffCorr[isp][ietaA]); + if (nCorrAB > 0) { + histos.fill(HIST("Prof2D_MeanpTSub_RecoEffCorr"), cent, ietaA, ietaC, (sumWiptiRecoEffCorr[ietaA] + sumWiptiRecoEffCorr[ietaC]) / nCorrAB); } } - } // end ietaA + if (sumWiTruth[ietaA] > 0) { + histos.fill(HIST("pmeanTru_nch_etabin"), multPV, ietaA, sumWiptiTruth[ietaA] / sumWiTruth[ietaA]); + histos.fill(HIST("pmeanMultTru_nch_etabin"), multPV, ietaA, sumWiTruth[ietaA]); + } + if (sumWiReco[ietaA] > 0) { + histos.fill(HIST("pmeanReco_nch_etabin"), multPV, ietaA, sumWiptiReco[ietaA] / sumWiReco[ietaA]); + histos.fill(HIST("pmeanMultReco_nch_etabin"), multPV, ietaA, sumWiReco[ietaA]); + } + if (sumWiRecoEffCorr[ietaA] > 0) { + histos.fill(HIST("pmeanRecoEffcorr_nch_etabin"), multPV, ietaA, sumWiptiRecoEffCorr[ietaA] / sumWiRecoEffCorr[ietaA]); + histos.fill(HIST("pmeanMultRecoEffcorr_nch_etabin"), multPV, ietaA, sumWiRecoEffCorr[ietaA]); + } + } + + // FT0 double amplFT0A = 0, amplFT0C = 0; if (mcCollision.has_foundFT0()) { const auto& ft0 = mcCollision.foundFT0(); @@ -2161,265 +1451,180 @@ struct RadialFlowDecorr { histos.fill(HIST("h3_cent_id_eta_FT0"), cent, globalId, eta, ampl); } } - histos.fill(HIST("pmeanFT0Amultpv"), multPV, amplFT0A); histos.fill(HIST("pmeanFT0A_cent"), cent, amplFT0A); histos.fill(HIST("pmeanFT0Cmultpv"), multPV, amplFT0C); histos.fill(HIST("pmeanFT0C_cent"), cent, amplFT0C); } - PROCESS_SWITCH(RadialFlowDecorr, processMCMean, "process MC to calculate mean pt and Eff Hists", cfgRunMCMean); + PROCESS_SWITCH(RadialFlowDecorr, processMCMean, "process MC to calculate mean pt", cfgRunMCMean); + // =========================================================================== + // MC: fluctuations (C2, subevent) at three levels + // =========================================================================== void processMCFluc(MyRun3MCCollisions::iterator const& mcCollision, FilteredTCs const& mcTracks, aod::FT0s const&, aod::McParticles const& mcParticles) { - if (!state.pmeanTruNchEtabinSpbinStep2 || !state.pmeanRecoNchEtabinSpbinStep2 || !state.pmeanRecoEffcorrNchEtabinSpbinStep2 || - !state.pmeanMultTruNchEtabinSpbinStep2 || !state.pmeanMultRecoNchEtabinSpbinStep2 || !state.pmeanMultRecoEffcorrNchEtabinSpbinStep2) { - LOGF(warning, "MC fluc: Unified Mean pT or Mult map missing"); + if (!state.pmeanTruNchEtabinStep2 || !state.pmeanRecoNchEtabinStep2 || !state.pmeanRecoEffcorrNchEtabinStep2 || + !state.pmeanMultTruNchEtabinStep2 || !state.pmeanMultRecoNchEtabinStep2 || !state.pmeanMultRecoEffcorrNchEtabinStep2) { + LOGF(warning, "MC fluc: mean pT or mult map missing"); return; } - double sumPmwkTru[KNsp][KNEta][KIntM][KIntK]{}; - double sumWkTru[KNsp][KNEta][KIntK]{}; - double sumPmwkReco[KNsp][KNEta][KIntM][KIntK]{}; - double sumWkReco[KNsp][KNEta][KIntK]{}; - double sumPmwkRecoEffCor[KNsp][KNEta][KIntM][KIntK]{}; - double sumWkRecoEffCor[KNsp][KNEta][KIntK]{}; - double meanTru[KNsp][KNEta]{}, c2Tru[KNsp][KNEta]{}; - double meanReco[KNsp][KNEta]{}, c2Reco[KNsp][KNEta]{}; - double meanRecoEffCor[KNsp][KNEta]{}, c2RecoEffCor[KNsp][KNEta]{}; + std::array, KIntM>, KNEtaMax> sumPmwkTru{}; + std::array, KNEtaMax> sumWkTru{}; + std::array, KIntM>, KNEtaMax> sumPmwkReco{}; + std::array, KNEtaMax> sumWkReco{}; + std::array, KIntM>, KNEtaMax> sumPmwkRecoEffCor{}; + std::array, KNEtaMax> sumWkRecoEffCor{}; - double meanTruMult[KNsp][KNEta]{}; - double meanRecoMult[KNsp][KNEta]{}; - double meanRecoEffCorMult[KNsp][KNEta]{}; + std::array meanTru{}, c2Tru{}; + std::array meanReco{}, c2Reco{}; + std::array meanRecoEffCor{}, c2RecoEffCor{}; - double p1kBarTru[KNsp][KNEta]{}, p1kBarReco[KNsp][KNEta]{}, p1kBarRecoEffCor[KNsp][KNEta]{}; - double p1kBarTruMult[KNsp][KNEta]{}, p1kBarRecoMult[KNsp][KNEta]{}, p1kBarRecoEffCorMult[KNsp][KNEta]{}; + std::array meanTruMult{}, meanRecoMult{}, meanRecoEffCorMult{}; + std::array p1kBarTru{}, p1kBarReco{}, p1kBarRecoEffCor{}; + std::array p1kBarTruMult{}, p1kBarRecoMult{}, p1kBarRecoEffCorMult{}; - if (!mcCollision.has_mcCollision() || !isEventSelected(mcCollision)) + if (!mcCollision.has_mcCollision() || !isEventSelected(mcCollision)) { return; + } float cent = getCentrality(mcCollision); - if (cent > KCentMax) + if (cent > KCentMax) { return; + } float multPV = mcCollision.multNTracksPV(); float vz = mcCollision.posZ(); - if (!isPassAddPileup(multPV, mcTracks.size(), cent)) + if (!isPassAddPileup(multPV, mcTracks.size(), cent)) { return; + } histos.fill(HIST("hVtxZ_after_sel"), mcCollision.posZ()); histos.fill(HIST("hCentrality"), cent); - histos.fill(HIST("Hist2D_globalTracks_PVTracks"), multPV, mcTracks.size()); histos.fill(HIST("Hist2D_cent_nch"), mcTracks.size(), cent); histos.fill(HIST("Hist2D_globalTracks_cent"), cent, mcTracks.size()); histos.fill(HIST("Hist2D_PVTracks_cent"), cent, multPV); - memset(sumPmwkTru, 0, sizeof(sumPmwkTru)); - memset(sumWkTru, 0, sizeof(sumWkTru)); - memset(sumPmwkReco, 0, sizeof(sumPmwkReco)); - memset(sumWkReco, 0, sizeof(sumWkReco)); - memset(sumPmwkRecoEffCor, 0, sizeof(sumPmwkRecoEffCor)); - memset(sumWkRecoEffCor, 0, sizeof(sumWkRecoEffCor)); - - memset(meanTru, 0, sizeof(meanTru)); - memset(c2Tru, 0, sizeof(c2Tru)); - memset(meanReco, 0, sizeof(meanReco)); - memset(c2Reco, 0, sizeof(c2Reco)); - memset(meanRecoEffCor, 0, sizeof(meanRecoEffCor)); - memset(c2RecoEffCor, 0, sizeof(c2RecoEffCor)); - - memset(meanTruMult, 0, sizeof(meanTruMult)); - memset(meanRecoMult, 0, sizeof(meanRecoMult)); - memset(meanRecoEffCorMult, 0, sizeof(meanRecoEffCorMult)); - - memset(p1kBarTru, 0, sizeof(p1kBarTru)); - memset(p1kBarReco, 0, sizeof(p1kBarReco)); - memset(p1kBarRecoEffCor, 0, sizeof(p1kBarRecoEffCor)); - - memset(p1kBarTruMult, 0, sizeof(p1kBarTruMult)); - memset(p1kBarRecoMult, 0, sizeof(p1kBarRecoMult)); - memset(p1kBarRecoEffCorMult, 0, sizeof(p1kBarRecoEffCorMult)); - double p1kBarFt0A = 0.0, p1kBarFt0C = 0.0; + // --- truth sums --- for (const auto& particle : mcParticles) { - if (particle.mcCollisionId() != mcCollision.mcCollisionId()) + if (particle.mcCollisionId() != mcCollision.mcCollisionId()) { continue; - - if (!isParticleSelected(particle) || !particle.isPhysicalPrimary()) + } + if (!isParticleSelected(particle) || !particle.isPhysicalPrimary()) { continue; - + } float pt = particle.pt(); - if (pt <= cfgPtMin || pt > cfgPtMax) + if (pt <= cfgPtMin || pt > cfgPtMax) { continue; + } float eta = particle.eta(); - int pdgCode = particle.pdgCode(); - int absPdg = std::abs(pdgCode); - - bool isSpecies[KNsp] = { - true, // kInclusiveIdx - pdgCode == -KPiPlus, // kPiMinusIdx - pdgCode == KPiPlus, // kPiPlusIdx - absPdg == KPiPlus, // kPiAllIdx - pdgCode == -KKPlus, // kKaMinusIdx - pdgCode == KKPlus, // kKaPlusIdx - absPdg == KKPlus, // kKaAllIdx - pdgCode == -KProton, // kAntiPrIdx - pdgCode == KProton, // kPrIdx - absPdg == KProton // kPrAllIdx - }; - - for (int ieta = 0; ieta < KNEta; ++ieta) { - if (eta <= etaLw[ieta] || eta > etaUp[ieta]) + for (int ieta = 0; ieta < nEta; ++ieta) { + if (eta <= etaLw[ieta] || eta > etaUp[ieta]) { continue; - for (int isp = 0; isp < KNsp; ++isp) { - if (isSpecies[isp]) { - for (int k = 0; k < KIntK; ++k) { - for (int m = 0; m < KIntM; ++m) { - sumPmwkTru[isp][ieta][m][k] += std::pow(pt, m); - } - sumWkTru[isp][ieta][k]++; - } + } + for (int k = 0; k < KIntK; ++k) { + for (int m = 0; m < KIntM; ++m) { + sumPmwkTru[ieta][m][k] += std::pow(pt, m); } + sumWkTru[ieta][k]++; } } - } // end truth loop + } + // --- reco sums --- for (const auto& track : mcTracks) { - if (track.collisionId() != mcCollision.index()) + if (track.collisionId() != mcCollision.index()) { continue; - - if (!isTrackSelected(track)) + } + if (!isTrackSelected(track)) { continue; - + } float pt = track.pt(); - if (pt <= cfgPtMin || pt > cfgPtMax) + if (pt <= cfgPtMin || pt > cfgPtMax) { continue; - float eta = track.eta(); - float phi = track.phi(); + } + float eta = track.eta(), phi = track.phi(); auto sign = track.sign(); histos.fill(HIST("hPt"), pt); histos.fill(HIST("hEta"), eta); histos.fill(HIST("hPhi"), phi); - int id = identifyTrack(track, cent); - bool isPi = (id == KPidPionOne); - bool isKa = (id == KPidKaonTwo); - bool isPr = (id == KPidProtonThree); - bool isSpecies[KNsp] = { - true, - isPi && sign < 0, isPi && sign > 0, isPi, - isKa && sign < 0, isKa && sign > 0, isKa, - isPr && sign < 0, isPr && sign > 0, isPr}; - - for (int isp = 0; isp < KNsp; ++isp) { - if (!isSpecies[isp]) - continue; - float eff = getEfficiency(multPV, pt, eta, static_cast(isp), 0, cfgEff); - float fake = getEfficiency(multPV, pt, eta, static_cast(isp), 1, cfgEff); - float flatW = getFlatteningWeight(vz, sign, pt, eta, phi, static_cast(isp), cfgFlat); - float w = flatW * (1.0 - fake) / eff; + float eff = getEfficiency(multPV, pt, eta, 0, cfgEff); + float fake = getEfficiency(multPV, pt, eta, 1, cfgEff); + float flatW = getFlatteningWeight(vz, sign, pt, eta, phi, cfgFlat); + float w = flatW * (1.0 - fake) / eff; + if (!std::isfinite(w) || w <= 0.f || eff <= KFloatEpsilon) { + continue; + } - if (!std::isfinite(w) || w <= 0.f || eff <= KFloatEpsilon) + for (int ieta = 0; ieta < nEta; ++ieta) { + if (eta <= etaLw[ieta] || eta > etaUp[ieta]) { continue; - - for (int ieta = 0; ieta < KNEta; ++ieta) { - if (eta <= etaLw[ieta] || eta > etaUp[ieta]) - continue; - for (int k = 0; k < KIntK; ++k) { - for (int m = 0; m < KIntM; ++m) { - sumPmwkReco[isp][ieta][m][k] += std::pow(1.0, k) * std::pow(pt, m); - sumPmwkRecoEffCor[isp][ieta][m][k] += std::pow(w, k) * std::pow(pt, m); - } - sumWkReco[isp][ieta][k] += std::pow(1.0, k); - sumWkRecoEffCor[isp][ieta][k] += std::pow(w, k); - } } - - if (isp == kInclusiveIdx) { - histos.fill(HIST("hEtaPhiReco"), vz, sign, pt, eta, phi); - histos.fill(HIST("hEtaPhiRecoWtd"), vz, sign, pt, eta, phi, w); - histos.fill(HIST("hEtaPhiRecoEffWtd"), vz, sign, pt, eta, phi, (1.0 - fake) / eff); - } else if (isp == kPiMinusIdx) { - histos.fill(HIST("hEtaPhiReco_PiMinus"), vz, sign, pt, eta, phi); - histos.fill(HIST("hEtaPhiRecoWtd_PiMinus"), vz, sign, pt, eta, phi, w); - histos.fill(HIST("hEtaPhiRecoEffWtd_PiMinus"), vz, sign, pt, eta, phi, (1.0 - fake) / eff); - } else if (isp == kPiPlusIdx) { - histos.fill(HIST("hEtaPhiReco_PiPlus"), vz, sign, pt, eta, phi); - histos.fill(HIST("hEtaPhiRecoWtd_PiPlus"), vz, sign, pt, eta, phi, w); - histos.fill(HIST("hEtaPhiRecoEffWtd_PiPlus"), vz, sign, pt, eta, phi, (1.0 - fake) / eff); - } else if (isp == kPiAllIdx) { - histos.fill(HIST("hEtaPhiReco_PiAll"), vz, sign, pt, eta, phi); - histos.fill(HIST("hEtaPhiRecoWtd_PiAll"), vz, sign, pt, eta, phi, w); - histos.fill(HIST("hEtaPhiRecoEffWtd_PiAll"), vz, sign, pt, eta, phi, (1.0 - fake) / eff); - } else if (isp == kKaMinusIdx) { - histos.fill(HIST("hEtaPhiReco_KaMinus"), vz, sign, pt, eta, phi); - histos.fill(HIST("hEtaPhiRecoWtd_KaMinus"), vz, sign, pt, eta, phi, w); - histos.fill(HIST("hEtaPhiRecoEffWtd_KaMinus"), vz, sign, pt, eta, phi, (1.0 - fake) / eff); - } else if (isp == kKaPlusIdx) { - histos.fill(HIST("hEtaPhiReco_KaPlus"), vz, sign, pt, eta, phi); - histos.fill(HIST("hEtaPhiRecoWtd_KaPlus"), vz, sign, pt, eta, phi, w); - histos.fill(HIST("hEtaPhiRecoEffWtd_KaPlus"), vz, sign, pt, eta, phi, (1.0 - fake) / eff); - } else if (isp == kKaAllIdx) { - histos.fill(HIST("hEtaPhiReco_KaAll"), vz, sign, pt, eta, phi); - histos.fill(HIST("hEtaPhiRecoWtd_KaAll"), vz, sign, pt, eta, phi, w); - histos.fill(HIST("hEtaPhiRecoEffWtd_KaAll"), vz, sign, pt, eta, phi, (1.0 - fake) / eff); - } else if (isp == kPrIdx) { - histos.fill(HIST("hEtaPhiReco_Pr"), vz, sign, pt, eta, phi); - histos.fill(HIST("hEtaPhiRecoWtd_Pr"), vz, sign, pt, eta, phi, w); - histos.fill(HIST("hEtaPhiRecoEffWtd_Pr"), vz, sign, pt, eta, phi, (1.0 - fake) / eff); - } else if (isp == kAntiPrIdx) { - histos.fill(HIST("hEtaPhiReco_AntiPr"), vz, sign, pt, eta, phi); - histos.fill(HIST("hEtaPhiRecoWtd_AntiPr"), vz, sign, pt, eta, phi, w); - histos.fill(HIST("hEtaPhiRecoEffWtd_AntiPr"), vz, sign, pt, eta, phi, (1.0 - fake) / eff); - } else if (isp == kPrAllIdx) { - histos.fill(HIST("hEtaPhiReco_PrAll"), vz, sign, pt, eta, phi); - histos.fill(HIST("hEtaPhiRecoWtd_PrAll"), vz, sign, pt, eta, phi, w); - histos.fill(HIST("hEtaPhiRecoEffWtd_PrAll"), vz, sign, pt, eta, phi, (1.0 - fake) / eff); + for (int k = 0; k < KIntK; ++k) { + for (int m = 0; m < KIntM; ++m) { + sumPmwkReco[ieta][m][k] += std::pow(1.0, k) * std::pow(pt, m); + sumPmwkRecoEffCor[ieta][m][k] += std::pow(w, k) * std::pow(pt, m); + } + sumWkReco[ieta][k] += std::pow(1.0, k); + sumWkRecoEffCor[ieta][k] += std::pow(w, k); } } - } // trkslice - if (!hasMinTracksInAllEtaBins(sumWkTru) || !hasMinTracksInAllEtaBins(sumWkReco)) + histos.fill(HIST("hEtaPhiReco"), vz, sign, pt, eta, phi); + histos.fill(HIST("hEtaPhiRecoWtd"), vz, sign, pt, eta, phi, w); + histos.fill(HIST("hEtaPhiRecoEffWtd"), vz, sign, pt, eta, phi, (1.0 - fake) / eff); + } + + if (!hasMinTracksInAllEtaBins(sumWkTru) || !hasMinTracksInAllEtaBins(sumWkReco)) { return; + } - for (int ieta = 0; ieta < KNEta; ++ieta) { - const int ibx = state.pmeanTruNchEtabinSpbinStep2->GetXaxis()->FindBin(mcCollision.multNTracksPV()); + for (int ieta = 0; ieta < nEta; ++ieta) { + const int ibx = state.pmeanTruNchEtabinStep2->GetXaxis()->FindBin(multPV); const int iby = ieta + 1; - for (int isp = 0; isp < KNsp; ++isp) { - const int ibz = isp + 1; - - meanTruMult[isp][ieta] = sumWkTru[isp][ieta][1]; - meanRecoMult[isp][ieta] = sumWkReco[isp][ieta][1]; - meanRecoEffCorMult[isp][ieta] = sumWkRecoEffCor[isp][ieta][1]; + meanTruMult[ieta] = sumWkTru[ieta][1]; + meanRecoMult[ieta] = sumWkReco[ieta][1]; + meanRecoEffCorMult[ieta] = sumWkRecoEffCor[ieta][1]; - float mmptTru = state.pmeanTruNchEtabinSpbinStep2->GetBinContent(ibx, iby, ibz); - float mmptReco = state.pmeanRecoNchEtabinSpbinStep2->GetBinContent(ibx, iby, ibz); - float mmptRecoEffCor = state.pmeanRecoEffcorrNchEtabinSpbinStep2->GetBinContent(ibx, iby, ibz); + float mmptTru = state.pmeanTruNchEtabinStep2->GetBinContent(ibx, iby); + float mmptReco = state.pmeanRecoNchEtabinStep2->GetBinContent(ibx, iby); + float mmptRecoEffCor = state.pmeanRecoEffcorrNchEtabinStep2->GetBinContent(ibx, iby); - float mmMultTru = state.pmeanMultTruNchEtabinSpbinStep2->GetBinContent(ibx, iby, ibz); - float mmMultReco = state.pmeanMultRecoNchEtabinSpbinStep2->GetBinContent(ibx, iby, ibz); - float mmMultRecoEffCor = state.pmeanMultRecoEffcorrNchEtabinSpbinStep2->GetBinContent(ibx, iby, ibz); + float mmMultTru = state.pmeanMultTruNchEtabinStep2->GetBinContent(ibx, iby); + float mmMultReco = state.pmeanMultRecoNchEtabinStep2->GetBinContent(ibx, iby); + float mmMultRecoEffCor = state.pmeanMultRecoEffcorrNchEtabinStep2->GetBinContent(ibx, iby); - if (std::isfinite(mmptTru)) - std::tie(meanTru[isp][ieta], c2Tru[isp][ieta]) = calculateMeanAndC2FromSums(sumPmwkTru[isp][ieta], sumWkTru[isp][ieta], mmptTru); - if (std::isfinite(mmptReco)) - std::tie(meanReco[isp][ieta], c2Reco[isp][ieta]) = calculateMeanAndC2FromSums(sumPmwkReco[isp][ieta], sumWkReco[isp][ieta], mmptReco); - if (std::isfinite(mmptRecoEffCor)) - std::tie(meanRecoEffCor[isp][ieta], c2RecoEffCor[isp][ieta]) = calculateMeanAndC2FromSums(sumPmwkRecoEffCor[isp][ieta], sumWkRecoEffCor[isp][ieta], mmptRecoEffCor); + if (std::isfinite(mmptTru)) { + std::tie(meanTru[ieta], c2Tru[ieta]) = calculateMeanAndC2FromSums(sumPmwkTru[ieta], sumWkTru[ieta], mmptTru); + } + if (std::isfinite(mmptReco)) { + std::tie(meanReco[ieta], c2Reco[ieta]) = calculateMeanAndC2FromSums(sumPmwkReco[ieta], sumWkReco[ieta], mmptReco); + } + if (std::isfinite(mmptRecoEffCor)) { + std::tie(meanRecoEffCor[ieta], c2RecoEffCor[ieta]) = calculateMeanAndC2FromSums(sumPmwkRecoEffCor[ieta], sumWkRecoEffCor[ieta], mmptRecoEffCor); + } - if (mmptTru != 0.0f) - p1kBarTru[isp][ieta] = meanTru[isp][ieta] - mmptTru; - if (mmptReco != 0.0f) - p1kBarReco[isp][ieta] = meanReco[isp][ieta] - mmptReco; - if (mmptRecoEffCor != 0.0f) - p1kBarRecoEffCor[isp][ieta] = meanRecoEffCor[isp][ieta] - mmptRecoEffCor; + if (mmptTru != 0.0f) { + p1kBarTru[ieta] = meanTru[ieta] - mmptTru; + } + if (mmptReco != 0.0f) { + p1kBarReco[ieta] = meanReco[ieta] - mmptReco; + } + if (mmptRecoEffCor != 0.0f) { + p1kBarRecoEffCor[ieta] = meanRecoEffCor[ieta] - mmptRecoEffCor; + } - if (mmMultTru != 0.0f) - p1kBarTruMult[isp][ieta] = meanTruMult[isp][ieta] - mmMultTru; - if (mmMultReco != 0.0f) - p1kBarRecoMult[isp][ieta] = meanRecoMult[isp][ieta] - mmMultReco; - if (mmMultRecoEffCor != 0.0f) - p1kBarRecoEffCorMult[isp][ieta] = meanRecoEffCorMult[isp][ieta] - mmMultRecoEffCor; + if (mmMultTru != 0.0f) { + p1kBarTruMult[ieta] = meanTruMult[ieta] - mmMultTru; + } + if (mmMultReco != 0.0f) { + p1kBarRecoMult[ieta] = meanRecoMult[ieta] - mmMultReco; + } + if (mmMultRecoEffCor != 0.0f) { + p1kBarRecoEffCorMult[ieta] = meanRecoEffCorMult[ieta] - mmMultRecoEffCor; } } @@ -2427,586 +1632,229 @@ struct RadialFlowDecorr { if (mcCollision.has_foundFT0()) { const auto& ft0 = mcCollision.foundFT0(); for (std::size_t iCh = 0; iCh < ft0.channelA().size(); iCh++) { - float ampl = ft0.amplitudeA()[iCh]; - amplFT0A += ampl; + amplFT0A += ft0.amplitudeA()[iCh]; } for (std::size_t iCh = 0; iCh < ft0.channelC().size(); iCh++) { - float ampl = ft0.amplitudeC()[iCh]; - amplFT0C += ampl; + amplFT0C += ft0.amplitudeC()[iCh]; } } + p1kBarFt0A = amplFT0A - state.pmeanFT0AmultpvStep2->GetBinContent(state.pmeanFT0AmultpvStep2->GetXaxis()->FindBin(multPV)); + p1kBarFt0C = amplFT0C - state.pmeanFT0CmultpvStep2->GetBinContent(state.pmeanFT0CmultpvStep2->GetXaxis()->FindBin(multPV)); - for (int isp = 0; isp < KNsp; ++isp) { - for (int ieta = 0; ieta < KNEta; ++ieta) { - histos.fill(HIST("MCGen/Prof_Cent_NEta_Nsp_Nchrec"), cent, ieta, isp, sumWkTru[isp][ieta][1]); - histos.fill(HIST("MCGen/Prof_Mult_NEta_Nsp_Nchrec"), multPV, ieta, isp, sumWkTru[isp][ieta][1]); + // per-eta counts & means + for (int ieta = 0; ieta < nEta; ++ieta) { + histos.fill(HIST("MCGen/Prof_Cent_NEta_Nchrec"), cent, ieta, sumWkTru[ieta][1]); + histos.fill(HIST("MCGen/Prof_Mult_NEta_Nchrec"), multPV, ieta, sumWkTru[ieta][1]); + histos.fill(HIST("MCReco/Prof_Cent_NEta_Nchrec"), cent, ieta, sumWkReco[ieta][1]); + histos.fill(HIST("MCReco/Prof_Mult_NEta_Nchrec"), multPV, ieta, sumWkReco[ieta][1]); + histos.fill(HIST("MCRecoEffCorr/Prof_Cent_NEta_Nchrec"), cent, ieta, sumWkRecoEffCor[ieta][1]); + histos.fill(HIST("MCRecoEffCorr/Prof_Mult_NEta_Nchrec"), multPV, ieta, sumWkRecoEffCor[ieta][1]); + + if (sumWkTru[ieta][1] > 1.0f) { + histos.fill(HIST("MCGen/Prof_Cent_NEta_MeanpT"), cent, ieta, meanTru[ieta]); + histos.fill(HIST("MCGen/Prof_Mult_NEta_MeanpT"), multPV, ieta, meanTru[ieta]); + } + if (sumWkReco[ieta][1] > 1.0f) { + histos.fill(HIST("MCReco/Prof_Cent_NEta_MeanpT"), cent, ieta, meanReco[ieta]); + histos.fill(HIST("MCReco/Prof_Mult_NEta_MeanpT"), multPV, ieta, meanReco[ieta]); + } + if (sumWkRecoEffCor[ieta][1] > 1.0f) { + histos.fill(HIST("MCRecoEffCorr/Prof_Cent_NEta_MeanpT"), cent, ieta, meanRecoEffCor[ieta]); + histos.fill(HIST("MCRecoEffCorr/Prof_Mult_NEta_MeanpT"), multPV, ieta, meanRecoEffCor[ieta]); + } + } - histos.fill(HIST("MCReco/Prof_Cent_NEta_Nsp_Nchrec"), cent, ieta, isp, sumWkReco[isp][ieta][1]); - histos.fill(HIST("MCReco/Prof_Mult_NEta_Nsp_Nchrec"), multPV, ieta, isp, sumWkReco[isp][ieta][1]); + // meanpT & C2 vs eta bin + for (int ieta = 0; ieta < nEta; ++ieta) { + if (std::isfinite(meanTru[ieta])) { + histos.fill(HIST("MCGen/Prof_MeanpT_Cent_etabin"), cent, ieta, meanTru[ieta]); + histos.fill(HIST("MCGen/Prof_MeanpT_Mult_etabin"), multPV, ieta, meanTru[ieta]); + } + if (std::isfinite(c2Tru[ieta])) { + histos.fill(HIST("MCGen/Prof_C2_Cent_etabin"), cent, ieta, c2Tru[ieta]); + histos.fill(HIST("MCGen/Prof_C2_Mult_etabin"), multPV, ieta, c2Tru[ieta]); + } + if (std::isfinite(meanReco[ieta])) { + histos.fill(HIST("MCReco/Prof_MeanpT_Cent_etabin"), cent, ieta, meanReco[ieta]); + histos.fill(HIST("MCReco/Prof_MeanpT_Mult_etabin"), multPV, ieta, meanReco[ieta]); + } + if (std::isfinite(c2Reco[ieta])) { + histos.fill(HIST("MCReco/Prof_C2_Cent_etabin"), cent, ieta, c2Reco[ieta]); + histos.fill(HIST("MCReco/Prof_C2_Mult_etabin"), multPV, ieta, c2Reco[ieta]); + } + if (std::isfinite(meanRecoEffCor[ieta])) { + histos.fill(HIST("MCRecoEffCorr/Prof_MeanpT_Cent_etabin"), cent, ieta, meanRecoEffCor[ieta]); + histos.fill(HIST("MCRecoEffCorr/Prof_MeanpT_Mult_etabin"), multPV, ieta, meanRecoEffCor[ieta]); + } + if (std::isfinite(c2RecoEffCor[ieta])) { + histos.fill(HIST("MCRecoEffCorr/Prof_C2_Cent_etabin"), cent, ieta, c2RecoEffCor[ieta]); + histos.fill(HIST("MCRecoEffCorr/Prof_C2_Mult_etabin"), multPV, ieta, c2RecoEffCor[ieta]); + } + } - histos.fill(HIST("MCRecoEffCorr/Prof_Cent_NEta_Nsp_Nchrec"), cent, ieta, isp, sumWkRecoEffCor[isp][ieta][1]); - histos.fill(HIST("MCRecoEffCorr/Prof_Mult_NEta_Nsp_Nchrec"), multPV, ieta, isp, sumWkRecoEffCor[isp][ieta][1]); + // mirror-pair subevent (C2Sub) & covariances vs eta bin + for (int ietaA = 1; ietaA <= (nEta - 1) / 2; ++ietaA) { + int ietaC = nEta - ietaA; - if (sumWkTru[isp][ieta][1] > 1.0f) { - histos.fill(HIST("MCGen/Prof_Cent_NEta_Nsp_MeanpT"), cent, ieta, isp, meanTru[isp][ieta]); - histos.fill(HIST("MCGen/Prof_Mult_NEta_Nsp_MeanpT"), multPV, ieta, isp, meanTru[isp][ieta]); - } - if (sumWkReco[isp][ieta][1] > 1.0f) { - histos.fill(HIST("MCReco/Prof_Cent_NEta_Nsp_MeanpT"), cent, ieta, isp, meanReco[isp][ieta]); - histos.fill(HIST("MCReco/Prof_Mult_NEta_Nsp_MeanpT"), multPV, ieta, isp, meanReco[isp][ieta]); - } - if (sumWkRecoEffCor[isp][ieta][1] > 1.0f) { - histos.fill(HIST("MCRecoEffCorr/Prof_Cent_NEta_Nsp_MeanpT"), cent, ieta, isp, meanRecoEffCor[isp][ieta]); - histos.fill(HIST("MCRecoEffCorr/Prof_Mult_NEta_Nsp_MeanpT"), multPV, ieta, isp, meanRecoEffCor[isp][ieta]); - } + float c2SubTru = p1kBarTru[ietaA] * p1kBarTru[ietaC]; + float c2SubReco = p1kBarReco[ietaA] * p1kBarReco[ietaC]; + float c2SubRecoEffCor = p1kBarRecoEffCor[ietaA] * p1kBarRecoEffCor[ietaC]; + + float covTru = p1kBarTruMult[ietaA] * p1kBarTru[ietaC]; + float covReco = p1kBarRecoMult[ietaA] * p1kBarReco[ietaC]; + float covRecoEffCor = p1kBarRecoEffCorMult[ietaA] * p1kBarRecoEffCor[ietaC]; + + if (std::isfinite(c2SubTru)) { + histos.fill(HIST("MCGen/Prof_C2Sub_Cent_etabin"), cent, ietaA, c2SubTru); + histos.fill(HIST("MCGen/Prof_C2Sub_Mult_etabin"), multPV, ietaA, c2SubTru); + } + if (std::isfinite(c2SubReco)) { + histos.fill(HIST("MCReco/Prof_C2Sub_Cent_etabin"), cent, ietaA, c2SubReco); + histos.fill(HIST("MCReco/Prof_C2Sub_Mult_etabin"), multPV, ietaA, c2SubReco); + } + if (std::isfinite(c2SubRecoEffCor)) { + histos.fill(HIST("MCRecoEffCorr/Prof_C2Sub_Cent_etabin"), cent, ietaA, c2SubRecoEffCor); + histos.fill(HIST("MCRecoEffCorr/Prof_C2Sub_Mult_etabin"), multPV, ietaA, c2SubRecoEffCor); + } + if (std::isfinite(covTru)) { + histos.fill(HIST("MCGen/Prof_Cov_Cent_etabin"), cent, ietaA, covTru); + histos.fill(HIST("MCGen/Prof_Cov_Mult_etabin"), multPV, ietaA, covTru); + } + if (std::isfinite(covReco)) { + histos.fill(HIST("MCReco/Prof_Cov_Cent_etabin"), cent, ietaA, covReco); + histos.fill(HIST("MCReco/Prof_Cov_Mult_etabin"), multPV, ietaA, covReco); + } + if (std::isfinite(covRecoEffCor)) { + histos.fill(HIST("MCRecoEffCorr/Prof_Cov_Cent_etabin"), cent, ietaA, covRecoEffCor); + histos.fill(HIST("MCRecoEffCorr/Prof_Cov_Mult_etabin"), multPV, ietaA, covRecoEffCor); } } - for (int ieta = 0; ieta < KNEta; ++ieta) { - for (int isp = 0; isp < KNsp; ++isp) { - if (std::isfinite(meanTru[isp][ieta])) { - histos.fill(HIST("MCGen/Prof_MeanpT_Cent_etabin_spbin"), cent, ieta, isp, meanTru[isp][ieta]); - histos.fill(HIST("MCGen/Prof_MeanpT_Mult_etabin_spbin"), multPV, ieta, isp, meanTru[isp][ieta]); - } - if (std::isfinite(c2Tru[isp][ieta])) { - histos.fill(HIST("MCGen/Prof_C2_Cent_etabin_spbin"), cent, ieta, isp, c2Tru[isp][ieta]); - histos.fill(HIST("MCGen/Prof_C2_Mult_etabin_spbin"), multPV, ieta, isp, c2Tru[isp][ieta]); - } - if (std::isfinite(meanReco[isp][ieta])) { - histos.fill(HIST("MCReco/Prof_MeanpT_Cent_etabin_spbin"), cent, ieta, isp, meanReco[isp][ieta]); - histos.fill(HIST("MCReco/Prof_MeanpT_Mult_etabin_spbin"), multPV, ieta, isp, meanReco[isp][ieta]); - } - if (std::isfinite(c2Reco[isp][ieta])) { - histos.fill(HIST("MCReco/Prof_C2_Cent_etabin_spbin"), cent, ieta, isp, c2Reco[isp][ieta]); - histos.fill(HIST("MCReco/Prof_C2_Mult_etabin_spbin"), multPV, ieta, isp, c2Reco[isp][ieta]); - } - if (std::isfinite(meanRecoEffCor[isp][ieta])) { - histos.fill(HIST("MCRecoEffCorr/Prof_MeanpT_Cent_etabin_spbin"), cent, ieta, isp, meanRecoEffCor[isp][ieta]); - histos.fill(HIST("MCRecoEffCorr/Prof_MeanpT_Mult_etabin_spbin"), multPV, ieta, isp, meanRecoEffCor[isp][ieta]); - } - if (std::isfinite(c2RecoEffCor[isp][ieta])) { - histos.fill(HIST("MCRecoEffCorr/Prof_C2_Cent_etabin_spbin"), cent, ieta, isp, c2RecoEffCor[isp][ieta]); - histos.fill(HIST("MCRecoEffCorr/Prof_C2_Mult_etabin_spbin"), multPV, ieta, isp, c2RecoEffCor[isp][ieta]); - } + // FT0 covariance vs narrow eta bin (full range, indexed by the actual pT bin) + for (int ieta = 1; ieta < nEta; ++ieta) { + float covFT0ATru = p1kBarFt0A * p1kBarTru[ieta]; + float covFT0AReco = p1kBarFt0A * p1kBarReco[ieta]; + float covFT0ARecoEffCor = p1kBarFt0A * p1kBarRecoEffCor[ieta]; + float covFT0CTru = p1kBarFt0C * p1kBarTru[ieta]; + float covFT0CReco = p1kBarFt0C * p1kBarReco[ieta]; + float covFT0CRecoEffCor = p1kBarFt0C * p1kBarRecoEffCor[ieta]; + + if (std::isfinite(covFT0ATru)) { + histos.fill(HIST("MCGen/Prof_CovFT0A_Cent_etabin"), cent, ieta, covFT0ATru); + histos.fill(HIST("MCGen/Prof_CovFT0A_Mult_etabin"), multPV, ieta, covFT0ATru); + } + if (std::isfinite(covFT0AReco)) { + histos.fill(HIST("MCReco/Prof_CovFT0A_Cent_etabin"), cent, ieta, covFT0AReco); + histos.fill(HIST("MCReco/Prof_CovFT0A_Mult_etabin"), multPV, ieta, covFT0AReco); + } + if (std::isfinite(covFT0ARecoEffCor)) { + histos.fill(HIST("MCRecoEffCorr/Prof_CovFT0A_Cent_etabin"), cent, ieta, covFT0ARecoEffCor); + histos.fill(HIST("MCRecoEffCorr/Prof_CovFT0A_Mult_etabin"), multPV, ieta, covFT0ARecoEffCor); + } + if (std::isfinite(covFT0CTru)) { + histos.fill(HIST("MCGen/Prof_CovFT0C_Cent_etabin"), cent, ieta, covFT0CTru); + histos.fill(HIST("MCGen/Prof_CovFT0C_Mult_etabin"), multPV, ieta, covFT0CTru); + } + if (std::isfinite(covFT0CReco)) { + histos.fill(HIST("MCReco/Prof_CovFT0C_Cent_etabin"), cent, ieta, covFT0CReco); + histos.fill(HIST("MCReco/Prof_CovFT0C_Mult_etabin"), multPV, ieta, covFT0CReco); + } + if (std::isfinite(covFT0CRecoEffCor)) { + histos.fill(HIST("MCRecoEffCorr/Prof_CovFT0C_Cent_etabin"), cent, ieta, covFT0CRecoEffCor); + histos.fill(HIST("MCRecoEffCorr/Prof_CovFT0C_Mult_etabin"), multPV, ieta, covFT0CRecoEffCor); } } - p1kBarFt0A = amplFT0A - state.pmeanFT0AmultpvStep2->GetBinContent(state.pmeanFT0AmultpvStep2->GetXaxis()->FindBin(multPV)); - p1kBarFt0C = amplFT0C - state.pmeanFT0CmultpvStep2->GetBinContent(state.pmeanFT0CmultpvStep2->GetXaxis()->FindBin(multPV)); + // full 2D subevent map + for (int ietaA = 1; ietaA < nEta; ++ietaA) { + for (int ietaC = 1; ietaC < nEta; ++ietaC) { + float etaValA = (etaLw[ietaA] + etaUp[ietaA]) / 2.0f; + float etaValB = (etaLw[ietaC] + etaUp[ietaC]) / 2.0f; + float gap = etaValA - etaValB; + float sum = (etaValA + etaValB); - for (int ietaA = 1; ietaA <= (KNEta - 1) / 2; ++ietaA) { - int ietaC = KNEta - ietaA; - for (int isp = 0; isp < KNsp; ++isp) { - float c2SubTru = p1kBarTru[isp][ietaA] * p1kBarTru[isp][ietaC]; - float c2SubReco = p1kBarReco[isp][ietaA] * p1kBarReco[isp][ietaC]; - float c2SubRecoEffCor = p1kBarRecoEffCor[isp][ietaA] * p1kBarRecoEffCor[isp][ietaC]; + float c2SubTru = (ietaA == ietaC) ? static_cast(c2Tru[ietaA]) : p1kBarTru[ietaA] * p1kBarTru[ietaC]; + float c2SubReco = (ietaA == ietaC) ? static_cast(c2Reco[ietaA]) : p1kBarReco[ietaA] * p1kBarReco[ietaC]; + float c2SubRecoEffCor = (ietaA == ietaC) ? static_cast(c2RecoEffCor[ietaA]) : p1kBarRecoEffCor[ietaA] * p1kBarRecoEffCor[ietaC]; - float covTru = p1kBarTruMult[isp][ietaA] * p1kBarTru[isp][ietaC]; - float covReco = p1kBarRecoMult[isp][ietaA] * p1kBarReco[isp][ietaC]; - float covRecoEffCor = p1kBarRecoEffCorMult[isp][ietaA] * p1kBarRecoEffCor[isp][ietaC]; + float covTru = p1kBarTruMult[ietaA] * p1kBarTru[ietaC]; + float covReco = p1kBarRecoMult[ietaA] * p1kBarReco[ietaC]; + float covRecoEffCor = p1kBarRecoEffCorMult[ietaA] * p1kBarRecoEffCor[ietaC]; - float covFT0ATru = p1kBarFt0A * p1kBarTru[isp][ietaC]; - float covFT0AReco = p1kBarFt0A * p1kBarReco[isp][ietaC]; - float covFT0ARecoEffCor = p1kBarFt0A * p1kBarRecoEffCor[isp][ietaC]; + float covFT0ATru = p1kBarFt0A * p1kBarTru[ietaC]; + float covFT0AReco = p1kBarFt0A * p1kBarReco[ietaC]; + float covFT0ARecoEffCor = p1kBarFt0A * p1kBarRecoEffCor[ietaC]; - float covFT0CTru = p1kBarFt0C * p1kBarTru[isp][ietaA]; - float covFT0CReco = p1kBarFt0C * p1kBarReco[isp][ietaA]; - float covFT0CRecoEffCor = p1kBarFt0C * p1kBarRecoEffCor[isp][ietaA]; + float covFT0CTru = p1kBarFt0C * p1kBarTru[ietaA]; + float covFT0CReco = p1kBarFt0C * p1kBarReco[ietaA]; + float covFT0CRecoEffCor = p1kBarFt0C * p1kBarRecoEffCor[ietaA]; if (std::isfinite(c2SubTru)) { - histos.fill(HIST("MCGen/Prof_C2Sub_Cent_etabin_spbin"), cent, ietaA, isp, c2SubTru); - histos.fill(HIST("MCGen/Prof_C2Sub_Mult_etabin_spbin"), multPV, ietaA, isp, c2SubTru); + histos.fill(HIST("MCGen/Prof_C2Sub2D_Cent_etaA_etaC"), cent, etaValA, etaValB, c2SubTru); + histos.fill(HIST("MCGen/Prof_GapSum2D"), cent, gap, sum, c2SubTru); } if (std::isfinite(c2SubReco)) { - histos.fill(HIST("MCReco/Prof_C2Sub_Cent_etabin_spbin"), cent, ietaA, isp, c2SubReco); - histos.fill(HIST("MCReco/Prof_C2Sub_Mult_etabin_spbin"), multPV, ietaA, isp, c2SubReco); + histos.fill(HIST("MCReco/Prof_C2Sub2D_Cent_etaA_etaC"), cent, etaValA, etaValB, c2SubReco); + histos.fill(HIST("MCReco/Prof_GapSum2D"), cent, gap, sum, c2SubReco); } if (std::isfinite(c2SubRecoEffCor)) { - histos.fill(HIST("MCRecoEffCorr/Prof_C2Sub_Cent_etabin_spbin"), cent, ietaA, isp, c2SubRecoEffCor); - histos.fill(HIST("MCRecoEffCorr/Prof_C2Sub_Mult_etabin_spbin"), multPV, ietaA, isp, c2SubRecoEffCor); + histos.fill(HIST("MCRecoEffCorr/Prof_C2Sub2D_Cent_etaA_etaC"), cent, etaValA, etaValB, c2SubRecoEffCor); + histos.fill(HIST("MCRecoEffCorr/Prof_GapSum2D"), cent, gap, sum, c2SubRecoEffCor); } + if (std::isfinite(covTru)) { - histos.fill(HIST("MCGen/Prof_Cov_Cent_etabin_spbin"), cent, ietaA, isp, covTru); - histos.fill(HIST("MCGen/Prof_Cov_Mult_etabin_spbin"), multPV, ietaA, isp, covTru); + histos.fill(HIST("MCGen/Prof_Cov2D_Cent_etaA_etaC"), cent, etaValA, etaValB, covTru); } if (std::isfinite(covReco)) { - histos.fill(HIST("MCReco/Prof_Cov_Cent_etabin_spbin"), cent, ietaA, isp, covReco); - histos.fill(HIST("MCReco/Prof_Cov_Mult_etabin_spbin"), multPV, ietaA, isp, covReco); + histos.fill(HIST("MCReco/Prof_Cov2D_Cent_etaA_etaC"), cent, etaValA, etaValB, covReco); } if (std::isfinite(covRecoEffCor)) { - histos.fill(HIST("MCRecoEffCorr/Prof_Cov_Cent_etabin_spbin"), cent, ietaA, isp, covRecoEffCor); - histos.fill(HIST("MCRecoEffCorr/Prof_Cov_Mult_etabin_spbin"), multPV, ietaA, isp, covRecoEffCor); + histos.fill(HIST("MCRecoEffCorr/Prof_Cov2D_Cent_etaA_etaC"), cent, etaValA, etaValB, covRecoEffCor); } if (std::isfinite(covFT0ATru)) { - histos.fill(HIST("MCGen/Prof_CovFT0A_Cent_etabin_spbin"), cent, ietaA, isp, covFT0ATru); - histos.fill(HIST("MCGen/Prof_CovFT0A_Mult_etabin_spbin"), multPV, ietaA, isp, covFT0ATru); + histos.fill(HIST("MCGen/Prof_CovFT0A2D_Cent_etaA_etaC"), cent, etaValA, etaValB, covFT0ATru); } if (std::isfinite(covFT0AReco)) { - histos.fill(HIST("MCReco/Prof_CovFT0A_Cent_etabin_spbin"), cent, ietaA, isp, covFT0AReco); - histos.fill(HIST("MCReco/Prof_CovFT0A_Mult_etabin_spbin"), multPV, ietaA, isp, covFT0AReco); + histos.fill(HIST("MCReco/Prof_CovFT0A2D_Cent_etaA_etaC"), cent, etaValA, etaValB, covFT0AReco); } if (std::isfinite(covFT0ARecoEffCor)) { - histos.fill(HIST("MCRecoEffCorr/Prof_CovFT0A_Cent_etabin_spbin"), cent, ietaA, isp, covFT0ARecoEffCor); - histos.fill(HIST("MCRecoEffCorr/Prof_CovFT0A_Mult_etabin_spbin"), multPV, ietaA, isp, covFT0ARecoEffCor); + histos.fill(HIST("MCRecoEffCorr/Prof_CovFT0A2D_Cent_etaA_etaC"), cent, etaValA, etaValB, covFT0ARecoEffCor); } if (std::isfinite(covFT0CTru)) { - histos.fill(HIST("MCGen/Prof_CovFT0C_Cent_etabin_spbin"), cent, ietaA, isp, covFT0CTru); - histos.fill(HIST("MCGen/Prof_CovFT0C_Mult_etabin_spbin"), multPV, ietaA, isp, covFT0CTru); + histos.fill(HIST("MCGen/Prof_CovFT0C2D_Cent_etaA_etaC"), cent, etaValA, etaValB, covFT0CTru); } if (std::isfinite(covFT0CReco)) { - histos.fill(HIST("MCReco/Prof_CovFT0C_Cent_etabin_spbin"), cent, ietaA, isp, covFT0CReco); - histos.fill(HIST("MCReco/Prof_CovFT0C_Mult_etabin_spbin"), multPV, ietaA, isp, covFT0CReco); + histos.fill(HIST("MCReco/Prof_CovFT0C2D_Cent_etaA_etaC"), cent, etaValA, etaValB, covFT0CReco); } if (std::isfinite(covFT0CRecoEffCor)) { - histos.fill(HIST("MCRecoEffCorr/Prof_CovFT0C_Cent_etabin_spbin"), cent, ietaA, isp, covFT0CRecoEffCor); - histos.fill(HIST("MCRecoEffCorr/Prof_CovFT0C_Mult_etabin_spbin"), multPV, ietaA, isp, covFT0CRecoEffCor); - } - } - } - - for (int ietaA = 1; ietaA < KNEta; ++ietaA) { - for (int ietaC = 1; ietaC < KNEta; ++ietaC) { - - float etaValA = (etaLw[ietaA] + etaUp[ietaA]) / 2.0f; - float etaValB = (etaLw[ietaC] + etaUp[ietaC]) / 2.0f; - float gap = etaValA - etaValB; - float sum = (etaValA + etaValB); - for (int isp = 0; isp < KNsp; ++isp) { - - float c2SubTru = (ietaA == ietaC) ? static_cast(c2Tru[isp][ietaA]) - : p1kBarTru[isp][ietaA] * p1kBarTru[isp][ietaC]; - float c2SubReco = (ietaA == ietaC) ? static_cast(c2Reco[isp][ietaA]) - : p1kBarReco[isp][ietaA] * p1kBarReco[isp][ietaC]; - float c2SubRecoEffCor = (ietaA == ietaC) ? static_cast(c2RecoEffCor[isp][ietaA]) - : p1kBarRecoEffCor[isp][ietaA] * p1kBarRecoEffCor[isp][ietaC]; - - float covTru = p1kBarTruMult[isp][ietaA] * p1kBarTru[isp][ietaC]; - float covReco = p1kBarRecoMult[isp][ietaA] * p1kBarReco[isp][ietaC]; - float covRecoEffCor = p1kBarRecoEffCorMult[isp][ietaA] * p1kBarRecoEffCor[isp][ietaC]; - - float covFT0ATru = p1kBarFt0A * p1kBarTru[isp][ietaC]; - float covFT0AReco = p1kBarFt0A * p1kBarReco[isp][ietaC]; - float covFT0ARecoEffCor = p1kBarFt0A * p1kBarRecoEffCor[isp][ietaC]; - - float covFT0CTru = p1kBarFt0C * p1kBarTru[isp][ietaA]; - float covFT0CReco = p1kBarFt0C * p1kBarReco[isp][ietaA]; - float covFT0CRecoEffCor = p1kBarFt0C * p1kBarRecoEffCor[isp][ietaA]; - - if (isp == kInclusiveIdx) { - if (std::isfinite(c2SubTru)) { - histos.fill(HIST("MCGen/Prof_C2Sub2D_Cent_etaA_etaC"), cent, etaValA, etaValB, c2SubTru); - histos.fill(HIST("MCGen/Prof_GapSum2D"), cent, gap, sum, c2SubTru); - } - if (std::isfinite(c2SubReco)) { - histos.fill(HIST("MCReco/Prof_C2Sub2D_Cent_etaA_etaC"), cent, etaValA, etaValB, c2SubReco); - histos.fill(HIST("MCReco/Prof_GapSum2D"), cent, gap, sum, c2SubReco); - } - if (std::isfinite(c2SubRecoEffCor)) { - histos.fill(HIST("MCRecoEffCorr/Prof_C2Sub2D_Cent_etaA_etaC"), cent, etaValA, etaValB, c2SubRecoEffCor); - histos.fill(HIST("MCRecoEffCorr/Prof_GapSum2D"), cent, gap, sum, c2SubRecoEffCor); - } - - if (std::isfinite(covTru)) - histos.fill(HIST("MCGen/Prof_Cov2D_Cent_etaA_etaC"), cent, etaValA, etaValB, covTru); - if (std::isfinite(covReco)) - histos.fill(HIST("MCReco/Prof_Cov2D_Cent_etaA_etaC"), cent, etaValA, etaValB, covReco); - if (std::isfinite(covRecoEffCor)) - histos.fill(HIST("MCRecoEffCorr/Prof_Cov2D_Cent_etaA_etaC"), cent, etaValA, etaValB, covRecoEffCor); - - if (std::isfinite(covFT0ATru)) - histos.fill(HIST("MCGen/Prof_CovFT0A2D_Cent_etaA_etaC"), cent, etaValA, etaValB, covFT0ATru); - if (std::isfinite(covFT0AReco)) - histos.fill(HIST("MCReco/Prof_CovFT0A2D_Cent_etaA_etaC"), cent, etaValA, etaValB, covFT0AReco); - if (std::isfinite(covFT0ARecoEffCor)) - histos.fill(HIST("MCRecoEffCorr/Prof_CovFT0A2D_Cent_etaA_etaC"), cent, etaValA, etaValB, covFT0ARecoEffCor); - - if (std::isfinite(covFT0CTru)) - histos.fill(HIST("MCGen/Prof_CovFT0C2D_Cent_etaA_etaC"), cent, etaValA, etaValB, covFT0CTru); - if (std::isfinite(covFT0CReco)) - histos.fill(HIST("MCReco/Prof_CovFT0C2D_Cent_etaA_etaC"), cent, etaValA, etaValB, covFT0CReco); - if (std::isfinite(covFT0CRecoEffCor)) - histos.fill(HIST("MCRecoEffCorr/Prof_CovFT0C2D_Cent_etaA_etaC"), cent, etaValA, etaValB, covFT0CRecoEffCor); - - } else if (isp == kPiMinusIdx) { - if (std::isfinite(c2SubTru)) { - histos.fill(HIST("MCGen/Prof_C2Sub2D_Cent_etaA_etaC_PiMinus"), cent, etaValA, etaValB, c2SubTru); - histos.fill(HIST("MCGen/Prof_GapSum2D_PiMinus"), cent, gap, sum, c2SubTru); - } - if (std::isfinite(c2SubReco)) { - histos.fill(HIST("MCReco/Prof_C2Sub2D_Cent_etaA_etaC_PiMinus"), cent, etaValA, etaValB, c2SubReco); - histos.fill(HIST("MCReco/Prof_GapSum2D_PiMinus"), cent, gap, sum, c2SubReco); - } - if (std::isfinite(c2SubRecoEffCor)) { - histos.fill(HIST("MCRecoEffCorr/Prof_C2Sub2D_Cent_etaA_etaC_PiMinus"), cent, etaValA, etaValB, c2SubRecoEffCor); - histos.fill(HIST("MCRecoEffCorr/Prof_GapSum2D_PiMinus"), cent, gap, sum, c2SubRecoEffCor); - } - - if (std::isfinite(covTru)) - histos.fill(HIST("MCGen/Prof_Cov2D_Cent_etaA_etaC_PiMinus"), cent, etaValA, etaValB, covTru); - if (std::isfinite(covReco)) - histos.fill(HIST("MCReco/Prof_Cov2D_Cent_etaA_etaC_PiMinus"), cent, etaValA, etaValB, covReco); - if (std::isfinite(covRecoEffCor)) - histos.fill(HIST("MCRecoEffCorr/Prof_Cov2D_Cent_etaA_etaC_PiMinus"), cent, etaValA, etaValB, covRecoEffCor); - - if (std::isfinite(covFT0ATru)) - histos.fill(HIST("MCGen/Prof_CovFT0A2D_Cent_etaA_etaC_PiMinus"), cent, etaValA, etaValB, covFT0ATru); - if (std::isfinite(covFT0AReco)) - histos.fill(HIST("MCReco/Prof_CovFT0A2D_Cent_etaA_etaC_PiMinus"), cent, etaValA, etaValB, covFT0AReco); - if (std::isfinite(covFT0ARecoEffCor)) - histos.fill(HIST("MCRecoEffCorr/Prof_CovFT0A2D_Cent_etaA_etaC_PiMinus"), cent, etaValA, etaValB, covFT0ARecoEffCor); - - if (std::isfinite(covFT0CTru)) - histos.fill(HIST("MCGen/Prof_CovFT0C2D_Cent_etaA_etaC_PiMinus"), cent, etaValA, etaValB, covFT0CTru); - if (std::isfinite(covFT0CReco)) - histos.fill(HIST("MCReco/Prof_CovFT0C2D_Cent_etaA_etaC_PiMinus"), cent, etaValA, etaValB, covFT0CReco); - if (std::isfinite(covFT0CRecoEffCor)) - histos.fill(HIST("MCRecoEffCorr/Prof_CovFT0C2D_Cent_etaA_etaC_PiMinus"), cent, etaValA, etaValB, covFT0CRecoEffCor); - - } else if (isp == kPiPlusIdx) { - if (std::isfinite(c2SubTru)) { - histos.fill(HIST("MCGen/Prof_C2Sub2D_Cent_etaA_etaC_PiPlus"), cent, etaValA, etaValB, c2SubTru); - histos.fill(HIST("MCGen/Prof_GapSum2D_PiPlus"), cent, gap, sum, c2SubTru); - } - if (std::isfinite(c2SubReco)) { - histos.fill(HIST("MCReco/Prof_C2Sub2D_Cent_etaA_etaC_PiPlus"), cent, etaValA, etaValB, c2SubReco); - histos.fill(HIST("MCReco/Prof_GapSum2D_PiPlus"), cent, gap, sum, c2SubReco); - } - if (std::isfinite(c2SubRecoEffCor)) { - histos.fill(HIST("MCRecoEffCorr/Prof_C2Sub2D_Cent_etaA_etaC_PiPlus"), cent, etaValA, etaValB, c2SubRecoEffCor); - histos.fill(HIST("MCRecoEffCorr/Prof_GapSum2D_PiPlus"), cent, gap, sum, c2SubRecoEffCor); - } - - if (std::isfinite(covTru)) - histos.fill(HIST("MCGen/Prof_Cov2D_Cent_etaA_etaC_PiPlus"), cent, etaValA, etaValB, covTru); - if (std::isfinite(covReco)) - histos.fill(HIST("MCReco/Prof_Cov2D_Cent_etaA_etaC_PiPlus"), cent, etaValA, etaValB, covReco); - if (std::isfinite(covRecoEffCor)) - histos.fill(HIST("MCRecoEffCorr/Prof_Cov2D_Cent_etaA_etaC_PiPlus"), cent, etaValA, etaValB, covRecoEffCor); - - if (std::isfinite(covFT0ATru)) - histos.fill(HIST("MCGen/Prof_CovFT0A2D_Cent_etaA_etaC_PiPlus"), cent, etaValA, etaValB, covFT0ATru); - if (std::isfinite(covFT0AReco)) - histos.fill(HIST("MCReco/Prof_CovFT0A2D_Cent_etaA_etaC_PiPlus"), cent, etaValA, etaValB, covFT0AReco); - if (std::isfinite(covFT0ARecoEffCor)) - histos.fill(HIST("MCRecoEffCorr/Prof_CovFT0A2D_Cent_etaA_etaC_PiPlus"), cent, etaValA, etaValB, covFT0ARecoEffCor); - - if (std::isfinite(covFT0CTru)) - histos.fill(HIST("MCGen/Prof_CovFT0C2D_Cent_etaA_etaC_PiPlus"), cent, etaValA, etaValB, covFT0CTru); - if (std::isfinite(covFT0CReco)) - histos.fill(HIST("MCReco/Prof_CovFT0C2D_Cent_etaA_etaC_PiPlus"), cent, etaValA, etaValB, covFT0CReco); - if (std::isfinite(covFT0CRecoEffCor)) - histos.fill(HIST("MCRecoEffCorr/Prof_CovFT0C2D_Cent_etaA_etaC_PiPlus"), cent, etaValA, etaValB, covFT0CRecoEffCor); - - } else if (isp == kPiAllIdx) { - if (std::isfinite(c2SubTru)) { - histos.fill(HIST("MCGen/Prof_C2Sub2D_Cent_etaA_etaC_PiAll"), cent, etaValA, etaValB, c2SubTru); - histos.fill(HIST("MCGen/Prof_GapSum2D_PiAll"), cent, gap, sum, c2SubTru); - } - if (std::isfinite(c2SubReco)) { - histos.fill(HIST("MCReco/Prof_C2Sub2D_Cent_etaA_etaC_PiAll"), cent, etaValA, etaValB, c2SubReco); - histos.fill(HIST("MCReco/Prof_GapSum2D_PiAll"), cent, gap, sum, c2SubReco); - } - if (std::isfinite(c2SubRecoEffCor)) { - histos.fill(HIST("MCRecoEffCorr/Prof_C2Sub2D_Cent_etaA_etaC_PiAll"), cent, etaValA, etaValB, c2SubRecoEffCor); - histos.fill(HIST("MCRecoEffCorr/Prof_GapSum2D_PiAll"), cent, gap, sum, c2SubRecoEffCor); - } - - if (std::isfinite(covTru)) - histos.fill(HIST("MCGen/Prof_Cov2D_Cent_etaA_etaC_PiAll"), cent, etaValA, etaValB, covTru); - if (std::isfinite(covReco)) - histos.fill(HIST("MCReco/Prof_Cov2D_Cent_etaA_etaC_PiAll"), cent, etaValA, etaValB, covReco); - if (std::isfinite(covRecoEffCor)) - histos.fill(HIST("MCRecoEffCorr/Prof_Cov2D_Cent_etaA_etaC_PiAll"), cent, etaValA, etaValB, covRecoEffCor); - - if (std::isfinite(covFT0ATru)) - histos.fill(HIST("MCGen/Prof_CovFT0A2D_Cent_etaA_etaC_PiAll"), cent, etaValA, etaValB, covFT0ATru); - if (std::isfinite(covFT0AReco)) - histos.fill(HIST("MCReco/Prof_CovFT0A2D_Cent_etaA_etaC_PiAll"), cent, etaValA, etaValB, covFT0AReco); - if (std::isfinite(covFT0ARecoEffCor)) - histos.fill(HIST("MCRecoEffCorr/Prof_CovFT0A2D_Cent_etaA_etaC_PiAll"), cent, etaValA, etaValB, covFT0ARecoEffCor); - - if (std::isfinite(covFT0CTru)) - histos.fill(HIST("MCGen/Prof_CovFT0C2D_Cent_etaA_etaC_PiAll"), cent, etaValA, etaValB, covFT0CTru); - if (std::isfinite(covFT0CReco)) - histos.fill(HIST("MCReco/Prof_CovFT0C2D_Cent_etaA_etaC_PiAll"), cent, etaValA, etaValB, covFT0CReco); - if (std::isfinite(covFT0CRecoEffCor)) - histos.fill(HIST("MCRecoEffCorr/Prof_CovFT0C2D_Cent_etaA_etaC_PiAll"), cent, etaValA, etaValB, covFT0CRecoEffCor); - - } else if (isp == kKaMinusIdx) { - if (std::isfinite(c2SubTru)) { - histos.fill(HIST("MCGen/Prof_C2Sub2D_Cent_etaA_etaC_KaMinus"), cent, etaValA, etaValB, c2SubTru); - histos.fill(HIST("MCGen/Prof_GapSum2D_KaMinus"), cent, gap, sum, c2SubTru); - } - if (std::isfinite(c2SubReco)) { - histos.fill(HIST("MCReco/Prof_C2Sub2D_Cent_etaA_etaC_KaMinus"), cent, etaValA, etaValB, c2SubReco); - histos.fill(HIST("MCReco/Prof_GapSum2D_KaMinus"), cent, gap, sum, c2SubReco); - } - if (std::isfinite(c2SubRecoEffCor)) { - histos.fill(HIST("MCRecoEffCorr/Prof_C2Sub2D_Cent_etaA_etaC_KaMinus"), cent, etaValA, etaValB, c2SubRecoEffCor); - histos.fill(HIST("MCRecoEffCorr/Prof_GapSum2D_KaMinus"), cent, gap, sum, c2SubRecoEffCor); - } - - if (std::isfinite(covTru)) - histos.fill(HIST("MCGen/Prof_Cov2D_Cent_etaA_etaC_KaMinus"), cent, etaValA, etaValB, covTru); - if (std::isfinite(covReco)) - histos.fill(HIST("MCReco/Prof_Cov2D_Cent_etaA_etaC_KaMinus"), cent, etaValA, etaValB, covReco); - if (std::isfinite(covRecoEffCor)) - histos.fill(HIST("MCRecoEffCorr/Prof_Cov2D_Cent_etaA_etaC_KaMinus"), cent, etaValA, etaValB, covRecoEffCor); - - if (std::isfinite(covFT0ATru)) - histos.fill(HIST("MCGen/Prof_CovFT0A2D_Cent_etaA_etaC_KaMinus"), cent, etaValA, etaValB, covFT0ATru); - if (std::isfinite(covFT0AReco)) - histos.fill(HIST("MCReco/Prof_CovFT0A2D_Cent_etaA_etaC_KaMinus"), cent, etaValA, etaValB, covFT0AReco); - if (std::isfinite(covFT0ARecoEffCor)) - histos.fill(HIST("MCRecoEffCorr/Prof_CovFT0A2D_Cent_etaA_etaC_KaMinus"), cent, etaValA, etaValB, covFT0ARecoEffCor); - - if (std::isfinite(covFT0CTru)) - histos.fill(HIST("MCGen/Prof_CovFT0C2D_Cent_etaA_etaC_KaMinus"), cent, etaValA, etaValB, covFT0CTru); - if (std::isfinite(covFT0CReco)) - histos.fill(HIST("MCReco/Prof_CovFT0C2D_Cent_etaA_etaC_KaMinus"), cent, etaValA, etaValB, covFT0CReco); - if (std::isfinite(covFT0CRecoEffCor)) - histos.fill(HIST("MCRecoEffCorr/Prof_CovFT0C2D_Cent_etaA_etaC_KaMinus"), cent, etaValA, etaValB, covFT0CRecoEffCor); - - } else if (isp == kKaPlusIdx) { - if (std::isfinite(c2SubTru)) { - histos.fill(HIST("MCGen/Prof_C2Sub2D_Cent_etaA_etaC_KaPlus"), cent, etaValA, etaValB, c2SubTru); - histos.fill(HIST("MCGen/Prof_GapSum2D_KaPlus"), cent, gap, sum, c2SubTru); - } - if (std::isfinite(c2SubReco)) { - histos.fill(HIST("MCReco/Prof_C2Sub2D_Cent_etaA_etaC_KaPlus"), cent, etaValA, etaValB, c2SubReco); - histos.fill(HIST("MCReco/Prof_GapSum2D_KaPlus"), cent, gap, sum, c2SubReco); - } - if (std::isfinite(c2SubRecoEffCor)) { - histos.fill(HIST("MCRecoEffCorr/Prof_C2Sub2D_Cent_etaA_etaC_KaPlus"), cent, etaValA, etaValB, c2SubRecoEffCor); - histos.fill(HIST("MCRecoEffCorr/Prof_GapSum2D_KaPlus"), cent, gap, sum, c2SubRecoEffCor); - } - - if (std::isfinite(covTru)) - histos.fill(HIST("MCGen/Prof_Cov2D_Cent_etaA_etaC_KaPlus"), cent, etaValA, etaValB, covTru); - if (std::isfinite(covReco)) - histos.fill(HIST("MCReco/Prof_Cov2D_Cent_etaA_etaC_KaPlus"), cent, etaValA, etaValB, covReco); - if (std::isfinite(covRecoEffCor)) - histos.fill(HIST("MCRecoEffCorr/Prof_Cov2D_Cent_etaA_etaC_KaPlus"), cent, etaValA, etaValB, covRecoEffCor); - - if (std::isfinite(covFT0ATru)) - histos.fill(HIST("MCGen/Prof_CovFT0A2D_Cent_etaA_etaC_KaPlus"), cent, etaValA, etaValB, covFT0ATru); - if (std::isfinite(covFT0AReco)) - histos.fill(HIST("MCReco/Prof_CovFT0A2D_Cent_etaA_etaC_KaPlus"), cent, etaValA, etaValB, covFT0AReco); - if (std::isfinite(covFT0ARecoEffCor)) - histos.fill(HIST("MCRecoEffCorr/Prof_CovFT0A2D_Cent_etaA_etaC_KaPlus"), cent, etaValA, etaValB, covFT0ARecoEffCor); - - if (std::isfinite(covFT0CTru)) - histos.fill(HIST("MCGen/Prof_CovFT0C2D_Cent_etaA_etaC_KaPlus"), cent, etaValA, etaValB, covFT0CTru); - if (std::isfinite(covFT0CReco)) - histos.fill(HIST("MCReco/Prof_CovFT0C2D_Cent_etaA_etaC_KaPlus"), cent, etaValA, etaValB, covFT0CReco); - if (std::isfinite(covFT0CRecoEffCor)) - histos.fill(HIST("MCRecoEffCorr/Prof_CovFT0C2D_Cent_etaA_etaC_KaPlus"), cent, etaValA, etaValB, covFT0CRecoEffCor); - - } else if (isp == kKaAllIdx) { - if (std::isfinite(c2SubTru)) { - histos.fill(HIST("MCGen/Prof_C2Sub2D_Cent_etaA_etaC_KaAll"), cent, etaValA, etaValB, c2SubTru); - histos.fill(HIST("MCGen/Prof_GapSum2D_KaAll"), cent, gap, sum, c2SubTru); - } - if (std::isfinite(c2SubReco)) { - histos.fill(HIST("MCReco/Prof_C2Sub2D_Cent_etaA_etaC_KaAll"), cent, etaValA, etaValB, c2SubReco); - histos.fill(HIST("MCReco/Prof_GapSum2D_KaAll"), cent, gap, sum, c2SubReco); - } - if (std::isfinite(c2SubRecoEffCor)) { - histos.fill(HIST("MCRecoEffCorr/Prof_C2Sub2D_Cent_etaA_etaC_KaAll"), cent, etaValA, etaValB, c2SubRecoEffCor); - histos.fill(HIST("MCRecoEffCorr/Prof_GapSum2D_KaAll"), cent, gap, sum, c2SubRecoEffCor); - } - - if (std::isfinite(covTru)) - histos.fill(HIST("MCGen/Prof_Cov2D_Cent_etaA_etaC_KaAll"), cent, etaValA, etaValB, covTru); - if (std::isfinite(covReco)) - histos.fill(HIST("MCReco/Prof_Cov2D_Cent_etaA_etaC_KaAll"), cent, etaValA, etaValB, covReco); - if (std::isfinite(covRecoEffCor)) - histos.fill(HIST("MCRecoEffCorr/Prof_Cov2D_Cent_etaA_etaC_KaAll"), cent, etaValA, etaValB, covRecoEffCor); - - if (std::isfinite(covFT0ATru)) - histos.fill(HIST("MCGen/Prof_CovFT0A2D_Cent_etaA_etaC_KaAll"), cent, etaValA, etaValB, covFT0ATru); - if (std::isfinite(covFT0AReco)) - histos.fill(HIST("MCReco/Prof_CovFT0A2D_Cent_etaA_etaC_KaAll"), cent, etaValA, etaValB, covFT0AReco); - if (std::isfinite(covFT0ARecoEffCor)) - histos.fill(HIST("MCRecoEffCorr/Prof_CovFT0A2D_Cent_etaA_etaC_KaAll"), cent, etaValA, etaValB, covFT0ARecoEffCor); - - if (std::isfinite(covFT0CTru)) - histos.fill(HIST("MCGen/Prof_CovFT0C2D_Cent_etaA_etaC_KaAll"), cent, etaValA, etaValB, covFT0CTru); - if (std::isfinite(covFT0CReco)) - histos.fill(HIST("MCReco/Prof_CovFT0C2D_Cent_etaA_etaC_KaAll"), cent, etaValA, etaValB, covFT0CReco); - if (std::isfinite(covFT0CRecoEffCor)) - histos.fill(HIST("MCRecoEffCorr/Prof_CovFT0C2D_Cent_etaA_etaC_KaAll"), cent, etaValA, etaValB, covFT0CRecoEffCor); - - } else if (isp == kPrIdx) { - if (std::isfinite(c2SubTru)) { - histos.fill(HIST("MCGen/Prof_C2Sub2D_Cent_etaA_etaC_Pr"), cent, etaValA, etaValB, c2SubTru); - histos.fill(HIST("MCGen/Prof_GapSum2D_Pr"), cent, gap, sum, c2SubTru); - } - if (std::isfinite(c2SubReco)) { - histos.fill(HIST("MCReco/Prof_C2Sub2D_Cent_etaA_etaC_Pr"), cent, etaValA, etaValB, c2SubReco); - histos.fill(HIST("MCReco/Prof_GapSum2D_Pr"), cent, gap, sum, c2SubReco); - } - if (std::isfinite(c2SubRecoEffCor)) { - histos.fill(HIST("MCRecoEffCorr/Prof_C2Sub2D_Cent_etaA_etaC_Pr"), cent, etaValA, etaValB, c2SubRecoEffCor); - histos.fill(HIST("MCRecoEffCorr/Prof_GapSum2D_Pr"), cent, gap, sum, c2SubRecoEffCor); - } - - if (std::isfinite(covTru)) - histos.fill(HIST("MCGen/Prof_Cov2D_Cent_etaA_etaC_Pr"), cent, etaValA, etaValB, covTru); - if (std::isfinite(covReco)) - histos.fill(HIST("MCReco/Prof_Cov2D_Cent_etaA_etaC_Pr"), cent, etaValA, etaValB, covReco); - if (std::isfinite(covRecoEffCor)) - histos.fill(HIST("MCRecoEffCorr/Prof_Cov2D_Cent_etaA_etaC_Pr"), cent, etaValA, etaValB, covRecoEffCor); - - if (std::isfinite(covFT0ATru)) - histos.fill(HIST("MCGen/Prof_CovFT0A2D_Cent_etaA_etaC_Pr"), cent, etaValA, etaValB, covFT0ATru); - if (std::isfinite(covFT0AReco)) - histos.fill(HIST("MCReco/Prof_CovFT0A2D_Cent_etaA_etaC_Pr"), cent, etaValA, etaValB, covFT0AReco); - if (std::isfinite(covFT0ARecoEffCor)) - histos.fill(HIST("MCRecoEffCorr/Prof_CovFT0A2D_Cent_etaA_etaC_Pr"), cent, etaValA, etaValB, covFT0ARecoEffCor); - - if (std::isfinite(covFT0CTru)) - histos.fill(HIST("MCGen/Prof_CovFT0C2D_Cent_etaA_etaC_Pr"), cent, etaValA, etaValB, covFT0CTru); - if (std::isfinite(covFT0CReco)) - histos.fill(HIST("MCReco/Prof_CovFT0C2D_Cent_etaA_etaC_Pr"), cent, etaValA, etaValB, covFT0CReco); - if (std::isfinite(covFT0CRecoEffCor)) - histos.fill(HIST("MCRecoEffCorr/Prof_CovFT0C2D_Cent_etaA_etaC_Pr"), cent, etaValA, etaValB, covFT0CRecoEffCor); - - } else if (isp == kAntiPrIdx) { - if (std::isfinite(c2SubTru)) { - histos.fill(HIST("MCGen/Prof_C2Sub2D_Cent_etaA_etaC_AntiPr"), cent, etaValA, etaValB, c2SubTru); - histos.fill(HIST("MCGen/Prof_GapSum2D_AntiPr"), cent, gap, sum, c2SubTru); - } - if (std::isfinite(c2SubReco)) { - histos.fill(HIST("MCReco/Prof_C2Sub2D_Cent_etaA_etaC_AntiPr"), cent, etaValA, etaValB, c2SubReco); - histos.fill(HIST("MCReco/Prof_GapSum2D_AntiPr"), cent, gap, sum, c2SubReco); - } - if (std::isfinite(c2SubRecoEffCor)) { - histos.fill(HIST("MCRecoEffCorr/Prof_C2Sub2D_Cent_etaA_etaC_AntiPr"), cent, etaValA, etaValB, c2SubRecoEffCor); - histos.fill(HIST("MCRecoEffCorr/Prof_GapSum2D_AntiPr"), cent, gap, sum, c2SubRecoEffCor); - } - - if (std::isfinite(covTru)) - histos.fill(HIST("MCGen/Prof_Cov2D_Cent_etaA_etaC_AntiPr"), cent, etaValA, etaValB, covTru); - if (std::isfinite(covReco)) - histos.fill(HIST("MCReco/Prof_Cov2D_Cent_etaA_etaC_AntiPr"), cent, etaValA, etaValB, covReco); - if (std::isfinite(covRecoEffCor)) - histos.fill(HIST("MCRecoEffCorr/Prof_Cov2D_Cent_etaA_etaC_AntiPr"), cent, etaValA, etaValB, covRecoEffCor); - - if (std::isfinite(covFT0ATru)) - histos.fill(HIST("MCGen/Prof_CovFT0A2D_Cent_etaA_etaC_AntiPr"), cent, etaValA, etaValB, covFT0ATru); - if (std::isfinite(covFT0AReco)) - histos.fill(HIST("MCReco/Prof_CovFT0A2D_Cent_etaA_etaC_AntiPr"), cent, etaValA, etaValB, covFT0AReco); - if (std::isfinite(covFT0ARecoEffCor)) - histos.fill(HIST("MCRecoEffCorr/Prof_CovFT0A2D_Cent_etaA_etaC_AntiPr"), cent, etaValA, etaValB, covFT0ARecoEffCor); - - if (std::isfinite(covFT0CTru)) - histos.fill(HIST("MCGen/Prof_CovFT0C2D_Cent_etaA_etaC_AntiPr"), cent, etaValA, etaValB, covFT0CTru); - if (std::isfinite(covFT0CReco)) - histos.fill(HIST("MCReco/Prof_CovFT0C2D_Cent_etaA_etaC_AntiPr"), cent, etaValA, etaValB, covFT0CReco); - if (std::isfinite(covFT0CRecoEffCor)) - histos.fill(HIST("MCRecoEffCorr/Prof_CovFT0C2D_Cent_etaA_etaC_AntiPr"), cent, etaValA, etaValB, covFT0CRecoEffCor); - - } else if (isp == kPrAllIdx) { - if (std::isfinite(c2SubTru)) { - histos.fill(HIST("MCGen/Prof_C2Sub2D_Cent_etaA_etaC_PrAll"), cent, etaValA, etaValB, c2SubTru); - histos.fill(HIST("MCGen/Prof_GapSum2D_PrAll"), cent, gap, sum, c2SubTru); - } - if (std::isfinite(c2SubReco)) { - histos.fill(HIST("MCReco/Prof_C2Sub2D_Cent_etaA_etaC_PrAll"), cent, etaValA, etaValB, c2SubReco); - histos.fill(HIST("MCReco/Prof_GapSum2D_PrAll"), cent, gap, sum, c2SubReco); - } - if (std::isfinite(c2SubRecoEffCor)) { - histos.fill(HIST("MCRecoEffCorr/Prof_C2Sub2D_Cent_etaA_etaC_PrAll"), cent, etaValA, etaValB, c2SubRecoEffCor); - histos.fill(HIST("MCRecoEffCorr/Prof_GapSum2D_PrAll"), cent, gap, sum, c2SubRecoEffCor); - } - - if (std::isfinite(covTru)) - histos.fill(HIST("MCGen/Prof_Cov2D_Cent_etaA_etaC_PrAll"), cent, etaValA, etaValB, covTru); - if (std::isfinite(covReco)) - histos.fill(HIST("MCReco/Prof_Cov2D_Cent_etaA_etaC_PrAll"), cent, etaValA, etaValB, covReco); - if (std::isfinite(covRecoEffCor)) - histos.fill(HIST("MCRecoEffCorr/Prof_Cov2D_Cent_etaA_etaC_PrAll"), cent, etaValA, etaValB, covRecoEffCor); - - if (std::isfinite(covFT0ATru)) - histos.fill(HIST("MCGen/Prof_CovFT0A2D_Cent_etaA_etaC_PrAll"), cent, etaValA, etaValB, covFT0ATru); - if (std::isfinite(covFT0AReco)) - histos.fill(HIST("MCReco/Prof_CovFT0A2D_Cent_etaA_etaC_PrAll"), cent, etaValA, etaValB, covFT0AReco); - if (std::isfinite(covFT0ARecoEffCor)) - histos.fill(HIST("MCRecoEffCorr/Prof_CovFT0A2D_Cent_etaA_etaC_PrAll"), cent, etaValA, etaValB, covFT0ARecoEffCor); - - if (std::isfinite(covFT0CTru)) - histos.fill(HIST("MCGen/Prof_CovFT0C2D_Cent_etaA_etaC_PrAll"), cent, etaValA, etaValB, covFT0CTru); - if (std::isfinite(covFT0CReco)) - histos.fill(HIST("MCReco/Prof_CovFT0C2D_Cent_etaA_etaC_PrAll"), cent, etaValA, etaValB, covFT0CReco); - if (std::isfinite(covFT0CRecoEffCor)) - histos.fill(HIST("MCRecoEffCorr/Prof_CovFT0C2D_Cent_etaA_etaC_PrAll"), cent, etaValA, etaValB, covFT0CRecoEffCor); - } + histos.fill(HIST("MCRecoEffCorr/Prof_CovFT0C2D_Cent_etaA_etaC"), cent, etaValA, etaValB, covFT0CRecoEffCor); } } } } PROCESS_SWITCH(RadialFlowDecorr, processMCFluc, "process MC to calculate pt fluc", cfgRunMCFluc); - void processDataGetNSig(AodCollisionsSel::iterator const& coll, BCsRun3 const& /*bcs*/, aod::Zdcs const& /*zdcsData*/, AodTracksSel const& tracks) - { - histos.fill(HIST("hVtxZ"), coll.posZ()); - if (!isEventSelected(coll)) - return; - float cent = getCentrality(coll); - if (cent > KCentMax) - return; - histos.fill(HIST("hVtxZ_after_sel"), coll.posZ()); - histos.fill(HIST("hCentrality"), cent); - - histos.fill(HIST("Hist2D_globalTracks_PVTracks"), coll.multNTracksPV(), tracks.size()); - histos.fill(HIST("Hist2D_cent_nch"), tracks.size(), cent); - histos.fill(HIST("Hist2D_globalTracks_cent"), cent, tracks.size()); - histos.fill(HIST("Hist2D_PVTracks_cent"), cent, coll.multNTracksPV()); - - int ntrk = 0; - for (const auto& track : tracks) { - if (!isTrackSelected(track)) - continue; - float pt = track.pt(); - if (pt <= cfgPtMin || pt > cfgPtMax) - continue; - float eta = track.eta(); - if (eta > etaLw[0] && eta < etaUp[0]) - ntrk++; - fillNSigmaBefCut(track, cent); - } - - if (cfgZDC) { - const auto& foundBC = coll.foundBC_as(); - if (!foundBC.has_zdc()) { - return; - } - auto zdc = foundBC.zdc(); - auto zdcAmp = zdc.energyCommonZNA() + zdc.energyCommonZNC(); - histos.fill(HIST("hnTrkPVZDC"), coll.multNTracksPV(), zdcAmp); - histos.fill(HIST("hNchZDC"), ntrk, zdcAmp); - } - } - PROCESS_SWITCH(RadialFlowDecorr, processDataGetNSig, "process data to Get Nsigma cuts", cfgRunDataGetNSig); - + // =========================================================================== + // DATA: build flattening maps (inclusive) + // =========================================================================== void processGetDataFlat(AodCollisionsSel::iterator const& coll, BCsRun3 const& /*bcs*/, aod::Zdcs const& /*zdcsData*/, AodTracksSel const& tracks) { histos.fill(HIST("hVtxZ"), coll.posZ()); - if (!isEventSelected(coll)) + if (!isEventSelected(coll)) { return; + } float cent = getCentrality(coll); - if (cent > KCentMax) + if (cent > KCentMax) { return; - - if (!isPassAddPileup(coll.multNTracksPV(), tracks.size(), cent)) + } + if (!isPassAddPileup(coll.multNTracksPV(), tracks.size(), cent)) { return; + } histos.fill(HIST("hVtxZ_after_sel"), coll.posZ()); histos.fill(HIST("hCentrality"), cent); - histos.fill(HIST("Hist2D_globalTracks_PVTracks"), coll.multNTracksPV(), tracks.size()); histos.fill(HIST("Hist2D_cent_nch"), tracks.size(), cent); histos.fill(HIST("Hist2D_globalTracks_cent"), cent, tracks.size()); @@ -3016,88 +1864,36 @@ struct RadialFlowDecorr { float vz = coll.posZ(); for (const auto& track : tracks) { - if (!isTrackSelected(track)) + if (!isTrackSelected(track)) { continue; + } float pt = track.pt(); - if (pt <= cfgPtMin || pt > cfgPtMax) + if (pt <= cfgPtMin || pt > cfgPtMax) { continue; - float eta = track.eta(); - float phi = track.phi(); + } + float eta = track.eta(), phi = track.phi(); auto sign = track.sign(); histos.fill(HIST("hPt"), pt); histos.fill(HIST("hEta"), eta); histos.fill(HIST("hPhi"), phi); - - if (eta > etaLw[0] && eta < etaUp[0]) + if (eta > etaLw[0] && eta < etaUp[0]) { ntrk++; - fillNSigmaBefCut(track, cent); - int id = identifyTrack(track, cent); - bool isPi = (id == KPidPionOne); - bool isKa = (id == KPidKaonTwo); - bool isPr = (id == KPidProtonThree); - bool isSpecies[KNsp] = { - true, - isPi && sign < 0, isPi && sign > 0, isPi, - isKa && sign < 0, isKa && sign > 0, isKa, - isPr && sign < 0, isPr && sign > 0, isPr}; - - fillNSigmaAftCut(track, cent, isSpecies); - for (int isp = 0; isp < KNsp; ++isp) { - if (!isSpecies[isp]) - continue; - float eff = getEfficiency(coll.multNTracksPV(), pt, eta, static_cast(isp), 0, cfgEff); - if (eff <= KFloatEpsilon) - continue; - - float fake = getEfficiency(coll.multNTracksPV(), pt, eta, static_cast(isp), 1, cfgEff); - float w = (1.0f - fake) / eff; - - if (!std::isfinite(w) || w <= 0.f) - continue; + } - if (isp == kInclusiveIdx) { - histos.fill(HIST("hEtaPhiReco"), vz, sign, pt, eta, phi); - histos.fill(HIST("hEtaPhiRecoEffWtd"), vz, sign, pt, eta, phi, (1.0f - fake) / eff); - histos.fill(HIST("hEtaPhiRecoWtd"), vz, sign, pt, eta, phi, w); - } else if (isp == kPiMinusIdx) { - histos.fill(HIST("hEtaPhiReco_PiMinus"), vz, sign, pt, eta, phi); - histos.fill(HIST("hEtaPhiRecoEffWtd_PiMinus"), vz, sign, pt, eta, phi, (1.0f - fake) / eff); - histos.fill(HIST("hEtaPhiRecoWtd_PiMinus"), vz, sign, pt, eta, phi, w); - } else if (isp == kPiPlusIdx) { - histos.fill(HIST("hEtaPhiReco_PiPlus"), vz, sign, pt, eta, phi); - histos.fill(HIST("hEtaPhiRecoEffWtd_PiPlus"), vz, sign, pt, eta, phi, (1.0f - fake) / eff); - histos.fill(HIST("hEtaPhiRecoWtd_PiPlus"), vz, sign, pt, eta, phi, w); - } else if (isp == kPiAllIdx) { - histos.fill(HIST("hEtaPhiReco_PiAll"), vz, sign, pt, eta, phi); - histos.fill(HIST("hEtaPhiRecoEffWtd_PiAll"), vz, sign, pt, eta, phi, (1.0f - fake) / eff); - histos.fill(HIST("hEtaPhiRecoWtd_PiAll"), vz, sign, pt, eta, phi, w); - } else if (isp == kKaMinusIdx) { - histos.fill(HIST("hEtaPhiReco_KaMinus"), vz, sign, pt, eta, phi); - histos.fill(HIST("hEtaPhiRecoEffWtd_KaMinus"), vz, sign, pt, eta, phi, (1.0f - fake) / eff); - histos.fill(HIST("hEtaPhiRecoWtd_KaMinus"), vz, sign, pt, eta, phi, w); - } else if (isp == kKaPlusIdx) { - histos.fill(HIST("hEtaPhiReco_KaPlus"), vz, sign, pt, eta, phi); - histos.fill(HIST("hEtaPhiRecoEffWtd_KaPlus"), vz, sign, pt, eta, phi, (1.0f - fake) / eff); - histos.fill(HIST("hEtaPhiRecoWtd_KaPlus"), vz, sign, pt, eta, phi, w); - } else if (isp == kKaAllIdx) { - histos.fill(HIST("hEtaPhiReco_KaAll"), vz, sign, pt, eta, phi); - histos.fill(HIST("hEtaPhiRecoEffWtd_KaAll"), vz, sign, pt, eta, phi, (1.0f - fake) / eff); - histos.fill(HIST("hEtaPhiRecoWtd_KaAll"), vz, sign, pt, eta, phi, w); - } else if (isp == kPrIdx) { - histos.fill(HIST("hEtaPhiReco_Pr"), vz, sign, pt, eta, phi); - histos.fill(HIST("hEtaPhiRecoEffWtd_Pr"), vz, sign, pt, eta, phi, (1.0f - fake) / eff); - histos.fill(HIST("hEtaPhiRecoWtd_Pr"), vz, sign, pt, eta, phi, w); - } else if (isp == kAntiPrIdx) { - histos.fill(HIST("hEtaPhiReco_AntiPr"), vz, sign, pt, eta, phi); - histos.fill(HIST("hEtaPhiRecoEffWtd_AntiPr"), vz, sign, pt, eta, phi, (1.0f - fake) / eff); - histos.fill(HIST("hEtaPhiRecoWtd_AntiPr"), vz, sign, pt, eta, phi, w); - } else if (isp == kPrAllIdx) { - histos.fill(HIST("hEtaPhiReco_PrAll"), vz, sign, pt, eta, phi); - histos.fill(HIST("hEtaPhiRecoEffWtd_PrAll"), vz, sign, pt, eta, phi, (1.0f - fake) / eff); - histos.fill(HIST("hEtaPhiRecoWtd_PrAll"), vz, sign, pt, eta, phi, w); - } + float eff = getEfficiency(coll.multNTracksPV(), pt, eta, 0, cfgEff); + if (eff <= KFloatEpsilon) { + continue; + } + float fake = getEfficiency(coll.multNTracksPV(), pt, eta, 1, cfgEff); + float w = (1.0f - fake) / eff; + if (!std::isfinite(w) || w <= 0.f) { + continue; } + + histos.fill(HIST("hEtaPhiReco"), vz, sign, pt, eta, phi); + histos.fill(HIST("hEtaPhiRecoEffWtd"), vz, sign, pt, eta, phi, (1.0f - fake) / eff); + histos.fill(HIST("hEtaPhiRecoWtd"), vz, sign, pt, eta, phi, w); } if (cfgZDC) { @@ -3113,23 +1909,26 @@ struct RadialFlowDecorr { } PROCESS_SWITCH(RadialFlowDecorr, processGetDataFlat, "process data to calculate Flattening maps", cfgRunGetDataFlat); + // =========================================================================== + // DATA: mean pT + // =========================================================================== void processDataMean(AodCollisionsSel::iterator const& coll, BCsRun3 const& /*bcs*/, aod::Zdcs const& /*zdcsData*/, aod::FT0s const&, AodTracksSel const& tracks) { - double sumWi[KNsp][KNEta]{}, sumWipti[KNsp][KNEta]{}; + std::array sumWi{}, sumWipti{}; - if (!isEventSelected(coll)) + if (!isEventSelected(coll)) { return; - + } float cent = getCentrality(coll); - if (cent > KCentMax) + if (cent > KCentMax) { return; - - if (!isPassAddPileup(coll.multNTracksPV(), tracks.size(), cent)) + } + if (!isPassAddPileup(coll.multNTracksPV(), tracks.size(), cent)) { return; + } histos.fill(HIST("hVtxZ_after_sel"), coll.posZ()); histos.fill(HIST("hCentrality"), cent); - histos.fill(HIST("Hist2D_globalTracks_PVTracks"), coll.multNTracksPV(), tracks.size()); histos.fill(HIST("Hist2D_cent_nch"), tracks.size(), cent); histos.fill(HIST("Hist2D_globalTracks_cent"), cent, tracks.size()); @@ -3138,154 +1937,88 @@ struct RadialFlowDecorr { float vz = coll.posZ(); for (const auto& track : tracks) { - if (!isTrackSelected(track)) + if (!isTrackSelected(track)) { continue; - + } float p = track.p(); float pt = track.pt(); - float eta = track.eta(); - float phi = track.phi(); + float eta = track.eta(), phi = track.phi(); auto sign = track.sign(); - - if (p < KFloatEpsilon) + if (p < KFloatEpsilon) { continue; - - if (pt <= cfgPtMin || pt > cfgPtMax) + } + if (pt <= cfgPtMin || pt > cfgPtMax) { continue; + } histos.fill(HIST("hP"), p); histos.fill(HIST("hPt"), pt); histos.fill(HIST("hEta"), eta); histos.fill(HIST("hPhi"), phi); - int id = identifyTrack(track, cent); - bool isPi = (id == KPidPionOne); - bool isKa = (id == KPidKaonTwo); - bool isPr = (id == KPidProtonThree); - bool isSpecies[KNsp] = { - true, - isPi && sign < 0, isPi && sign > 0, isPi, - isKa && sign < 0, isKa && sign > 0, isKa, - isPr && sign < 0, isPr && sign > 0, isPr}; - - for (int isp = 0; isp < KNsp; ++isp) { - if (!isSpecies[isp]) - continue; - float eff = getEfficiency(coll.multNTracksPV(), pt, eta, static_cast(isp), 0, cfgEff); - if (eff <= KFloatEpsilon) - continue; - float fake = getEfficiency(coll.multNTracksPV(), pt, eta, static_cast(isp), 1, cfgEff); - float flatWeight = getFlatteningWeight(vz, sign, pt, eta, phi, static_cast(isp), cfgFlat); - float w = flatWeight * (1.0f - fake) / eff; + float eff = getEfficiency(coll.multNTracksPV(), pt, eta, 0, cfgEff); + float fake = getEfficiency(coll.multNTracksPV(), pt, eta, 1, cfgEff); + float flatWeight = getFlatteningWeight(vz, sign, pt, eta, phi, cfgFlat); - if (!std::isfinite(w) || w <= 0.f) - continue; + histos.fill(HIST("pEffWeight_pt_eta_cent"), pt, eta, cent, eff); + histos.fill(HIST("pFakeWeight_pt_eta_cent"), pt, eta, cent, fake); + histos.fill(HIST("pFlatWeight_pt_eta_cent"), pt, eta, cent, flatWeight); - if (isp == kInclusiveIdx) { - histos.fill(HIST("hEtaPhiReco"), vz, sign, pt, eta, phi); - histos.fill(HIST("hEtaPhiRecoEffWtd"), vz, sign, pt, eta, phi, (1.0f - fake) / eff); - histos.fill(HIST("hEtaPhiRecoWtd"), vz, sign, pt, eta, phi, w); - } else if (isp == kPiMinusIdx) { - histos.fill(HIST("hEtaPhiReco_PiMinus"), vz, sign, pt, eta, phi); - histos.fill(HIST("hEtaPhiRecoEffWtd_PiMinus"), vz, sign, pt, eta, phi, (1.0f - fake) / eff); - histos.fill(HIST("hEtaPhiRecoWtd_PiMinus"), vz, sign, pt, eta, phi, w); - } else if (isp == kPiPlusIdx) { - histos.fill(HIST("hEtaPhiReco_PiPlus"), vz, sign, pt, eta, phi); - histos.fill(HIST("hEtaPhiRecoEffWtd_PiPlus"), vz, sign, pt, eta, phi, (1.0f - fake) / eff); - histos.fill(HIST("hEtaPhiRecoWtd_PiPlus"), vz, sign, pt, eta, phi, w); - } else if (isp == kPiAllIdx) { - histos.fill(HIST("hEtaPhiReco_PiAll"), vz, sign, pt, eta, phi); - histos.fill(HIST("hEtaPhiRecoEffWtd_PiAll"), vz, sign, pt, eta, phi, (1.0f - fake) / eff); - histos.fill(HIST("hEtaPhiRecoWtd_PiAll"), vz, sign, pt, eta, phi, w); - } else if (isp == kKaMinusIdx) { - histos.fill(HIST("hEtaPhiReco_KaMinus"), vz, sign, pt, eta, phi); - histos.fill(HIST("hEtaPhiRecoEffWtd_KaMinus"), vz, sign, pt, eta, phi, (1.0f - fake) / eff); - histos.fill(HIST("hEtaPhiRecoWtd_KaMinus"), vz, sign, pt, eta, phi, w); - } else if (isp == kKaPlusIdx) { - histos.fill(HIST("hEtaPhiReco_KaPlus"), vz, sign, pt, eta, phi); - histos.fill(HIST("hEtaPhiRecoEffWtd_KaPlus"), vz, sign, pt, eta, phi, (1.0f - fake) / eff); - histos.fill(HIST("hEtaPhiRecoWtd_KaPlus"), vz, sign, pt, eta, phi, w); - } else if (isp == kKaAllIdx) { - histos.fill(HIST("hEtaPhiReco_KaAll"), vz, sign, pt, eta, phi); - histos.fill(HIST("hEtaPhiRecoEffWtd_KaAll"), vz, sign, pt, eta, phi, (1.0f - fake) / eff); - histos.fill(HIST("hEtaPhiRecoWtd_KaAll"), vz, sign, pt, eta, phi, w); - } else if (isp == kPrIdx) { - histos.fill(HIST("hEtaPhiReco_Pr"), vz, sign, pt, eta, phi); - histos.fill(HIST("hEtaPhiRecoEffWtd_Pr"), vz, sign, pt, eta, phi, (1.0f - fake) / eff); - histos.fill(HIST("hEtaPhiRecoWtd_Pr"), vz, sign, pt, eta, phi, w); - } else if (isp == kAntiPrIdx) { - histos.fill(HIST("hEtaPhiReco_AntiPr"), vz, sign, pt, eta, phi); - histos.fill(HIST("hEtaPhiRecoEffWtd_AntiPr"), vz, sign, pt, eta, phi, (1.0f - fake) / eff); - histos.fill(HIST("hEtaPhiRecoWtd_AntiPr"), vz, sign, pt, eta, phi, w); - } else if (isp == kPrAllIdx) { - histos.fill(HIST("hEtaPhiReco_PrAll"), vz, sign, pt, eta, phi); - histos.fill(HIST("hEtaPhiRecoEffWtd_PrAll"), vz, sign, pt, eta, phi, (1.0f - fake) / eff); - histos.fill(HIST("hEtaPhiRecoWtd_PrAll"), vz, sign, pt, eta, phi, w); - } + if (eff <= KFloatEpsilon) { + continue; + } - for (int ieta = 0; ieta < KNEta; ++ieta) { - if (eta <= etaLw[ieta] || eta > etaUp[ieta]) - continue; - sumWi[isp][ieta] += w; - sumWipti[isp][ieta] += w * pt; + float w = flatWeight * (1.0f - fake) / eff; + if (!std::isfinite(w) || w <= 0.f) { + continue; + } + + histos.fill(HIST("hEtaPhiReco"), vz, sign, pt, eta, phi); + histos.fill(HIST("hEtaPhiRecoEffWtd"), vz, sign, pt, eta, phi, (1.0f - fake) / eff); + histos.fill(HIST("hEtaPhiRecoWtd"), vz, sign, pt, eta, phi, w); + + for (int ieta = 0; ieta < nEta; ++ieta) { + if (eta <= etaLw[ieta] || eta > etaUp[ieta]) { + continue; } + sumWi[ieta] += w; + sumWipti[ieta] += w * pt; } } - if (!hasMinTracksInAllEtaBins(sumWi)) + if (!hasMinTracksInAllEtaBins(sumWi)) { return; + } - for (int isp = 0; isp < KNsp; ++isp) { - if (sumWi[isp][0] < 1.0f) - continue; - histos.fill(HIST("Prof_Cent_Nsp_Nchrec"), cent, isp, sumWi[isp][0]); - histos.fill(HIST("Prof_Mult_Nsp_Nchrec"), coll.multNTracksPV(), isp, sumWi[isp][0]); - histos.fill(HIST("Prof_Cent_Nsp_MeanpT"), cent, isp, sumWipti[isp][0] / sumWi[isp][0]); - histos.fill(HIST("Prof_Mult_Nsp_MeanpT"), coll.multNTracksPV(), isp, sumWipti[isp][0] / sumWi[isp][0]); - } - - for (int ietaA = 0; ietaA < KNEta; ++ietaA) { - for (int ietaC = 0; ietaC < KNEta; ++ietaC) { - for (int isp = 0; isp < KNsp; ++isp) { - if ((sumWi[isp][ietaA] < 1.0f) || (sumWi[isp][ietaC] < 1.0f)) - continue; - - double wCorrAB = sumWi[isp][ietaA] + sumWi[isp][ietaC]; - if (wCorrAB > 0) { - float mptsub = (sumWipti[isp][ietaA] + sumWipti[isp][ietaC]) / wCorrAB; - if (isp == kInclusiveIdx) - histos.fill(HIST("Prof2D_MeanpTSub"), cent, ietaA, ietaC, mptsub); - else if (isp == kPiMinusIdx) - histos.fill(HIST("Prof2D_MeanpTSub_PiMinus"), cent, ietaA, ietaC, mptsub); - else if (isp == kPiPlusIdx) - histos.fill(HIST("Prof2D_MeanpTSub_PiPlus"), cent, ietaA, ietaC, mptsub); - else if (isp == kPiAllIdx) - histos.fill(HIST("Prof2D_MeanpTSub_PiAll"), cent, ietaA, ietaC, mptsub); - else if (isp == kKaMinusIdx) - histos.fill(HIST("Prof2D_MeanpTSub_KaMinus"), cent, ietaA, ietaC, mptsub); - else if (isp == kKaPlusIdx) - histos.fill(HIST("Prof2D_MeanpTSub_KaPlus"), cent, ietaA, ietaC, mptsub); - else if (isp == kKaAllIdx) - histos.fill(HIST("Prof2D_MeanpTSub_KaAll"), cent, ietaA, ietaC, mptsub); - else if (isp == kPrIdx) - histos.fill(HIST("Prof2D_MeanpTSub_Pr"), cent, ietaA, ietaC, mptsub); - else if (isp == kAntiPrIdx) - histos.fill(HIST("Prof2D_MeanpTSub_AntiPr"), cent, ietaA, ietaC, mptsub); - else if (isp == kPrAllIdx) - histos.fill(HIST("Prof2D_MeanpTSub_PrAll"), cent, ietaA, ietaC, mptsub); - } - if (ietaA == ietaC) { - double mpt = sumWipti[isp][ietaA] / sumWi[isp][ietaA]; - if (sumWi[isp][ietaA] >= 1.0f && std::isfinite(mpt)) { - histos.fill(HIST("pmean_nch_etabin_spbin"), coll.multNTracksPV(), ietaA, isp, mpt); - histos.fill(HIST("pmeanMult_nch_etabin_spbin"), coll.multNTracksPV(), ietaA, isp, sumWi[isp][ietaA]); - histos.fill(HIST("pmean_cent_etabin_spbin"), cent, ietaA, isp, mpt); - histos.fill(HIST("pmeanMult_cent_etabin_spbin"), cent, ietaA, isp, sumWi[isp][ietaA]); - } - } + if (sumWi[0] >= 1.0f) { + histos.fill(HIST("Prof_Cent_Nchrec"), cent, sumWi[0]); + histos.fill(HIST("Prof_Mult_Nchrec"), coll.multNTracksPV(), sumWi[0]); + histos.fill(HIST("Prof_Cent_MeanpT"), cent, sumWipti[0] / sumWi[0]); + histos.fill(HIST("Prof_Mult_MeanpT"), coll.multNTracksPV(), sumWipti[0] / sumWi[0]); + } + + for (int ietaA = 0; ietaA < nEta; ++ietaA) { + for (int ietaC = 0; ietaC < nEta; ++ietaC) { + if ((sumWi[ietaA] < 1.0f) || (sumWi[ietaC] < 1.0f)) { + continue; + } + double wCorrAB = sumWi[ietaA] + sumWi[ietaC]; + if (wCorrAB > 0) { + float mptsub = (sumWipti[ietaA] + sumWipti[ietaC]) / wCorrAB; + histos.fill(HIST("Prof2D_MeanpTSub"), cent, ietaA, ietaC, mptsub); + } + } + if (sumWi[ietaA] >= 1.0f) { + double mpt = sumWipti[ietaA] / sumWi[ietaA]; + if (std::isfinite(mpt)) { + histos.fill(HIST("pmean_nch_etabin"), coll.multNTracksPV(), ietaA, mpt); + histos.fill(HIST("pmeanMult_nch_etabin"), coll.multNTracksPV(), ietaA, sumWi[ietaA]); + histos.fill(HIST("pmean_cent_etabin"), cent, ietaA, mpt); + histos.fill(HIST("pmeanMult_cent_etabin"), cent, ietaA, sumWi[ietaA]); } } } + double amplFT0A = 0, amplFT0C = 0; if (coll.has_foundFT0()) { const auto& ft0 = coll.foundFT0(); @@ -3301,13 +2034,12 @@ struct RadialFlowDecorr { auto chanelid = ft0.channelC()[iCh]; auto globalId = chanelid + KnFt0cCell; float ampl = ft0.amplitudeC()[iCh]; - auto eta = getEtaFT0(globalId, 1); amplFT0C += ampl; + auto eta = getEtaFT0(globalId, 1); histos.fill(HIST("pmean_cent_id_eta_FT0"), cent, globalId, eta, ampl); histos.fill(HIST("h3_cent_id_eta_FT0"), cent, globalId, eta, ampl); } } - histos.fill(HIST("pmeanFT0Amultpv"), coll.multNTracksPV(), amplFT0A); histos.fill(HIST("pmeanFT0A_cent"), cent, amplFT0A); histos.fill(HIST("pmeanFT0Cmultpv"), coll.multNTracksPV(), amplFT0C); @@ -3315,310 +2047,233 @@ struct RadialFlowDecorr { } PROCESS_SWITCH(RadialFlowDecorr, processDataMean, "process data to calculate mean pT", cfgRunDataMean); + // =========================================================================== + // DATA: fluctuations (C2, subevent) + Poisson bootstrap for the base run + // =========================================================================== void processDataFluc(AodCollisionsSel::iterator const& coll, BCsRun3 const& /*bcs*/, aod::Zdcs const& /*zdcsData*/, aod::FT0s const&, AodTracksSel const& tracks) { - if (!isEventSelected(coll)) + if (!isEventSelected(coll)) { return; + } float cent = getCentrality(coll); - if (cent > KCentMax) + if (cent > KCentMax) { return; - - if (!isPassAddPileup(coll.multNTracksPV(), tracks.size(), cent)) + } + if (!isPassAddPileup(coll.multNTracksPV(), tracks.size(), cent)) { return; + } histos.fill(HIST("hVtxZ_after_sel"), coll.posZ()); histos.fill(HIST("hCentrality"), cent); - histos.fill(HIST("Hist2D_globalTracks_PVTracks"), coll.multNTracksPV(), tracks.size()); histos.fill(HIST("Hist2D_cent_nch"), tracks.size(), cent); histos.fill(HIST("Hist2D_globalTracks_cent"), cent, tracks.size()); histos.fill(HIST("Hist2D_PVTracks_cent"), cent, coll.multNTracksPV()); - if (!state.pmeanNchEtabinSpbinStep2 || !state.pmeanMultNchEtabinSpbinStep2) { - LOGF(warning, "Data fluc: Mean pT or Mult map missing"); + if (!state.pmeanNchEtabinStep2 || !state.pmeanMultNchEtabinStep2) { + LOGF(warning, "Data fluc: mean pT or mult map missing"); return; } - - for (int isp = 0; isp < KNsp; ++isp) { - auto pid = static_cast(isp); - if (!state.hEff[pid] || !state.hFake[pid] || !state.hFlatWeight[pid]) { - LOGF(warning, "Data fluc: Correction maps (Eff, Fake, or Flat) are null for species index %d", isp); - return; - } + if ((cfgEff || cfgFlat) && (!state.hEff || !state.hFake || !state.hFlatWeight)) { + LOGF(warning, "Data fluc: correction maps requested but not all present."); + return; } - double sumpmwk[KNsp][KNEta][KIntM][KIntK]{}; - double sumwk[KNsp][KNEta][KIntK]{}; - - double mean[KNsp][KNEta]{}, c2[KNsp][KNEta]{}; - double p1kBar[KNsp][KNEta]{}; - double meanMult[KNsp][KNEta]{}, p1kBarMult[KNsp][KNEta]{}; + std::array, KIntM>, KNEtaMax> sumpmwk{}; + std::array, KNEtaMax> sumwk{}; + std::array mean{}, c2{}, p1kBar{}; + std::array meanMult{}, p1kBarMult{}; float vz = coll.posZ(); for (const auto& track : tracks) { - if (!isTrackSelected(track)) + if (!isTrackSelected(track)) { continue; - + } float p = track.p(); float pt = track.pt(); - float eta = track.eta(); - float phi = track.phi(); + float eta = track.eta(), phi = track.phi(); auto sign = track.sign(); - - if (p < KFloatEpsilon) + if (p < KFloatEpsilon) { continue; - - if (pt <= cfgPtMin || pt > cfgPtMax) + } + if (pt <= cfgPtMin || pt > cfgPtMax) { continue; - int id = identifyTrack(track, cent); - bool isPi = (id == KPidPionOne); - bool isKa = (id == KPidKaonTwo); - bool isPr = (id == KPidProtonThree); - bool isSpecies[KNsp] = { - true, - isPi && sign < 0, isPi && sign > 0, isPi, - isKa && sign < 0, isKa && sign > 0, isKa, - isPr && sign < 0, isPr && sign > 0, isPr}; - - for (int isp = 0; isp < KNsp; ++isp) { - if (!isSpecies[isp]) - continue; - float eff = getEfficiency(coll.multNTracksPV(), pt, eta, static_cast(isp), 0, cfgEff); - if (eff <= KFloatEpsilon) - continue; + } - float fake = getEfficiency(coll.multNTracksPV(), pt, eta, static_cast(isp), 1, cfgEff); - float flatWeight = getFlatteningWeight(vz, sign, pt, eta, phi, static_cast(isp), cfgFlat); - float w = flatWeight * (1.0f - fake) / eff; + float eff = getEfficiency(coll.multNTracksPV(), pt, eta, 0, cfgEff); + if (eff <= KFloatEpsilon) { + continue; + } + float fake = getEfficiency(coll.multNTracksPV(), pt, eta, 1, cfgEff); + float flatWeight = getFlatteningWeight(vz, sign, pt, eta, phi, cfgFlat); + float w = flatWeight * (1.0f - fake) / eff; + if (!std::isfinite(w) || w <= 0.f) { + continue; + } - if (!std::isfinite(w) || w <= 0.f) + for (int ieta = 0; ieta < nEta; ++ieta) { + if (eta <= etaLw[ieta] || eta > etaUp[ieta]) { continue; - - for (int ieta = 0; ieta < KNEta; ++ieta) { - if (eta <= etaLw[ieta] || eta > etaUp[ieta]) - continue; - for (int k = 0; k < KIntK; ++k) { - for (int m = 0; m < KIntM; ++m) { - sumpmwk[isp][ieta][m][k] += std::pow(w, k) * std::pow(pt, m); - } - sumwk[isp][ieta][k] += std::pow(w, k); + } + for (int k = 0; k < KIntK; ++k) { + for (int m = 0; m < KIntM; ++m) { + sumpmwk[ieta][m][k] += std::pow(w, k) * std::pow(pt, m); } + sumwk[ieta][k] += std::pow(w, k); } } } - if (!hasMinTracksInAllEtaBins(sumwk)) + if (!hasMinTracksInAllEtaBins(sumwk)) { return; + } double amplFT0A = 0, amplFT0C = 0; if (coll.has_foundFT0()) { const auto& ft0 = coll.foundFT0(); for (std::size_t iCh = 0; iCh < ft0.channelA().size(); iCh++) { - float ampl = ft0.amplitudeA()[iCh]; - amplFT0A += ampl; + amplFT0A += ft0.amplitudeA()[iCh]; } for (std::size_t iCh = 0; iCh < ft0.channelC().size(); iCh++) { - float ampl = ft0.amplitudeC()[iCh]; - amplFT0C += ampl; + amplFT0C += ft0.amplitudeC()[iCh]; } } double p1kBarFt0A = amplFT0A - state.pmeanFT0AmultpvStep2->GetBinContent(state.pmeanFT0AmultpvStep2->GetXaxis()->FindBin(coll.multNTracksPV())); double p1kBarFt0C = amplFT0C - state.pmeanFT0CmultpvStep2->GetBinContent(state.pmeanFT0CmultpvStep2->GetXaxis()->FindBin(coll.multNTracksPV())); - for (int ieta = 0; ieta < KNEta; ++ieta) { - const int ibx = state.pmeanNchEtabinSpbinStep2->GetXaxis()->FindBin(coll.multNTracksPV()); + for (int ieta = 0; ieta < nEta; ++ieta) { + const int ibx = state.pmeanNchEtabinStep2->GetXaxis()->FindBin(coll.multNTracksPV()); const int iby = ieta + 1; - for (int isp = 0; isp < KNsp; ++isp) { - const int ibz = isp + 1; - - float mmpt = state.pmeanNchEtabinSpbinStep2->GetBinContent(ibx, iby, ibz); - float mmMult = state.pmeanMultNchEtabinSpbinStep2->GetBinContent(ibx, iby, ibz); + float mmpt = state.pmeanNchEtabinStep2->GetBinContent(ibx, iby); + float mmMult = state.pmeanMultNchEtabinStep2->GetBinContent(ibx, iby); - mean[isp][ieta] = sumpmwk[isp][ieta][1][1] / sumwk[isp][ieta][1]; - meanMult[isp][ieta] = sumwk[isp][ieta][1]; + mean[ieta] = sumpmwk[ieta][1][1] / sumwk[ieta][1]; + meanMult[ieta] = sumwk[ieta][1]; - if (std::isfinite(mmpt)) { - std::tie(mean[isp][ieta], c2[isp][ieta]) = calculateMeanAndC2FromSums(sumpmwk[isp][ieta], sumwk[isp][ieta], mmpt); - p1kBar[isp][ieta] = mean[isp][ieta] - mmpt; - } - p1kBarMult[isp][ieta] = meanMult[isp][ieta] - mmMult; + if (std::isfinite(mmpt) && mmpt != 0) { + std::tie(mean[ieta], c2[ieta]) = calculateMeanAndC2FromSums(sumpmwk[ieta], sumwk[ieta], mmpt); + p1kBar[ieta] = mean[ieta] - mmpt; } + p1kBarMult[ieta] = meanMult[ieta] - mmMult; } - for (int ieta = 0; ieta < KNEta; ++ieta) { - for (int isp = 0; isp < KNsp; ++isp) { - if (std::isfinite(mean[isp][ieta])) { - histos.fill(HIST("Prof_MeanpT_Cent_etabin_spbin"), cent, ieta, isp, mean[isp][ieta]); - histos.fill(HIST("Prof_MeanpT_Mult_etabin_spbin"), coll.multNTracksPV(), ieta, isp, mean[isp][ieta]); - } - if (std::isfinite(c2[isp][ieta])) { - histos.fill(HIST("Prof_C2_Cent_etabin_spbin"), cent, ieta, isp, c2[isp][ieta]); - histos.fill(HIST("Prof_C2_Mult_etabin_spbin"), coll.multNTracksPV(), ieta, isp, c2[isp][ieta]); - } + + // --- Poisson bootstrap: one weight per sample per event --- + std::array poisW{}; + if (doBoot) { + for (int s = 0; s < nBoot; ++s) { + poisW[s] = rng.Poisson(1.0); } } + auto fillBS2D = [&](std::array, KMaxBoot>& arr, double x, double y, double val) { + if (!doBoot) { + return; + } + for (int s = 0; s < nBoot; ++s) { + arr[s]->Fill(x, y, val, poisW[s]); + } + }; + auto fillBS3D = [&](std::array, KMaxBoot>& arr, double x, double y, double z, double val) { + if (!doBoot) { + return; + } + for (int s = 0; s < nBoot; ++s) { + arr[s]->Fill(x, y, z, val, poisW[s]); + } + }; - for (int ietaA = 1; ietaA <= (KNEta - 1) / 2; ++ietaA) { - int ietaC = KNEta - ietaA; - for (int isp = 0; isp < KNsp; ++isp) { - float c2Sub = p1kBar[isp][ietaA] * p1kBar[isp][ietaC]; - float covAC = p1kBarMult[isp][ietaA] * p1kBar[isp][ietaC]; - float covCA = p1kBar[isp][ietaA] * p1kBarMult[isp][ietaC]; + // meanpT & C2 vs eta bin + for (int ieta = 0; ieta < nEta; ++ieta) { + if (std::isfinite(mean[ieta])) { + histos.fill(HIST("Prof_MeanpT_Cent_etabin"), cent, ieta, mean[ieta]); + histos.fill(HIST("Prof_MeanpT_Mult_etabin"), coll.multNTracksPV(), ieta, mean[ieta]); + fillBS2D(bs.meanpTCent, cent, ieta, mean[ieta]); + fillBS2D(bs.meanpTMult, coll.multNTracksPV(), ieta, mean[ieta]); + } + if (std::isfinite(c2[ieta])) { + histos.fill(HIST("Prof_C2_Cent_etabin"), cent, ieta, c2[ieta]); + histos.fill(HIST("Prof_C2_Mult_etabin"), coll.multNTracksPV(), ieta, c2[ieta]); + fillBS2D(bs.c2Cent, cent, ieta, c2[ieta]); + fillBS2D(bs.c2Mult, coll.multNTracksPV(), ieta, c2[ieta]); + } + } - float covFT0A = p1kBarFt0A * p1kBar[isp][ietaC]; - float covFT0C = p1kBarFt0C * p1kBar[isp][ietaA]; + // mirror-pair subevent (C2Sub) & covariances vs eta bin + for (int ietaA = 1; ietaA <= (nEta - 1) / 2; ++ietaA) { + int ietaC = nEta - ietaA; + float c2Sub = p1kBar[ietaA] * p1kBar[ietaC]; + float covAC = p1kBarMult[ietaA] * p1kBar[ietaC]; + float covCA = p1kBar[ietaA] * p1kBarMult[ietaC]; - if (std::isfinite(c2Sub)) { - histos.fill(HIST("Prof_C2Sub_Cent_etabin_spbin"), cent, ietaA, isp, c2Sub); - histos.fill(HIST("Prof_C2Sub_Mult_etabin_spbin"), coll.multNTracksPV(), ietaA, isp, c2Sub); - } - if (std::isfinite(covAC)) { - histos.fill(HIST("Prof_Cov_Cent_etabin_spbin"), cent, ietaA, isp, covAC); - histos.fill(HIST("Prof_Cov_Mult_etabin_spbin"), coll.multNTracksPV(), ietaA, isp, covAC); - } - if (std::isfinite(covCA)) { - histos.fill(HIST("Prof_Cov_Cent_etabin_spbin"), cent, ietaC, isp, covCA); - histos.fill(HIST("Prof_Cov_Mult_etabin_spbin"), coll.multNTracksPV(), ietaC, isp, covCA); - } - if (std::isfinite(covFT0A)) { - histos.fill(HIST("Prof_CovFT0A_Cent_etabin_spbin"), cent, ietaA, isp, covFT0A); - histos.fill(HIST("Prof_CovFT0A_Mult_etabin_spbin"), coll.multNTracksPV(), ietaA, isp, covFT0A); - } - if (std::isfinite(covFT0C)) { - histos.fill(HIST("Prof_CovFT0C_Cent_etabin_spbin"), cent, ietaA, isp, covFT0C); - histos.fill(HIST("Prof_CovFT0C_Mult_etabin_spbin"), coll.multNTracksPV(), ietaA, isp, covFT0C); - } + if (std::isfinite(c2Sub)) { + histos.fill(HIST("Prof_C2Sub_Cent_etabin"), cent, ietaA, c2Sub); + histos.fill(HIST("Prof_C2Sub_Mult_etabin"), coll.multNTracksPV(), ietaA, c2Sub); + fillBS2D(bs.c2SubCent, cent, ietaA, c2Sub); + fillBS2D(bs.c2SubMult, coll.multNTracksPV(), ietaA, c2Sub); + } + if (std::isfinite(covAC)) { + histos.fill(HIST("Prof_Cov_Cent_etabin"), cent, ietaA, covAC); + histos.fill(HIST("Prof_Cov_Mult_etabin"), coll.multNTracksPV(), ietaA, covAC); + fillBS2D(bs.covCent, cent, ietaA, covAC); + fillBS2D(bs.covMult, coll.multNTracksPV(), ietaA, covAC); + } + if (std::isfinite(covCA)) { + histos.fill(HIST("Prof_Cov_Cent_etabin"), cent, ietaC, covCA); + histos.fill(HIST("Prof_Cov_Mult_etabin"), coll.multNTracksPV(), ietaC, covCA); + fillBS2D(bs.covCent, cent, ietaC, covCA); + fillBS2D(bs.covMult, coll.multNTracksPV(), ietaC, covCA); } } - for (int ietaA = 1; ietaA < KNEta; ++ietaA) { - for (int ietaC = 1; ietaC < KNEta; ++ietaC) { + // FT0 covariance vs narrow eta bin (full range, indexed by the actual pT bin) + for (int ieta = 1; ieta < nEta; ++ieta) { + float covFT0Aeta = p1kBarFt0A * p1kBar[ieta]; + float covFT0Ceta = p1kBarFt0C * p1kBar[ieta]; + if (std::isfinite(covFT0Aeta)) { + histos.fill(HIST("Prof_CovFT0A_Cent_etabin"), cent, ieta, covFT0Aeta); + histos.fill(HIST("Prof_CovFT0A_Mult_etabin"), coll.multNTracksPV(), ieta, covFT0Aeta); + fillBS2D(bs.covFT0ACent, cent, ieta, covFT0Aeta); + fillBS2D(bs.covFT0AMult, coll.multNTracksPV(), ieta, covFT0Aeta); + } + if (std::isfinite(covFT0Ceta)) { + histos.fill(HIST("Prof_CovFT0C_Cent_etabin"), cent, ieta, covFT0Ceta); + histos.fill(HIST("Prof_CovFT0C_Mult_etabin"), coll.multNTracksPV(), ieta, covFT0Ceta); + fillBS2D(bs.covFT0CCent, cent, ieta, covFT0Ceta); + fillBS2D(bs.covFT0CMult, coll.multNTracksPV(), ieta, covFT0Ceta); + } + } + // full 2D subevent map + for (int ietaA = 1; ietaA < nEta; ++ietaA) { + for (int ietaC = 1; ietaC < nEta; ++ietaC) { float etaValA = (etaLw[ietaA] + etaUp[ietaA]) / 2.0f; float etaValB = (etaLw[ietaC] + etaUp[ietaC]) / 2.0f; float gap = etaValA - etaValB; float sum = (etaValA + etaValB); - for (int isp = 0; isp < KNsp; ++isp) { - - float c2Sub = (ietaA == ietaC) ? static_cast(c2[isp][ietaA]) - : p1kBar[isp][ietaA] * p1kBar[isp][ietaC]; - float cov = p1kBarMult[isp][ietaA] * p1kBar[isp][ietaC]; - float covFT0A = p1kBarFt0A * p1kBar[isp][ietaC]; - float covFT0C = p1kBarFt0C * p1kBar[isp][ietaA]; + float c2Sub = (ietaA == ietaC) ? static_cast(c2[ietaA]) : p1kBar[ietaA] * p1kBar[ietaC]; + float cov = p1kBarMult[ietaA] * p1kBar[ietaC]; + float covFT0A = p1kBarFt0A * p1kBar[ietaC]; + float covFT0C = p1kBarFt0C * p1kBar[ietaA]; - if (isp == kInclusiveIdx) { - if (std::isfinite(c2Sub)) { - histos.fill(HIST("Prof_C2Sub2D_Cent_etaA_etaC"), cent, etaValA, etaValB, c2Sub); - histos.fill(HIST("Prof_GapSum2D"), cent, gap, sum, c2Sub); - } - if (std::isfinite(cov)) - histos.fill(HIST("Prof_Cov2D_Cent_etaA_etaC"), cent, etaValA, etaValB, cov); - if (std::isfinite(covFT0A)) - histos.fill(HIST("Prof_CovFT0A2D_Cent_etaA_etaC"), cent, etaValA, etaValB, covFT0A); - if (std::isfinite(covFT0C)) - histos.fill(HIST("Prof_CovFT0C2D_Cent_etaA_etaC"), cent, etaValA, etaValB, covFT0C); - } else if (isp == kPiMinusIdx) { - if (std::isfinite(c2Sub)) { - histos.fill(HIST("Prof_C2Sub2D_Cent_etaA_etaC_PiMinus"), cent, etaValA, etaValB, c2Sub); - histos.fill(HIST("Prof_GapSum2D_PiMinus"), cent, gap, sum, c2Sub); - } - if (std::isfinite(cov)) - histos.fill(HIST("Prof_Cov2D_Cent_etaA_etaC_PiMinus"), cent, etaValA, etaValB, cov); - if (std::isfinite(covFT0A)) - histos.fill(HIST("Prof_CovFT0A2D_Cent_etaA_etaC_PiMinus"), cent, etaValA, etaValB, covFT0A); - if (std::isfinite(covFT0C)) - histos.fill(HIST("Prof_CovFT0C2D_Cent_etaA_etaC_PiMinus"), cent, etaValA, etaValB, covFT0C); - } else if (isp == kPiPlusIdx) { - if (std::isfinite(c2Sub)) { - histos.fill(HIST("Prof_C2Sub2D_Cent_etaA_etaC_PiPlus"), cent, etaValA, etaValB, c2Sub); - histos.fill(HIST("Prof_GapSum2D_PiPlus"), cent, gap, sum, c2Sub); - } - if (std::isfinite(cov)) - histos.fill(HIST("Prof_Cov2D_Cent_etaA_etaC_PiPlus"), cent, etaValA, etaValB, cov); - if (std::isfinite(covFT0A)) - histos.fill(HIST("Prof_CovFT0A2D_Cent_etaA_etaC_PiPlus"), cent, etaValA, etaValB, covFT0A); - if (std::isfinite(covFT0C)) - histos.fill(HIST("Prof_CovFT0C2D_Cent_etaA_etaC_PiPlus"), cent, etaValA, etaValB, covFT0C); - } else if (isp == kPiAllIdx) { - if (std::isfinite(c2Sub)) { - histos.fill(HIST("Prof_C2Sub2D_Cent_etaA_etaC_PiAll"), cent, etaValA, etaValB, c2Sub); - histos.fill(HIST("Prof_GapSum2D_PiAll"), cent, gap, sum, c2Sub); - } - if (std::isfinite(cov)) - histos.fill(HIST("Prof_Cov2D_Cent_etaA_etaC_PiAll"), cent, etaValA, etaValB, cov); - if (std::isfinite(covFT0A)) - histos.fill(HIST("Prof_CovFT0A2D_Cent_etaA_etaC_PiAll"), cent, etaValA, etaValB, covFT0A); - if (std::isfinite(covFT0C)) - histos.fill(HIST("Prof_CovFT0C2D_Cent_etaA_etaC_PiAll"), cent, etaValA, etaValB, covFT0C); - } else if (isp == kKaMinusIdx) { - if (std::isfinite(c2Sub)) { - histos.fill(HIST("Prof_C2Sub2D_Cent_etaA_etaC_KaMinus"), cent, etaValA, etaValB, c2Sub); - histos.fill(HIST("Prof_GapSum2D_KaMinus"), cent, gap, sum, c2Sub); - } - if (std::isfinite(cov)) - histos.fill(HIST("Prof_Cov2D_Cent_etaA_etaC_KaMinus"), cent, etaValA, etaValB, cov); - if (std::isfinite(covFT0A)) - histos.fill(HIST("Prof_CovFT0A2D_Cent_etaA_etaC_KaMinus"), cent, etaValA, etaValB, covFT0A); - if (std::isfinite(covFT0C)) - histos.fill(HIST("Prof_CovFT0C2D_Cent_etaA_etaC_KaMinus"), cent, etaValA, etaValB, covFT0C); - } else if (isp == kKaPlusIdx) { - if (std::isfinite(c2Sub)) { - histos.fill(HIST("Prof_C2Sub2D_Cent_etaA_etaC_KaPlus"), cent, etaValA, etaValB, c2Sub); - histos.fill(HIST("Prof_GapSum2D_KaPlus"), cent, gap, sum, c2Sub); - } - if (std::isfinite(cov)) - histos.fill(HIST("Prof_Cov2D_Cent_etaA_etaC_KaPlus"), cent, etaValA, etaValB, cov); - if (std::isfinite(covFT0A)) - histos.fill(HIST("Prof_CovFT0A2D_Cent_etaA_etaC_KaPlus"), cent, etaValA, etaValB, covFT0A); - if (std::isfinite(covFT0C)) - histos.fill(HIST("Prof_CovFT0C2D_Cent_etaA_etaC_KaPlus"), cent, etaValA, etaValB, covFT0C); - } else if (isp == kKaAllIdx) { - if (std::isfinite(c2Sub)) { - histos.fill(HIST("Prof_C2Sub2D_Cent_etaA_etaC_KaAll"), cent, etaValA, etaValB, c2Sub); - histos.fill(HIST("Prof_GapSum2D_KaAll"), cent, gap, sum, c2Sub); - } - if (std::isfinite(cov)) - histos.fill(HIST("Prof_Cov2D_Cent_etaA_etaC_KaAll"), cent, etaValA, etaValB, cov); - if (std::isfinite(covFT0A)) - histos.fill(HIST("Prof_CovFT0A2D_Cent_etaA_etaC_KaAll"), cent, etaValA, etaValB, covFT0A); - if (std::isfinite(covFT0C)) - histos.fill(HIST("Prof_CovFT0C2D_Cent_etaA_etaC_KaAll"), cent, etaValA, etaValB, covFT0C); - } else if (isp == kPrIdx) { - if (std::isfinite(c2Sub)) { - histos.fill(HIST("Prof_C2Sub2D_Cent_etaA_etaC_Pr"), cent, etaValA, etaValB, c2Sub); - histos.fill(HIST("Prof_GapSum2D_Pr"), cent, gap, sum, c2Sub); - } - if (std::isfinite(cov)) - histos.fill(HIST("Prof_Cov2D_Cent_etaA_etaC_Pr"), cent, etaValA, etaValB, cov); - if (std::isfinite(covFT0A)) - histos.fill(HIST("Prof_CovFT0A2D_Cent_etaA_etaC_Pr"), cent, etaValA, etaValB, covFT0A); - if (std::isfinite(covFT0C)) - histos.fill(HIST("Prof_CovFT0C2D_Cent_etaA_etaC_Pr"), cent, etaValA, etaValB, covFT0C); - } else if (isp == kAntiPrIdx) { - if (std::isfinite(c2Sub)) { - histos.fill(HIST("Prof_C2Sub2D_Cent_etaA_etaC_AntiPr"), cent, etaValA, etaValB, c2Sub); - histos.fill(HIST("Prof_GapSum2D_AntiPr"), cent, gap, sum, c2Sub); - } - if (std::isfinite(cov)) - histos.fill(HIST("Prof_Cov2D_Cent_etaA_etaC_AntiPr"), cent, etaValA, etaValB, cov); - if (std::isfinite(covFT0A)) - histos.fill(HIST("Prof_CovFT0A2D_Cent_etaA_etaC_AntiPr"), cent, etaValA, etaValB, covFT0A); - if (std::isfinite(covFT0C)) - histos.fill(HIST("Prof_CovFT0C2D_Cent_etaA_etaC_AntiPr"), cent, etaValA, etaValB, covFT0C); - } else if (isp == kPrAllIdx) { - if (std::isfinite(c2Sub)) { - histos.fill(HIST("Prof_C2Sub2D_Cent_etaA_etaC_PrAll"), cent, etaValA, etaValB, c2Sub); - histos.fill(HIST("Prof_GapSum2D_PrAll"), cent, gap, sum, c2Sub); - } - if (std::isfinite(cov)) - histos.fill(HIST("Prof_Cov2D_Cent_etaA_etaC_PrAll"), cent, etaValA, etaValB, cov); - if (std::isfinite(covFT0A)) - histos.fill(HIST("Prof_CovFT0A2D_Cent_etaA_etaC_PrAll"), cent, etaValA, etaValB, covFT0A); - if (std::isfinite(covFT0C)) - histos.fill(HIST("Prof_CovFT0C2D_Cent_etaA_etaC_PrAll"), cent, etaValA, etaValB, covFT0C); - } + if (std::isfinite(c2Sub)) { + histos.fill(HIST("Prof_C2Sub2D_Cent_etaA_etaC"), cent, etaValA, etaValB, c2Sub); + histos.fill(HIST("Prof_GapSum2D"), cent, gap, sum, c2Sub); + fillBS3D(bs.c2Sub2D, cent, etaValA, etaValB, c2Sub); + fillBS3D(bs.gapSum2D, cent, gap, sum, c2Sub); + } + if (std::isfinite(cov)) { + histos.fill(HIST("Prof_Cov2D_Cent_etaA_etaC"), cent, etaValA, etaValB, cov); + fillBS3D(bs.cov2D, cent, etaValA, etaValB, cov); + } + if (std::isfinite(covFT0A)) { + histos.fill(HIST("Prof_CovFT0A2D_Cent_etaA_etaC"), cent, etaValA, etaValB, covFT0A); + fillBS3D(bs.covFT0A2D, cent, etaValA, etaValB, covFT0A); + } + if (std::isfinite(covFT0C)) { + histos.fill(HIST("Prof_CovFT0C2D_Cent_etaA_etaC"), cent, etaValA, etaValB, covFT0C); + fillBS3D(bs.covFT0C2D, cent, etaValA, etaValB, covFT0C); } } }