diff --git a/PWGCF/EbyEFluctuations/Tasks/radialFlowDecorr.cxx b/PWGCF/EbyEFluctuations/Tasks/radialFlowDecorr.cxx index 33b16fb7c06..4f710a0e788 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}; @@ -202,13 +202,20 @@ struct RadialFlowDecorr { 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, "Linear Pileup Cut Const"}; + Configurable cfgLinPupParam1{"cfgLinPupParam1", 3.0f, "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 +273,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,7 +297,6 @@ 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()); @@ -397,39 +406,62 @@ struct RadialFlowDecorr { // 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 +469,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 +476,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 +492,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,6 +542,29 @@ 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); return true; } @@ -771,20 +824,22 @@ struct RadialFlowDecorr { 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("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"); histos.add("hTrkCount", "Number of Tracks;; Count", kTH1F, {{11, 0, 11}}); histos.get(HIST("hTrkCount"))->GetXaxis()->SetBinLabel(1, "all Tracks"); @@ -1193,50 +1248,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 +1297,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); } } @@ -1446,15 +1462,16 @@ struct RadialFlowDecorr { if (cent > KCentMax) continue; float multPV = col.multNTracksPV(); + + 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& particle : partSlice) { if (!isParticleSelected(particle) || !particle.isPhysicalPrimary()) continue; - 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) { if (!isTrackSelected(track)) continue; @@ -1490,6 +1507,15 @@ struct RadialFlowDecorr { float multPV = col.multNTracksPV(); float vz = col.posZ(); + if (!isPassAddPileup(multPV, trackSlice.size(), cent)) + continue; + + 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& particle : partSlice) { if (!isParticleSelected(particle) || !particle.isPhysicalPrimary()) continue; @@ -1499,16 +1525,16 @@ struct RadialFlowDecorr { 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 + (absPdg == KPiPlus || absPdg == KKPlus || absPdg == KProton), // 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); @@ -1534,12 +1560,6 @@ 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) { if (!isTrackSelected(track)) continue; @@ -1548,14 +1568,14 @@ 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 isAny = (isPi || isKa || isPr); // Only true if it passed PID! bool isSpecies[KNsp] = { - true, + isAny, isPi && sign < 0, isPi && sign > 0, isPi, isKa && sign < 0, isKa && sign > 0, isKa, isPr && sign < 0, isPr && sign > 0, isPr}; @@ -1749,6 +1769,7 @@ struct RadialFlowDecorr { continue; for (const auto& col : colSlice) { + histos.fill(HIST("hVtxZ"), col.posZ()); if (!col.has_mcCollision() || !isEventSelected(col)) continue; @@ -1763,6 +1784,8 @@ struct RadialFlowDecorr { float multPV = col.multNTracksPV(); float vz = col.posZ(); + if (!isPassAddPileup(multPV, trackSlice.size(), cent)) + continue; histos.fill(HIST("hVtxZ_after_sel"), col.posZ()); histos.fill(HIST("hCentrality"), cent); @@ -1774,14 +1797,14 @@ struct RadialFlowDecorr { 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 isAny = (isPi || isKa || isPr); // Only true if it passed PID! bool isSpecies[KNsp] = { - true, + isAny, isPi && sign < 0, isPi && sign > 0, isPi, isKa && sign < 0, isKa && sign > 0, isKa, isPr && sign < 0, isPr && sign > 0, isPr}; @@ -1856,6 +1879,7 @@ struct RadialFlowDecorr { continue; for (const auto& col : colSlice) { + histos.fill(HIST("hVtxZ"), col.posZ()); if (!col.has_mcCollision() || !isEventSelected(col)) continue; @@ -1870,6 +1894,15 @@ struct RadialFlowDecorr { float multPV = col.multNTracksPV(); float vz = col.posZ(); + if (!isPassAddPileup(multPV, trackSlice.size(), cent)) + continue; + + 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); + memset(sumWiTruth, 0, sizeof(sumWiTruth)); memset(sumWiptiTruth, 0, sizeof(sumWiptiTruth)); memset(sumWiReco, 0, sizeof(sumWiReco)); @@ -1887,16 +1920,16 @@ struct RadialFlowDecorr { 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 + (absPdg == KPiPlus || absPdg == KKPlus || absPdg == KProton), // 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) { @@ -1921,12 +1954,6 @@ 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) { if (!isTrackSelected(track)) continue; @@ -1934,14 +1961,13 @@ struct RadialFlowDecorr { 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 isAny = (isPi || isKa || isPr); // Only true if it passed PID! bool isSpecies[KNsp] = { - true, + isAny, isPi && sign < 0, isPi && sign > 0, isPi, isKa && sign < 0, isKa && sign > 0, isKa, isPr && sign < 0, isPr && sign > 0, isPr}; @@ -2220,6 +2246,15 @@ struct RadialFlowDecorr { continue; float multPV = col.multNTracksPV(); + if (!isPassAddPileup(multPV, trackSlice.size(), cent)) + continue; + + 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); + memset(sumPmwkTru, 0, sizeof(sumPmwkTru)); memset(sumWkTru, 0, sizeof(sumWkTru)); memset(sumPmwkReco, 0, sizeof(sumPmwkReco)); @@ -2260,16 +2295,16 @@ struct RadialFlowDecorr { 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 + (absPdg == KPiPlus || absPdg == KKPlus || absPdg == KProton), // 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) { @@ -2289,12 +2324,6 @@ 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) { if (!isTrackSelected(track)) continue; @@ -2305,14 +2334,14 @@ 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 isAny = (isPi || isKa || isPr); // Only true if it passed PID! bool isSpecies[KNsp] = { - true, + isAny, isPi && sign < 0, isPi && sign > 0, isPi, isKa && sign < 0, isKa && sign > 0, isKa, isPr && sign < 0, isPr && sign > 0, isPr}; @@ -3006,6 +3035,9 @@ 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); @@ -3028,14 +3060,14 @@ 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 isAny = (isPi || isKa || isPr); // Only true if it passed PID! bool isSpecies[KNsp] = { - true, + isAny, isPi && sign < 0, isPi && sign > 0, isPi, isKa && sign < 0, isKa && sign > 0, isKa, isPr && sign < 0, isPr && sign > 0, isPr}; @@ -3125,6 +3157,9 @@ 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); @@ -3153,14 +3188,14 @@ 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 isAny = (isPi || isKa || isPr); // Only true if it passed PID! bool isSpecies[KNsp] = { - true, + isAny, isPi && sign < 0, isPi && sign > 0, isPi, isKa && sign < 0, isKa && sign > 0, isKa, isPr && sign < 0, isPr && sign > 0, isPr}; @@ -3317,6 +3352,15 @@ 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); + if (!state.pmeanNchEtabinSpbinStep2 || !state.pmeanMultNchEtabinSpbinStep2) { LOGF(warning, "Data fluc: Mean pT or Mult map missing"); return; @@ -3354,14 +3398,14 @@ 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 isAny = (isPi || isKa || isPr); // Only true if it passed PID! bool isSpecies[KNsp] = { - true, + isAny, isPi && sign < 0, isPi && sign > 0, isPi, isKa && sign < 0, isKa && sign > 0, isKa, isPr && sign < 0, isPr && sign > 0, isPr};