From 1d06bca61a860088d1922c8d2b967ccafe05f24d Mon Sep 17 00:00:00 2001 From: Somadutta Bhatta Date: Fri, 13 Mar 2026 17:16:52 +0100 Subject: [PATCH 1/3] [PWGCF] update track selection --- .../Tasks/radialFlowDecorr.cxx | 735 +++++++++--------- 1 file changed, 346 insertions(+), 389 deletions(-) diff --git a/PWGCF/EbyEFluctuations/Tasks/radialFlowDecorr.cxx b/PWGCF/EbyEFluctuations/Tasks/radialFlowDecorr.cxx index 33b16fb7c06..bc7c5078594 100644 --- a/PWGCF/EbyEFluctuations/Tasks/radialFlowDecorr.cxx +++ b/PWGCF/EbyEFluctuations/Tasks/radialFlowDecorr.cxx @@ -130,7 +130,7 @@ struct RadialFlowDecorr { KNsp }; - const std::vector pidSuffix = {"", "_PiMinus", "_PiPlus", "_PiAll", "_KaMinus", "_KaPlus", "_KaAll", "_AntiPr", "_Pr", "_PrAll"}; + inline static const std::vector pidSuffix = {"", "_PiMinus", "_PiPlus", "_PiAll", "_KaMinus", "_KaPlus", "_KaAll", "_AntiPr", "_Pr", "_PrAll"}; struct PIDMeanSigmaMap { static constexpr int MaxCentBins = 100; @@ -155,10 +155,10 @@ struct RadialFlowDecorr { kpp = 4 }; static constexpr float KinvalidCentrality = -1.0f; - const std::vector etaLw = { + inline static const std::vector etaLw = { -0.8, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7}; - const std::vector etaUp = { + inline static const std::vector etaUp = { 0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8}; @@ -166,10 +166,6 @@ struct RadialFlowDecorr { Configurable cfgPtMin{"cfgPtMin", 0.2f, "min pT"}; Configurable cfgPtMax{"cfgPtMax", 5.0f, "max pT"}; Configurable cfgEtaCut{"cfgEtaCut", 0.8f, "|η| cut"}; - Configurable cfgTPCClsMin{"cfgTPCClsMin", 70.f, "min TPC clusters"}; - Configurable cfgChi2TPCMax{"cfgChi2TPCMax", 4.0f, "max TPC χ²"}; - Configurable cfgCutTpcChi2NCl{"cfgCutTpcChi2NCl", 2.5f, "Maximum TPCchi2NCl"}; - Configurable cfgCutItsChi2NCl{"cfgCutItsChi2NCl", 36.0f, "Maximum ITSchi2NCl"}; Configurable cfgCutVertex{"cfgCutVertex", 10.0f, "Accepted z-vertex range"}; Configurable cfgCutTracKDcaMaxZ{"cfgCutTracKDcaMaxZ", 2.0f, "Maximum DcaZ"}; Configurable cfgCutTracKDcaMaxXY{"cfgCutTracKDcaMaxXY", 0.2f, "Maximum DcaZ"}; @@ -193,22 +189,26 @@ struct RadialFlowDecorr { Configurable cfgnSigmaCutTOF{"cfgnSigmaCutTOF", 2.0f, "PID nSigma cut for TOF"}; Configurable cfgnSigmaCutCombTPCTOF{"cfgnSigmaCutCombTPCTOF", 2.0f, "PID nSigma combined cut for TPC and TOF"}; - Configurable cfgTpcElRejCutMin{"cfgTpcElRejCutMin", -3.0f, "Electron Rejection Cut Minimum"}; - Configurable cfgTpcElRejCutMax{"cfgTpcElRejCutMax", 5.0f, "Electron Rejection Cut Maximum"}; - Configurable cfgTpcElRejCut{"cfgTpcElRejCut", 3.0f, "TPC Hadron Rejection Cut"}; - 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 cfgCutEta{"cfgCutEta", 0.8f, "absolute Eta cut"}; Configurable cfgNsubsample{"cfgNsubsample", 10, "Number of subsamples"}; - Configurable cfgCentralityChoice{"cfgCentralityChoice", 2, "Which centrality estimator? 1-->FT0C, 2-->FT0M, 3-->FDDM, 4-->FV0A"}; + 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"}; Configurable cfgIsGoodZvtxFT0VsPV{"cfgIsGoodZvtxFT0VsPV", true, "Good Vertexing cut"}; - Configurable cfgNchPbMax{"cfgNchPbMax", 6000, "Max Nch range for PbPb collisions"}; - Configurable cfgNchOMax{"cfgNchOMax", 600, "Max Nch range for OO collisions"}; + Configurable cfgPupnSig{"cfgPupnSig", 3.0f, "Additional Pileup Cut"}; + Configurable cfgApplySigPupCut{"cfgApplySigPupCut", 0, "nSig Pileup Cut"}; + Configurable cfgApplyLinPupCut{"cfgApplyLinPupCut", 0, "Lin Pileup Cut"}; + Configurable cfgLinPupParam0{"cfgLinPupParam0", 3.0f, "(Upper) Linear Pileup Cut Const"}; + Configurable cfgLinPupParam1{"cfgLinPupParam1", 3.0f, "(Upper) Linear Pileup Slope"}; + Configurable cfgLinPupParam2{"cfgLinPupParam2", 3.0f, "(Lower) Linear Pileup Cut Const"}; + Configurable cfgLinPupParam3{"cfgLinPupParam3", 3.0f, "(Lower) Linear Pileup Slope"}; + + Configurable cfgNchPbMax{"cfgNchPbMax", 10000, "Max Nch range for PbPb collisions"}; + Configurable cfgNchOMax{"cfgNchOMax", 1000, "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 cfgFlat{"cfgFlat", false, "Whether to use flattening weights"}; @@ -266,6 +266,9 @@ struct RadialFlowDecorr { std::array hFake{}; std::array hFlatWeight{}; + std::vector> mLimitsNchCent; + float mMinXNchCent = 0, mMaxXNchCent = 0; + TProfile3D* pmeanTruNchEtabinSpbinStep2 = nullptr; TProfile3D* pmeanRecoNchEtabinSpbinStep2 = nullptr; TProfile3D* pmeanRecoEffcorrNchEtabinSpbinStep2 = nullptr; @@ -287,25 +290,16 @@ struct RadialFlowDecorr { { float pt = track.pt(); auto sign = track.sign(); - // --- Generic Inclusive --- - histos.fill(HIST("h3DnsigmaTpcVsPtBefCut_Cent"), cent, pt, track.tpcNSigmaPi()); - histos.fill(HIST("h3DnsigmaTofVsPtBefCut_Cent"), cent, pt, track.tofNSigmaPi()); - histos.fill(HIST("h3DnsigmaTpcVsTofBefCut_Cent"), cent, track.tofNSigmaPi(), track.tpcNSigmaPi()); - - histos.fill(HIST("h3DnsigmaTpcVsPtBefCut_Cent"), cent, pt, track.tpcNSigmaKa()); - histos.fill(HIST("h3DnsigmaTofVsPtBefCut_Cent"), cent, pt, track.tofNSigmaKa()); - histos.fill(HIST("h3DnsigmaTpcVsTofBefCut_Cent"), cent, track.tofNSigmaKa(), track.tpcNSigmaKa()); - - histos.fill(HIST("h3DnsigmaTpcVsPtBefCut_Cent"), cent, pt, track.tpcNSigmaPr()); - histos.fill(HIST("h3DnsigmaTofVsPtBefCut_Cent"), cent, pt, track.tofNSigmaPr()); - histos.fill(HIST("h3DnsigmaTpcVsTofBefCut_Cent"), cent, track.tofNSigmaPr(), track.tpcNSigmaPr()); + 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()); @@ -313,9 +307,11 @@ struct RadialFlowDecorr { 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()); @@ -337,99 +333,122 @@ struct RadialFlowDecorr { void fillNSigmaAftCut(const T& track, float cent, bool isSpecies[]) { float pt = track.pt(); + float tpcPi = track.tpcNSigmaPi(); + float tofPi = track.tofNSigmaPi(); - if (isSpecies[kInclusiveIdx]) { - histos.fill(HIST("h3DnsigmaTpcVsPtAftCut_Cent"), cent, pt, track.tpcNSigmaPi()); - histos.fill(HIST("h3DnsigmaTofVsPtAftCut_Cent"), cent, pt, track.tofNSigmaPi()); - histos.fill(HIST("h3DnsigmaTpcVsTofAftCut_Cent"), cent, track.tofNSigmaPi(), track.tpcNSigmaPi()); + float tpcKa = track.tpcNSigmaKa(); + float tofKa = track.tofNSigmaKa(); - histos.fill(HIST("h3DnsigmaTpcVsPtAftCut_Cent"), cent, pt, track.tpcNSigmaKa()); - histos.fill(HIST("h3DnsigmaTofVsPtAftCut_Cent"), cent, pt, track.tofNSigmaKa()); - histos.fill(HIST("h3DnsigmaTpcVsTofAftCut_Cent"), cent, track.tofNSigmaKa(), track.tpcNSigmaKa()); + float tpcPr = track.tpcNSigmaPr(); + float tofPr = track.tofNSigmaPr(); - histos.fill(HIST("h3DnsigmaTpcVsPtAftCut_Cent"), cent, pt, track.tpcNSigmaPr()); - histos.fill(HIST("h3DnsigmaTofVsPtAftCut_Cent"), cent, pt, track.tofNSigmaPr()); - histos.fill(HIST("h3DnsigmaTpcVsTofAftCut_Cent"), cent, track.tofNSigmaPr(), track.tpcNSigmaPr()); - } if (isSpecies[kPiPlusIdx]) { - histos.fill(HIST("h3DnsigmaTpcVsPtAftCut_Cent_PiPlus"), cent, pt, track.tpcNSigmaPi()); - histos.fill(HIST("h3DnsigmaTofVsPtAftCut_Cent_PiPlus"), cent, pt, track.tofNSigmaPi()); - histos.fill(HIST("h3DnsigmaTpcVsTofAftCut_Cent_PiPlus"), cent, track.tofNSigmaPi(), track.tpcNSigmaPi()); - } else if (isSpecies[kPiMinusIdx]) { - histos.fill(HIST("h3DnsigmaTpcVsPtAftCut_Cent_PiMinus"), cent, pt, track.tpcNSigmaPi()); - histos.fill(HIST("h3DnsigmaTofVsPtAftCut_Cent_PiMinus"), cent, pt, track.tofNSigmaPi()); - histos.fill(HIST("h3DnsigmaTpcVsTofAftCut_Cent_PiMinus"), cent, track.tofNSigmaPi(), track.tpcNSigmaPi()); + 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, track.tpcNSigmaPi()); - histos.fill(HIST("h3DnsigmaTofVsPtAftCut_Cent_PiAll"), cent, pt, track.tofNSigmaPi()); - histos.fill(HIST("h3DnsigmaTpcVsTofAftCut_Cent_PiAll"), cent, track.tofNSigmaPi(), track.tpcNSigmaPi()); + 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, track.tpcNSigmaKa()); - histos.fill(HIST("h3DnsigmaTofVsPtAftCut_Cent_KaPlus"), cent, pt, track.tofNSigmaKa()); - histos.fill(HIST("h3DnsigmaTpcVsTofAftCut_Cent_KaPlus"), cent, track.tofNSigmaKa(), track.tpcNSigmaKa()); - } else if (isSpecies[kKaMinusIdx]) { - histos.fill(HIST("h3DnsigmaTpcVsPtAftCut_Cent_KaMinus"), cent, pt, track.tpcNSigmaKa()); - histos.fill(HIST("h3DnsigmaTofVsPtAftCut_Cent_KaMinus"), cent, pt, track.tofNSigmaKa()); - histos.fill(HIST("h3DnsigmaTpcVsTofAftCut_Cent_KaMinus"), cent, track.tofNSigmaKa(), track.tpcNSigmaKa()); + 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, track.tpcNSigmaKa()); - histos.fill(HIST("h3DnsigmaTofVsPtAftCut_Cent_KaAll"), cent, pt, track.tofNSigmaKa()); - histos.fill(HIST("h3DnsigmaTpcVsTofAftCut_Cent_KaAll"), cent, track.tofNSigmaKa(), track.tpcNSigmaKa()); + 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, track.tpcNSigmaPr()); - histos.fill(HIST("h3DnsigmaTofVsPtAftCut_Cent_Pr"), cent, pt, track.tofNSigmaPr()); - histos.fill(HIST("h3DnsigmaTpcVsTofAftCut_Cent_Pr"), cent, track.tofNSigmaPr(), track.tpcNSigmaPr()); - } else if (isSpecies[kAntiPrIdx]) { - histos.fill(HIST("h3DnsigmaTpcVsPtAftCut_Cent_AntiPr"), cent, pt, track.tpcNSigmaPr()); - histos.fill(HIST("h3DnsigmaTofVsPtAftCut_Cent_AntiPr"), cent, pt, track.tofNSigmaPr()); - histos.fill(HIST("h3DnsigmaTpcVsTofAftCut_Cent_AntiPr"), cent, track.tofNSigmaPr(), track.tpcNSigmaPr()); + 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, track.tpcNSigmaPr()); - histos.fill(HIST("h3DnsigmaTofVsPtAftCut_Cent_PrAll"), cent, pt, track.tofNSigmaPr()); - histos.fill(HIST("h3DnsigmaTpcVsTofAftCut_Cent_PrAll"), cent, track.tofNSigmaPr(), track.tpcNSigmaPr()); + 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 centBin) + int identifyTrack(const T& candidate, int cent) { - // Initial sanity checks if (!candidate.hasTPC()) return 0; + float pt = candidate.pt(); if (pt <= cfgCutPtLower || pt >= cfgCutPtUpperPID) return 0; // Out of bounds - // Determine species indices based on charge + + 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; - // Fetch Calibration Data (Means and Sigmas) + // TPC - float mPiTpc = pidMeanSigmaMap ? pidMeanSigmaMap->meanTPC[piIdx][centBin] : 0.f; - float sPiTpc = pidMeanSigmaMap ? pidMeanSigmaMap->sigmaTPC[piIdx][centBin] : 1.f; - float mKaTpc = pidMeanSigmaMap ? pidMeanSigmaMap->meanTPC[kaIdx][centBin] : 0.f; - float sKaTpc = pidMeanSigmaMap ? pidMeanSigmaMap->sigmaTPC[kaIdx][centBin] : 1.f; - float mPrTpc = pidMeanSigmaMap ? pidMeanSigmaMap->meanTPC[prIdx][centBin] : 0.f; - float sPrTpc = pidMeanSigmaMap ? pidMeanSigmaMap->sigmaTPC[prIdx][centBin] : 1.f; + 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 ? pidMeanSigmaMap->meanTOF[piIdx][centBin] : 0.f; - float sPiTof = pidMeanSigmaMap ? pidMeanSigmaMap->sigmaTOF[piIdx][centBin] : 1.f; - float mKaTof = pidMeanSigmaMap ? pidMeanSigmaMap->meanTOF[kaIdx][centBin] : 0.f; - float sKaTof = pidMeanSigmaMap ? pidMeanSigmaMap->sigmaTOF[kaIdx][centBin] : 1.f; - float mPrTof = pidMeanSigmaMap ? pidMeanSigmaMap->meanTOF[prIdx][centBin] : 0.f; - float sPrTof = pidMeanSigmaMap ? pidMeanSigmaMap->sigmaTOF[prIdx][centBin] : 1.f; + 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 rawTpcEl = candidate.tpcNSigmaEl(); float rawTofPi = 0.f, rawTofKa = 0.f, rawTofPr = 0.f; if (candidate.hasTOF()) { @@ -437,13 +456,6 @@ struct RadialFlowDecorr { rawTofKa = candidate.tofNSigmaKa(); rawTofPr = candidate.tofNSigmaPr(); } - // ELECTRON REJECTION: Reject if it is > cTpcRejCut from ALL hadron hypotheses AND falls within the electron band [-3, 5] - if (std::abs(rawTpcPi - mPiTpc) > cfgTpcElRejCut && - std::abs(rawTpcKa - mKaTpc) > cfgTpcElRejCut && - std::abs(rawTpcPr - mPrTpc) > cfgTpcElRejCut && - rawTpcEl > cfgTpcElRejCutMin && rawTpcEl < cfgTpcElRejCutMax) { - return 0; // It's an electron, reject it! - } // --- Low PT Regime --- if (pt <= cfgCutPtUpperTPC) { @@ -451,10 +463,12 @@ struct RadialFlowDecorr { 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; @@ -465,16 +479,19 @@ struct RadialFlowDecorr { 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) @@ -512,7 +529,32 @@ struct RadialFlowDecorr { 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) + { + auto checkLimits = [](float x, float y, const std::vector>& limits, float xM, float xMx) { + 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)) + return false; + return (y >= limits[bin].first && y <= limits[bin].second); + }; + if (cfgApplySigPupCut) { + if (!checkLimits(trksize, cent, state.mLimitsNchCent, state.mMinXNchCent, state.mMaxXNchCent)) + return false; + } + histos.fill(HIST("hEvtCount"), 6.5); + if (cfgApplyLinPupCut) { + if (trksize > (cfgLinPupParam0 + cfgLinPupParam1 * multPV)) + return false; + histos.fill(HIST("hEvtCount"), 7.5); + if (trksize < (cfgLinPupParam2 + cfgLinPupParam3 * multPV)) + return false; + histos.fill(HIST("hEvtCount"), 8.5); + } return true; } @@ -681,6 +723,10 @@ struct RadialFlowDecorr { return; } offsetFT0 = ccdb->getForTimeStamp>("FT0/Calib/Align", timestamp); + if (!offsetFT0) { + LOGF(fatal, "Failed to load valid FT0 alignment from CCDB!"); + 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()); @@ -754,37 +800,34 @@ struct RadialFlowDecorr { aod::CentFT0Cs, aod::CentFT0Ms, aod::CentFDDMs, aod::CentFV0As, aod::CentNGlobals, aod::McCollisionLabels>; - using MyMCTracks = soa::Join< - aod::Tracks, aod::TrackSelection, aod::TracksExtra, aod::TracksDCA, - aod::McTrackLabels, - aod::pidTPCFullPi, aod::pidTPCFullKa, aod::pidTPCFullPr, aod::pidTPCFullEl, - aod::pidTOFFullPi, aod::pidTOFFullKa, aod::pidTOFFullPr, aod::pidTOFFullEl>; - - PresliceUnsorted partPerMcCollision = aod::mcparticle::mcCollisionId; - PresliceUnsorted colPerMcCollision = aod::mccollisionlabel::mcCollisionId; - PresliceUnsorted trackPerMcParticle = aod::mctracklabel::mcParticleId; - Preslice perCollision = aod::track::collisionId; - Preslice trackPerCollision = aod::track::collisionId; + PresliceUnsorted colPerMcCollision = aod::mccollisionlabel::mcCollisionId; void declareCommonQA() { histos.add("hVtxZ_after_sel", ";z_{vtx} (cm)", kTH1F, {{KNbinsZvtx, KZvtxMin, KZvtxMax}}); histos.add("hVtxZ", ";z_{vtx} (cm)", kTH1F, {{KNbinsZvtx, KZvtxMin, KZvtxMax}}); histos.add("hCentrality", ";centrality (%)", kTH1F, {{centAxis1Per}}); - histos.add("Hist2D_globalTracks_PVTracks", ";N_{global};N_{PV}", kTH2F, {{nChAxis2}, {nChAxis2}}); - histos.add("Hist2D_cent_nch", ";N_{PV};cent (%)", kTH2F, {{nChAxis2}, {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}}); + histos.add("hP", ";p (GeV/c)", kTH1F, {{KNbinsPt, KPMin, KPMax}}); histos.add("hPt", ";p_{T} (GeV/c)", kTH1F, {{KNbinsPt, KPtMin, KPtMax}}); histos.add("hEta", ";#eta", kTH1F, {{KNbinsEtaFine, KEtaMin, KEtaMax}}); histos.add("hPhi", ";#phi", kTH1F, {{KNbinsPhi, KPhiMin, TwoPI}}); - histos.add("hEvtCount", "Number of Event;; Count", kTH1F, {{7, 0, 7}}); + histos.add("hEvtCount", "Number of Event;; Count", kTH1F, {{9, 0, 9}}); histos.get(HIST("hEvtCount"))->GetXaxis()->SetBinLabel(1, "all Events"); histos.get(HIST("hEvtCount"))->GetXaxis()->SetBinLabel(2, "after sel8"); histos.get(HIST("hEvtCount"))->GetXaxis()->SetBinLabel(3, "after VertexZ Cut"); histos.get(HIST("hEvtCount"))->GetXaxis()->SetBinLabel(4, "after kNoSameBunchPileup"); histos.get(HIST("hEvtCount"))->GetXaxis()->SetBinLabel(5, "after kIsGoodZvtxFT0vsPV"); histos.get(HIST("hEvtCount"))->GetXaxis()->SetBinLabel(6, "after kIsGoodITSLayersAll"); + histos.get(HIST("hEvtCount"))->GetXaxis()->SetBinLabel(7, "after PVTracksCent Pileup Cut"); + histos.get(HIST("hEvtCount"))->GetXaxis()->SetBinLabel(8, "after Linear Pileup Cut (Up)"); + histos.get(HIST("hEvtCount"))->GetXaxis()->SetBinLabel(9, "after Linear Pileup Cut (Lw)"); histos.add("hTrkCount", "Number of Tracks;; Count", kTH1F, {{11, 0, 11}}); histos.get(HIST("hTrkCount"))->GetXaxis()->SetBinLabel(1, "all Tracks"); @@ -941,8 +984,8 @@ struct RadialFlowDecorr { histos.add("pmeanFT0Cmultpv", "N_{PV}; AmplitudeA", kTProfile, {nChAxis}); histos.add("pmeanFT0C_cent", "cent; AmplitudeA", kTProfile, {centAxis1Per}); - histos.add("pmean_cent_id_eta_FT0", ";cent;channel id; #eta;amplitude", kTProfile3D, {{centAxis1Per}, {100, -0.5, 99.5}, {100, -5.0, 5.0}}); - histos.add("h3_cent_id_eta_FT0", ";cent;channel id; #eta", kTH3F, {{centAxis1Per}, {100, -0.5, 99.5}, {100, -5.0, 5.0}}); + 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}}); @@ -1180,7 +1223,9 @@ struct RadialFlowDecorr { if (hNumS && hNumF && hDenF) { state.hFake[pidType] = reinterpret_cast(hNumS->Clone(Form("hFake%s", suffix.c_str()))); state.hFake[pidType]->Add(hNumF); - state.hFake[pidType]->Add(hNumF2); + if (pidType != kInclusiveIdx && hNumF2) { + state.hFake[pidType]->Add(hNumF2); + } state.hFake[pidType]->SetDirectory(nullptr); state.hFake[pidType]->Divide(hDenF); } else { @@ -1193,50 +1238,48 @@ struct RadialFlowDecorr { } } - if (cfgRunGetEff || cfgRunGetMCFlat || cfgRunMCMean || cfgRunMCFluc) { - TList* pidList = ccdb->getForTimeStamp(pathNsig, now); + 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.", pathNsig.c_str()); + LOGF(warn, "nSigma maps required but CCDB list is null at %s! Using raw values.", currentPath.c_str()); } else { - if (!pidMeanSigmaMap) + 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); - // Initial Parameters: Amp, MeanTOF, SigmaTOF, MeanTPC, SigmaTPC - f2.SetParameter(0, h2->GetMaximum()); - f2.SetParameter(1, 0.0); - f2.SetParLimits(1, -2, 2); - - h2->Fit(&f2, "QRN"); + 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, " Bin %d: Mean TOF = %.3f, Mean TPC = %.3f", - iCent - 1, f2.GetParameter(1), f2.GetParameter(3)); + 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)); } } } @@ -1244,60 +1287,23 @@ struct RadialFlowDecorr { for (int i = 0; i < KNsp; ++i) { loadPIDMeans(static_cast(i)); } - } - } - - if (cfgRunGetDataFlat || cfgRunDataMean || cfgRunDataFluc) { - TList* pidList = ccdb->getForTimeStamp(pathDataNsig, now); - - if (!pidList) { - LOGF(warn, "nSigma maps required but CCDB list is null at %s! Using raw values.", pathDataNsig.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) { - TF2 f2("f2", "[0]*TMath::Gaus(x,[1],[2])*TMath::Gaus(y,[3],[4])", -3, 3, -3, 3); - // Initial Parameters: Amp, MeanTOF, SigmaTOF, MeanTPC, SigmaTPC - f2.SetParameter(0, h2->GetMaximum()); - f2.SetParameter(1, 0.0); - f2.SetParLimits(1, -2, 2); - - h2->Fit(&f2, "QRN"); - 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, " Bin %d: Mean TOF = %.3f, Mean TPC = %.3f", - iCent - 1, f2.GetParameter(1), f2.GetParameter(3)); - } - } + auto loadLimits = [&](const char* name, std::vector>& limits, float& xMin, float& xMax) { + auto* h2 = reinterpret_cast(pidList->FindObject(name)); + if (!h2) + return; // Skip if missing + int nBins = h2->GetXaxis()->GetNbins(); + xMin = h2->GetXaxis()->GetXmin(); + xMax = h2->GetXaxis()->GetXmax(); + limits.assign(nBins + 2, {-9999.f, 99999.f}); + for (int i = 1; i <= nBins; ++i) { + std::unique_ptr proj(h2->ProjectionY("_py", i, i)); + float m = proj->GetMean(); + float s = proj->GetRMS(); + limits[i] = {m - cfgPupnSig * s, m + cfgPupnSig * s}; } }; - for (int i = 0; i < KNsp; ++i) { - loadPIDMeans(static_cast(i)); - } + loadLimits("Hist2D_cent_nch", state.mLimitsNchCent, state.mMinXNchCent, state.mMaxXNchCent); } } @@ -1428,69 +1434,52 @@ struct RadialFlowDecorr { LOGF(info, "CCDB initialization complete for RadialFlowDecorr."); } - void processMCGetMeanNsig(aod::McCollisions const& mcColl, MyRun3MCCollisions const& collisions, TCs const& tracks, FilteredTCs const& /*filteredTracks*/, aod::McParticles const& mcParticles) + void processMCGetMeanNsig(MyRun3MCCollisions::iterator const& mcCollision, FilteredTCs const& mcTracks) { - for (const auto& mcCollision : mcColl) { - auto colSlice = collisions.sliceBy(colPerMcCollision, mcCollision.globalIndex()); - if (colSlice.size() != 1) - continue; - for (const auto& col : colSlice) { - histos.fill(HIST("hVtxZ"), col.posZ()); - if (!col.has_mcCollision() || !isEventSelected(col)) - continue; - auto trackSlice = tracks.sliceBy(trackPerCollision, col.globalIndex()); - auto partSlice = mcParticles.sliceBy(partPerMcCollision, mcCollision.globalIndex()); - if (trackSlice.size() < 1 || partSlice.size() < 1) - continue; - float cent = getCentrality(col); + histos.fill(HIST("hVtxZ"), mcCollision.posZ()); + if (!mcCollision.has_mcCollision() || !isEventSelected(mcCollision)) + return; + float cent = getCentrality(mcCollision); if (cent > KCentMax) - continue; - float multPV = col.multNTracksPV(); - for (const auto& particle : partSlice) { - if (!isParticleSelected(particle) || !particle.isPhysicalPrimary()) - continue; - histos.fill(HIST("hVtxZ_after_sel"), col.posZ()); - histos.fill(HIST("hCentrality"), cent); + return; + float multPV = mcCollision.multNTracksPV(); - histos.fill(HIST("Hist2D_globalTracks_PVTracks"), multPV, trackSlice.size()); - histos.fill(HIST("Hist2D_cent_nch"), trackSlice.size(), cent); + 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& track : trackSlice) { - if (!isTrackSelected(track)) - continue; - fillNSigmaBefCut(track, cent); - } + for (const auto& track : mcTracks) { + if (!isTrackSelected(track)) + continue; + fillNSigmaBefCut(track, cent); } } - } - } PROCESS_SWITCH(RadialFlowDecorr, processMCGetMeanNsig, "process MC to calculate Mean values of nSig Plots", cfgRunMCGetNSig); - void processGetEffHists(aod::McCollisions const& mcColl, MyRun3MCCollisions const& collisions, TCs const& tracks, FilteredTCs const& /*filteredTracks*/, aod::McParticles const& mcParticles) + void processGetEffHists(MyRun3MCCollisions::iterator const& mcCollision, FilteredTCs const& mcTracks, aod::McParticles const& mcParticles) { - for (const auto& mcCollision : mcColl) { - auto colSlice = collisions.sliceBy(colPerMcCollision, mcCollision.globalIndex()); - if (colSlice.size() != 1) - continue; - - for (const auto& col : colSlice) { - histos.fill(HIST("hVtxZ"), col.posZ()); - if (!col.has_mcCollision() || !isEventSelected(col)) - continue; - - auto trackSlice = tracks.sliceBy(trackPerCollision, col.globalIndex()); - auto partSlice = mcParticles.sliceBy(partPerMcCollision, mcCollision.globalIndex()); - if (trackSlice.size() < 1 || partSlice.size() < 1) - continue; - - float cent = getCentrality(col); - + histos.fill(HIST("hVtxZ"), mcCollision.posZ()); + if (!mcCollision.has_mcCollision() || !isEventSelected(mcCollision)) + return; + float cent = getCentrality(mcCollision); if (cent > KCentMax) - continue; - float multPV = col.multNTracksPV(); - float vz = col.posZ(); + 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"), 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 : partSlice) { + for (const auto& particle : mcParticles) { if (!isParticleSelected(particle) || !particle.isPhysicalPrimary()) continue; @@ -1534,13 +1523,7 @@ struct RadialFlowDecorr { histos.fill(HIST("h3_AllPrimary_PrAll"), multPV, pt, eta); } - histos.fill(HIST("hVtxZ_after_sel"), col.posZ()); - histos.fill(HIST("hCentrality"), cent); - - histos.fill(HIST("Hist2D_globalTracks_PVTracks"), multPV, trackSlice.size()); - histos.fill(HIST("Hist2D_cent_nch"), trackSlice.size(), cent); - - for (const auto& track : trackSlice) { + for (const auto& track : mcTracks) { if (!isTrackSelected(track)) continue; @@ -1548,41 +1531,33 @@ struct RadialFlowDecorr { auto sign = track.sign(); fillNSigmaBefCut(track, cent); - int centBin = static_cast(cent); - int id = identifyTrack(track, centBin); + 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; - fillNSigmaAftCut(track, cent, isSpecies); - if (isp == kInclusiveIdx) { histos.fill(HIST("h3_AllReco"), multPV, pt, eta); if (track.has_mcParticle()) { auto mcP = track.mcParticle(); if (mcP.isPhysicalPrimary()) { - int mcPdg = std::abs(mcP.pdgCode()); - if (mcPdg == KPiPlus || mcPdg == KKPlus || mcPdg == KProton) { - 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 { - // Misidentified! Reconstructed track, but true particle is not pi/K/P - histos.fill(HIST("h3_RecoMatchedToPrimary_MisID"), multPV, pt, eta); - } + 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); } @@ -1736,50 +1711,41 @@ struct RadialFlowDecorr { } } } - } - } } PROCESS_SWITCH(RadialFlowDecorr, processGetEffHists, "process MC to calculate EffWeights", cfgRunGetEff); - void processMCFlat(aod::McCollisions const& mcColl, MyRun3MCCollisions const& collisions, TCs const& tracks, FilteredTCs const& /*filteredTracks*/) + void processMCFlat(MyRun3MCCollisions::iterator const& mcCollision, FilteredTCs const& mcTracks) { - for (const auto& mcCollision : mcColl) { - auto colSlice = collisions.sliceBy(colPerMcCollision, mcCollision.globalIndex()); - if (colSlice.size() != 1) - continue; - - for (const auto& col : colSlice) { - if (!col.has_mcCollision() || !isEventSelected(col)) - continue; + histos.fill(HIST("hVtxZ"), mcCollision.posZ()); + if (!mcCollision.has_mcCollision() || !isEventSelected(mcCollision)) + return; - float cent = getCentrality(col); + float cent = getCentrality(mcCollision); if (cent > KCentMax) - continue; - - auto trackSlice = tracks.sliceBy(trackPerCollision, col.globalIndex()); - if (trackSlice.size() < 1) - continue; + return; - float multPV = col.multNTracksPV(); - float vz = col.posZ(); + float multPV = mcCollision.multNTracksPV(); + float vz = mcCollision.posZ(); - histos.fill(HIST("hVtxZ_after_sel"), col.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"), col.multNTracksPV(), trackSlice.size()); - histos.fill(HIST("Hist2D_cent_nch"), trackSlice.size(), cent); - for (const auto& track : trackSlice) { + 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 (!isTrackSelected(track)) continue; float pt = track.pt(), eta = track.eta(), phi = track.phi(); auto sign = track.sign(); - int centBin = static_cast(cent); - int id = identifyTrack(track, centBin); + 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, @@ -1839,36 +1805,32 @@ struct RadialFlowDecorr { } } } - } - } } PROCESS_SWITCH(RadialFlowDecorr, processMCFlat, "process MC to calculate FlatWeights", cfgRunGetMCFlat); - void processMCMean(aod::McCollisions const& mcColl, MyRun3MCCollisions const& collisions, TCs const& tracks, FilteredTCs const& /*filteredTracks*/, aod::FT0s const&, aod::McParticles const& mcParticles) + 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]{}; + 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; - for (const auto& mcCollision : mcColl) { - auto colSlice = collisions.sliceBy(colPerMcCollision, mcCollision.globalIndex()); - if (colSlice.size() != 1) - continue; - - for (const auto& col : colSlice) { - if (!col.has_mcCollision() || !isEventSelected(col)) - continue; - - auto trackSlice = tracks.sliceBy(trackPerCollision, col.globalIndex()); - auto partSlice = mcParticles.sliceBy(partPerMcCollision, mcCollision.globalIndex()); - if (trackSlice.size() < 1 || partSlice.size() < 1) - continue; + histos.fill(HIST("hVtxZ_after_sel"), mcCollision.posZ()); + histos.fill(HIST("hCentrality"), cent); - float cent = getCentrality(col); - if (cent > KCentMax) - continue; - float multPV = col.multNTracksPV(); - float vz = col.posZ(); + 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); memset(sumWiTruth, 0, sizeof(sumWiTruth)); memset(sumWiptiTruth, 0, sizeof(sumWiptiTruth)); @@ -1877,7 +1839,7 @@ struct RadialFlowDecorr { memset(sumWiRecoEffCorr, 0, sizeof(sumWiRecoEffCorr)); memset(sumWiptiRecoEffCorr, 0, sizeof(sumWiptiRecoEffCorr)); - for (const auto& particle : partSlice) { + for (const auto& particle : mcParticles) { if (!isParticleSelected(particle) || !particle.isPhysicalPrimary()) continue; float pt = particle.pt(), eta = particle.eta(); @@ -1921,25 +1883,17 @@ struct RadialFlowDecorr { } } - histos.fill(HIST("hVtxZ_after_sel"), col.posZ()); - histos.fill(HIST("hCentrality"), cent); - - histos.fill(HIST("Hist2D_globalTracks_PVTracks"), col.multNTracksPV(), trackSlice.size()); - histos.fill(HIST("Hist2D_cent_nch"), trackSlice.size(), cent); - - for (const auto& track : trackSlice) { + for (const auto& track : mcTracks) { if (!isTrackSelected(track)) continue; float pt = track.pt(), eta = track.eta(), phi = track.phi(); if (pt <= cfgPtMin || pt > cfgPtMax) continue; auto sign = track.sign(); - int centBin = static_cast(cent); - int id = identifyTrack(track, centBin); + 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, @@ -1982,7 +1936,6 @@ struct RadialFlowDecorr { 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); @@ -2146,8 +2099,8 @@ struct RadialFlowDecorr { } // end ietaA double amplFT0A = 0, amplFT0C = 0; - if (col.has_foundFT0()) { - const auto& ft0 = col.foundFT0(); + if (mcCollision.has_foundFT0()) { + const auto& ft0 = mcCollision.foundFT0(); for (std::size_t iCh = 0; iCh < ft0.channelA().size(); iCh++) { auto chanelid = ft0.channelA()[iCh]; float ampl = ft0.amplitudeA()[iCh]; @@ -2171,12 +2124,10 @@ struct RadialFlowDecorr { 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); - void processMCFluc(aod::McCollisions const& mcColl, MyRun3MCCollisions const& collisions, TCs const& tracks, FilteredTCs const& /*filteredTracks*/, aod::FT0s const&, aod::McParticles const& mcParticles) + 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) { @@ -2201,24 +2152,23 @@ struct RadialFlowDecorr { double p1kBarTru[KNsp][KNEta]{}, p1kBarReco[KNsp][KNEta]{}, p1kBarRecoEffCor[KNsp][KNEta]{}; double p1kBarTruMult[KNsp][KNEta]{}, p1kBarRecoMult[KNsp][KNEta]{}, p1kBarRecoEffCorMult[KNsp][KNEta]{}; - for (const auto& mcCollision : mcColl) { - auto partSlice = mcParticles.sliceBy(partPerMcCollision, mcCollision.globalIndex()); - auto colSlice = collisions.sliceBy(colPerMcCollision, mcCollision.globalIndex()); - if (colSlice.size() != 1) - continue; - - for (const auto& col : colSlice) { - if (!col.has_mcCollision() || !isEventSelected(col)) - continue; + 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; - auto trackSlice = tracks.sliceBy(trackPerCollision, col.globalIndex()); - if (trackSlice.size() < 1) - continue; + histos.fill(HIST("hVtxZ_after_sel"), mcCollision.posZ()); + histos.fill(HIST("hCentrality"), cent); - float cent = getCentrality(col); - if (cent > KCentMax) - continue; - float multPV = col.multNTracksPV(); + 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)); @@ -2248,7 +2198,9 @@ struct RadialFlowDecorr { double p1kBarFt0A = 0.0, p1kBarFt0C = 0.0; - for (const auto& particle : partSlice) { + + + for (const auto& particle : mcParticles) { if (!isParticleSelected(particle) || !particle.isPhysicalPrimary()) continue; @@ -2287,15 +2239,9 @@ struct RadialFlowDecorr { } } } // end truth loop - float vz = col.posZ(); - - histos.fill(HIST("hVtxZ_after_sel"), col.posZ()); - histos.fill(HIST("hCentrality"), cent); - histos.fill(HIST("Hist2D_globalTracks_PVTracks"), col.multNTracksPV(), trackSlice.size()); - histos.fill(HIST("Hist2D_cent_nch"), trackSlice.size(), cent); - for (const auto& track : trackSlice) { + for (const auto& track : mcTracks) { if (!isTrackSelected(track)) continue; @@ -2305,12 +2251,10 @@ struct RadialFlowDecorr { float eta = track.eta(); float phi = track.phi(); auto sign = track.sign(); - int centBin = static_cast(cent); - int id = identifyTrack(track, centBin); + 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, @@ -2320,8 +2264,8 @@ struct RadialFlowDecorr { for (int isp = 0; isp < KNsp; ++isp) { if (!isSpecies[isp]) continue; - float eff = getEfficiency(col.multNTracksPV(), pt, eta, static_cast(isp), 0, cfgEff); - float fake = getEfficiency(col.multNTracksPV(), pt, eta, static_cast(isp), 1, cfgEff); + 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; @@ -2386,7 +2330,7 @@ struct RadialFlowDecorr { } // trkslice for (int ieta = 0; ieta < KNEta; ++ieta) { - const int ibx = state.pmeanTruNchEtabinSpbinStep2->GetXaxis()->FindBin(col.multNTracksPV()); + const int ibx = state.pmeanTruNchEtabinSpbinStep2->GetXaxis()->FindBin(mcCollision.multNTracksPV()); const int iby = ieta + 1; for (int isp = 0; isp < KNsp; ++isp) { @@ -2428,8 +2372,8 @@ struct RadialFlowDecorr { } double amplFT0A = 0, amplFT0C = 0; - if (col.has_foundFT0()) { - const auto& ft0 = col.foundFT0(); + 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; @@ -2470,33 +2414,33 @@ struct RadialFlowDecorr { 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"), col.multNTracksPV(), 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"), col.multNTracksPV(), 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"), col.multNTracksPV(), 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"), col.multNTracksPV(), 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"), col.multNTracksPV(), 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"), col.multNTracksPV(), ieta, isp, c2RecoEffCor[isp][ieta]); + histos.fill(HIST("MCRecoEffCorr/Prof_C2_Mult_etabin_spbin"), multPV, ieta, isp, c2RecoEffCor[isp][ieta]); } } } - p1kBarFt0A = amplFT0A - state.pmeanFT0AmultpvStep2->GetBinContent(state.pmeanFT0AmultpvStep2->GetXaxis()->FindBin(col.multNTracksPV())); - p1kBarFt0C = amplFT0C - state.pmeanFT0CmultpvStep2->GetBinContent(state.pmeanFT0CmultpvStep2->GetXaxis()->FindBin(col.multNTracksPV())); + p1kBarFt0A = amplFT0A - state.pmeanFT0AmultpvStep2->GetBinContent(state.pmeanFT0AmultpvStep2->GetXaxis()->FindBin(multPV)); + p1kBarFt0C = amplFT0C - state.pmeanFT0CmultpvStep2->GetBinContent(state.pmeanFT0CmultpvStep2->GetXaxis()->FindBin(multPV)); for (int ietaA = 1; ietaA <= (KNEta - 1) / 2; ++ietaA) { int ietaC = KNEta - ietaA; @@ -2519,53 +2463,53 @@ struct RadialFlowDecorr { 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"), col.multNTracksPV(), ietaA, isp, c2SubTru); + histos.fill(HIST("MCGen/Prof_C2Sub_Mult_etabin_spbin"), multPV, ietaA, isp, 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"), col.multNTracksPV(), ietaA, isp, c2SubReco); + histos.fill(HIST("MCReco/Prof_C2Sub_Mult_etabin_spbin"), multPV, ietaA, isp, 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"), col.multNTracksPV(), ietaA, isp, c2SubRecoEffCor); + histos.fill(HIST("MCRecoEffCorr/Prof_C2Sub_Mult_etabin_spbin"), multPV, ietaA, isp, 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"), col.multNTracksPV(), ietaA, isp, covTru); + histos.fill(HIST("MCGen/Prof_Cov_Mult_etabin_spbin"), multPV, ietaA, isp, 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"), col.multNTracksPV(), ietaA, isp, covReco); + histos.fill(HIST("MCReco/Prof_Cov_Mult_etabin_spbin"), multPV, ietaA, isp, 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"), col.multNTracksPV(), ietaA, isp, covRecoEffCor); + histos.fill(HIST("MCRecoEffCorr/Prof_Cov_Mult_etabin_spbin"), multPV, ietaA, isp, 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"), col.multNTracksPV(), ietaA, isp, covFT0ATru); + histos.fill(HIST("MCGen/Prof_CovFT0A_Mult_etabin_spbin"), multPV, ietaA, isp, 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"), col.multNTracksPV(), ietaA, isp, covFT0AReco); + histos.fill(HIST("MCReco/Prof_CovFT0A_Mult_etabin_spbin"), multPV, ietaA, isp, 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"), col.multNTracksPV(), ietaA, isp, covFT0ARecoEffCor); + histos.fill(HIST("MCRecoEffCorr/Prof_CovFT0A_Mult_etabin_spbin"), multPV, ietaA, isp, 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"), col.multNTracksPV(), ietaA, isp, covFT0CTru); + histos.fill(HIST("MCGen/Prof_CovFT0C_Mult_etabin_spbin"), multPV, ietaA, isp, 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"), col.multNTracksPV(), ietaA, isp, covFT0CReco); + histos.fill(HIST("MCReco/Prof_CovFT0C_Mult_etabin_spbin"), multPV, ietaA, isp, 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"), col.multNTracksPV(), ietaA, isp, covFT0CRecoEffCor); + histos.fill(HIST("MCRecoEffCorr/Prof_CovFT0C_Mult_etabin_spbin"), multPV, ietaA, isp, covFT0CRecoEffCor); } } } @@ -2948,8 +2892,6 @@ struct RadialFlowDecorr { } } } - } // colSlice - } // mcColl LOGF(info, "FINISHED RUNNING processMCFluc"); } PROCESS_SWITCH(RadialFlowDecorr, processMCFluc, "process MC to calculate pt fluc", cfgRunMCFluc); @@ -2967,6 +2909,8 @@ struct RadialFlowDecorr { 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) { @@ -3006,11 +2950,16 @@ struct RadialFlowDecorr { if (cent > KCentMax) return; + 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()); int ntrk = 0; float vz = coll.posZ(); @@ -3028,22 +2977,20 @@ struct RadialFlowDecorr { if (eta > etaLw[0] && eta < etaUp[0]) ntrk++; fillNSigmaBefCut(track, cent); - int centBin = static_cast(cent); - int id = identifyTrack(track, centBin); + 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; - fillNSigmaAftCut(track, cent, isSpecies); float eff = getEfficiency(coll.multNTracksPV(), pt, eta, static_cast(isp), 0, cfgEff); if (eff <= KFloatEpsilon) continue; @@ -3125,11 +3072,16 @@ struct RadialFlowDecorr { if (cent > KCentMax) return; + 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()); float vz = coll.posZ(); @@ -3153,12 +3105,10 @@ struct RadialFlowDecorr { histos.fill(HIST("hPt"), pt); histos.fill(HIST("hEta"), eta); histos.fill(HIST("hPhi"), phi); - int centBin = static_cast(cent); - int id = identifyTrack(track, centBin); + 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, @@ -3169,8 +3119,6 @@ struct RadialFlowDecorr { if (!isSpecies[isp]) continue; float eff = getEfficiency(coll.multNTracksPV(), pt, eta, static_cast(isp), 0, cfgEff); - - // Safety check BEFORE dividing if (eff <= KFloatEpsilon) continue; @@ -3233,16 +3181,20 @@ struct RadialFlowDecorr { } 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]); - if (sumWi[isp][0] > 1.0f) - histos.fill(HIST("Prof_Cent_Nsp_MeanpT"), cent, isp, sumWipti[isp][0] / 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; @@ -3279,7 +3231,6 @@ struct RadialFlowDecorr { } } } - double amplFT0A = 0, amplFT0C = 0; if (coll.has_foundFT0()) { const auto& ft0 = coll.foundFT0(); @@ -3317,6 +3268,17 @@ struct RadialFlowDecorr { if (cent > KCentMax) return; + 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"); return; @@ -3354,12 +3316,10 @@ struct RadialFlowDecorr { if (pt <= cfgPtMin || pt > cfgPtMax) continue; - int centBin = static_cast(cent); - int id = identifyTrack(track, centBin); + 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, @@ -3370,8 +3330,6 @@ struct RadialFlowDecorr { if (!isSpecies[isp]) continue; float eff = getEfficiency(coll.multNTracksPV(), pt, eta, static_cast(isp), 0, cfgEff); - - // Safety check BEFORE dividing if (eff <= KFloatEpsilon) continue; @@ -3491,7 +3449,6 @@ struct RadialFlowDecorr { float covFT0A = p1kBarFt0A * p1kBar[isp][ietaC]; float covFT0C = p1kBarFt0C * p1kBar[isp][ietaA]; - // Updated enum checks here if (isp == kInclusiveIdx) { if (std::isfinite(c2Sub)) { histos.fill(HIST("Prof_C2Sub2D_Cent_etaA_etaC"), cent, etaValA, etaValB, c2Sub); From 5f57b16cbffec7f6dc3d9961f20ef2ecb41ec1d0 Mon Sep 17 00:00:00 2001 From: Somadutta Bhatta Date: Fri, 13 Mar 2026 17:19:22 +0100 Subject: [PATCH 2/3] Resolve merge conflict --- PWGCF/EbyEFluctuations/Tasks/radialFlowDecorr.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PWGCF/EbyEFluctuations/Tasks/radialFlowDecorr.cxx b/PWGCF/EbyEFluctuations/Tasks/radialFlowDecorr.cxx index bc7c5078594..761fe3e2d12 100644 --- a/PWGCF/EbyEFluctuations/Tasks/radialFlowDecorr.cxx +++ b/PWGCF/EbyEFluctuations/Tasks/radialFlowDecorr.cxx @@ -3121,7 +3121,7 @@ struct RadialFlowDecorr { 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; From d4c2cfd42543f5746f1f8ee6febf5e382707bc52 Mon Sep 17 00:00:00 2001 From: Somadutta Bhatta Date: Fri, 13 Mar 2026 17:22:17 +0100 Subject: [PATCH 3/3] Resolve merge conflicts --- PWGCF/EbyEFluctuations/Tasks/radialFlowDecorr.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PWGCF/EbyEFluctuations/Tasks/radialFlowDecorr.cxx b/PWGCF/EbyEFluctuations/Tasks/radialFlowDecorr.cxx index 761fe3e2d12..bc7c5078594 100644 --- a/PWGCF/EbyEFluctuations/Tasks/radialFlowDecorr.cxx +++ b/PWGCF/EbyEFluctuations/Tasks/radialFlowDecorr.cxx @@ -3121,7 +3121,7 @@ struct RadialFlowDecorr { 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;