From 893b14cc8d07f8c9f356d527ba3600b5b6ab7dc7 Mon Sep 17 00:00:00 2001 From: Gyula Bencedi Date: Fri, 11 Sep 2026 20:13:23 +0200 Subject: [PATCH] Split MC process + optimize histo filling --- .../GlobalEventProperties/flattenictyPikp.cxx | 691 +++++++----------- 1 file changed, 274 insertions(+), 417 deletions(-) diff --git a/PWGLF/Tasks/GlobalEventProperties/flattenictyPikp.cxx b/PWGLF/Tasks/GlobalEventProperties/flattenictyPikp.cxx index 25549595664..136f76845be 100644 --- a/PWGLF/Tasks/GlobalEventProperties/flattenictyPikp.cxx +++ b/PWGLF/Tasks/GlobalEventProperties/flattenictyPikp.cxx @@ -58,7 +58,6 @@ #include #include #include -#include #include @@ -84,9 +83,7 @@ using namespace o2::constants::math; using namespace o2::aod::rctsel; auto static constexpr CminCharge = 3.f; -static constexpr int CnullInt = 0; static constexpr float Cnull = 0.0f; -static constexpr int ConeInt = 1; static constexpr float Cone = 1.0f; // FV0 specific constants @@ -97,19 +94,11 @@ const std::array cFV0IndexPhi{0, 8, 16, 24, 32}; static constexpr float CmaxEtaFV0 = 5.1; static constexpr float CminEtaFV0 = 2.2; static constexpr float CdEtaFV0 = (CmaxEtaFV0 - CminEtaFV0) / CmaxRingsFV0; -auto static constexpr CminAccFT0A = 3.5f; -auto static constexpr CmaxAccFT0A = 4.9f; -auto static constexpr CminAccFT0C = -3.3f; -auto static constexpr CmaxAccFT0C = -2.1f; // PID names static constexpr int CprocessIdWeak = 4; -static constexpr int Ncharges = 2; static constexpr o2::track::PID::ID Npart = 5; -static constexpr o2::track::PID::ID NpartChrg = Npart * Ncharges; const std::array pDGs{11, 13, 211, 321, 2212}; -const std::array pIdSgn{11, 13, 211, 321, 2212, -11, -13, -211, -321, -2212}; -const std::array pID{"el", "mu", "pi", "ka", "pr"}; static constexpr std::array CspeciesAll{"El", "Mu", "Pi", "Ka", "Pr"}; // histogram naming @@ -218,7 +207,9 @@ enum EvtSel { enum class MCEventStatus { kMcEvtAll = 0, - kMcEvtSplit, + kEvtSelected, + kEvtHasMcColl, + kEvtSplitVtxRemoved, nEvtMcEventStatus }; @@ -253,7 +244,7 @@ struct FlattenictyPikp { o2::parameters::GRPMagField* grpmag = nullptr; struct : ConfigurableGroup { - Configurable multEst{"multEst", 1, "0: without multiplicity; 1: MultFT0C; 2: MultFT0M; 3: MultTPC"}; + Configurable multEst{"multEst", 2, "0: without multiplicity; 1: MultFT0C; 2: MultFT0M; 3: MultTPC"}; Configurable applyCalibGainFromCCDB{"applyCalibGainFromCCDB", false, "equalize detector amplitudes"}; Configurable applyCalibVtxFromCCDB{"applyCalibVtxFromCCDB", false, "equalize Amp vs vtx"}; Configurable applyCalibDeDx{"applyCalibDeDx", false, "calibration of dedx signal"}; @@ -268,8 +259,6 @@ struct FlattenictyPikp { Configurable fillV0Hist{"fillV0Hist", false, "fill V0 histograms"}; Configurable fillChrgType{"fillChrgType", false, "fill histograms per charge types"}; Configurable fillChrgTypeV0s{"fillChrgTypeV0s", false, "fill V0s histograms per charge types"}; - Configurable fillMCRecCheck{"fillMCRecCheck", false, "fill MC rec histograms"}; - Configurable fillMCRecDCA{"fillMCRecDCA", false, "fill MC rec DCA histograms"}; Configurable calibDeDxFunction{"calibDeDxFunction", "pol8", "Functional form for dEdx calibration"}; Configurable> paramsFuncMIPposEtaP{"paramsFuncMIPposEtaP", std::vector{-1.f}, "function parameters"}; Configurable> paramsFuncMIPnegEtaP{"paramsFuncMIPnegEtaP", std::vector{-1.f}, "function parameters"}; @@ -312,7 +301,7 @@ struct FlattenictyPikp { ConfigurableAxis axisPtV0s{"axisPtV0s", {VARIABLE_WIDTH, 0, 0.15, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1, 1.2, 1.4, 1.6, 1.8, 2, 2.2, 2.4, 2.6, 2.8, 3, 3.2, 3.4, 3.6, 3.8, 4, 4.4, 4.8, 5.2, 5.6, 6, 6.4, 6.8, 7.2, 7.6, 8, 8.4, 8.8, 9.2, 9.6, 10, 10.5, 11, 11.5, 12, 12.5, 13, 13.5, 14, 14.5, 15, 20}, "pT V0s binning"}; ConfigurableAxis axisPtFine{"axisPtFine", {250, 0.1f, 20.1f}, "pT fine binning"}; ConfigurableAxis axisFlatPerc{"axisFlatPerc", {102, -0.01, 1.01}, "Flattenicity percentiles binning"}; - ConfigurableAxis axisMultPerc{"axisMultPerc", {VARIABLE_WIDTH, 0., 5., 10., 20., 30., 40., 50., 60., 70., 80., 90., 100.}, "T0 percentiles binning"}; + ConfigurableAxis axisMultPerc{"axisMultPerc", {VARIABLE_WIDTH, 0., 1., 5., 10., 15., 20., 30., 40., 50., 70., 100.}, "Multiplicity percentiles binning"}; ConfigurableAxis axisVertexZ{"axisVertexZ", {80, -20., 20.}, "Vertex z binning"}; ConfigurableAxis axisMult{"axisMult", {301, -0.5, 300.5}, "Multiplicity binning"}; ConfigurableAxis axisDCAxy{"axisDCAxy", {210, -1.05f, 1.05f}, "DCAxy binning"}; @@ -325,7 +314,7 @@ struct FlattenictyPikp { ConfigurableAxis axisNsigmaTPC{"axisNsigmaTPC", {200, -10, 10}, "nsigmaTPC binning"}; ConfigurableAxis axisNsigmaTOF{"axisNsigmaTOF", {200, -10, 10}, "nsigmaTOF binning"}; ConfigurableAxis axisAmplFV0{"axisAmplFV0", {4096, 0, 4096}, "FV0 amplitude (ADC) binning"}; - ConfigurableAxis axisAmplFV0Sum{"axisAmplFV0Sum", {4096, 0, 4096 * 49}, "FV0 amplitude sum (ADC) binning"}; + ConfigurableAxis axisAmplFV0Sum{"axisAmplFV0Sum", {4096, 0, 4096}, "FV0 amplitude sum (ADC) binning"}; ConfigurableAxis axisChannelFV0{"axisChannelFV0", {49, 0., 49.}, "FV0 channel ID binning"}; } binOpt; @@ -334,9 +323,9 @@ struct FlattenictyPikp { Configurable rapMax{"rapMax", 0.5f, "Maximum range of rapidity for tracks"}; Configurable trkPtMin{"trkPtMin", 0.1f, "Minimum pT of tracks"}; Configurable applyNcl{"applyNcl", false, "Apply cut on TPC clusters"}; - Configurable nclTPCMin{"nclTPCMin", 135.0f, "Minimum of number of TPC found clusters"}; + Configurable nclTPCMin{"nclTPCMin", 130.0f, "Minimum of number of TPC found clusters"}; Configurable applyNclPID{"applyNclPID", true, "Apply cut on TPC PID clusters"}; - Configurable nclPidTPCMin{"nclPidTPCMin", 135.0f, "Minimum of number of TPC PID clusters"}; + Configurable nclPidTPCMin{"nclPidTPCMin", 130.0f, "Minimum of number of TPC PID clusters"}; Configurable phiCutPtMin{"phiCutPtMin", 2.0f, "Minimum pT for phi cut"}; Configurable tofBetaPion{"tofBetaPion", 1.0f, "Minimum beta for TOF pions"}; Configurable tofBetaPiMax{"tofBetaPiMax", 5E-5, "Maximum beta for TOF pion selection"}; @@ -362,7 +351,7 @@ struct FlattenictyPikp { struct : ConfigurableGroup { // common selection Configurable v0TypeSel{"v0TypeSel", 1, "select on a certain V0 type (leave negative if no selection desired)"}; - Configurable v0Ymax{"v0Ymax", 0.5f, "Maximum rapidity of V0s"}; + Configurable v0Ymax{"v0Ymax", 0.8f, "Maximum rapidity of V0s"}; Configurable rejectV0sAtTPCSector{"rejectV0sAtTPCSector", true, "Reject V0s close to the TPC sector boundaries"}; Configurable v0requireITS{"v0requireITS", true, "Additional cut on the ITS requirement"}; Configurable nsigmaElTPC{"nsigmaElTPC", 5.0, "max nsigma of TPC for electorn"}; @@ -378,9 +367,9 @@ struct FlattenictyPikp { Configurable v0etamax{"v0etamax", +0.8f, "max eta of V0s"}; Configurable v0minNCrossedRowsTPC{"v0minNCrossedRowsTPC", 70, "Additional cut on the minimum number of crossed rows in the TPC"}; Configurable applyV0sNclFound{"applyV0sNclFound", false, "Apply cut on TPC Found clusters"}; - Configurable v0NclTPCMin{"v0NclTPCMin", 135.0f, "Minimum of number of TPC found clusters"}; + Configurable v0NclTPCMin{"v0NclTPCMin", 130.0f, "Minimum of number of TPC found clusters"}; Configurable applyV0sNclPID{"applyV0sNclPID", true, "Apply cut on TPC PID clusters"}; - Configurable v0NclPidTPCMin{"v0NclPidTPCMin", 135.0f, "Minimum of number of TPC PID clusters"}; + Configurable v0NclPidTPCMin{"v0NclPidTPCMin", 130.0f, "Minimum of number of TPC PID clusters"}; Configurable v0maxChi2PerClusterTPC{"v0maxChi2PerClusterTPC", 4.f, "Additional cut on the maximum value of the chi2 per cluster in the TPC"}; Configurable v0maxChi2PerClusterITS{"v0maxChi2PerClusterITS", 36.f, "Additional cut on the maximum value of the chi2 per cluster in the ITS"}; Configurable v0minITSnClusters{"v0minITSnClusters", 4, "minimum number of found ITS clusters"}; @@ -464,7 +453,7 @@ struct FlattenictyPikp { Configurable maxChi2PerClusterITS{"maxChi2PerClusterITS", 36.f, "Additional cut on the maximum value of the chi2 per cluster in the ITS"}; Configurable minITSnClusters{"minITSnClusters", 5, "minimum number of found ITS clusters"}; Configurable maxDcaXYFactor{"maxDcaXYFactor", 1.f, "Multiplicative factor on the maximum value of the DCA xy"}; - Configurable maxDcaZ{"maxDcaZ", 2.f, "Additional cut on the maximum value of the DCA z"}; + Configurable maxDcaZ{"maxDcaZ", 0.1f, "Additional cut on the maximum value of the DCA z"}; TF1* fPhiCutLow = nullptr; TF1* fPhiCutHigh = nullptr; @@ -474,27 +463,6 @@ struct FlattenictyPikp { std::vector> fEDeDxVsEta; std::vector> vecParamsPLA; - std::array, NpartChrg> hPtGenRecEvt{}; - std::array, NpartChrg> hPtGenPrimRecEvt{}; - std::array, NpartChrg> hPtGenRecEvtGtZero{}; - std::array, NpartChrg> hPtGenPrimRecEvtGtZero{}; - std::array, NpartChrg> hPtEffRecGoodCollPrimSelTrk{}; - std::array, NpartChrg> hPtEffRecGoodCollPrim{}; - std::array, NpartChrg> hPtEffRecGoodCollWeak{}; - std::array, NpartChrg> hPtEffRecGoodCollMat{}; - std::array, NpartChrg> hPtEffGenPrim{}; - std::array, NpartChrg> hPtEffGenWeak{}; - std::array, NpartChrg> hPtEffGenMat{}; - std::array, NpartChrg> hPtEffGenPrimEvtSelGen{}; - std::array, NpartChrg> hPtEffGenWeakEvtSelGen{}; - std::array, NpartChrg> hPtEffGenMatEvtSelGen{}; - std::array, NpartChrg> hDCAxyRecBadCollPrim{}; - std::array, NpartChrg> hDCAxyRecBadCollWeak{}; - std::array, NpartChrg> hDCAxyRecBadCollMat{}; - std::array, NpartChrg> hPtVsDCAxyRecGoodCollPrim{}; - std::array, NpartChrg> hPtVsDCAxyRecGoodCollWeak{}; - std::array, NpartChrg> hPtVsDCAxyRecGoodCollMat{}; - using MyCollisions = soa::Join; using Colls = soa::Join; using CollsGen = soa::Join; @@ -795,18 +763,47 @@ struct FlattenictyPikp { registryData.print(); } - if (doprocessMC) { - registryMC.add("Events/hMCEventStatus", "Number of events; Cut; #Evt Passed Cut", {HistType::kTH1F, {{static_cast(MCEventStatus::nEvtMcEventStatus), -0.5, +static_cast(MCEventStatus::nEvtMcEventStatus) - 0.5}}}); + if (doprocessMcEfficiency || doprocessMcRecTrack) { + if (defOpt.fillNclVsPhiCutQaHist) { + registryMC.add("Tracks/postSel/hPtPhi", "; #it{p}_{T} (GeV/#it{c}); fmod(#varphi,#pi/9)", {kTH2F, {ptAxis, phiAxisMod}}); + registryMC.add("Tracks/postSel/hPtPhiNclTPC", "; #eta; #it{p}_{T} (GeV/#it{c}); fmod(#varphi,#pi/9); N_{cluster}", {kTHnSparseF, {etaAxis, ptAxis, phiAxisMod, clTpcAxis}}); + registryMC.add("Tracks/postSel/hPtPhiNclPIDTPC", "; #eta; #it{p}_{T} (GeV/#it{c}); fmod(#varphi,#pi/9); N_{PID cluster}", {kTHnSparseF, {etaAxis, ptAxis, phiAxisMod, clTpcAxis}}); + registryMC.add("Tracks/postSel/hPtNclTPC", "; #it{p}_{T} (GeV/#it{c}); N_{cluster}", {kTH2F, {ptAxis, clTpcAxis}}); + registryMC.add("Tracks/postSel/pPtNclTPC", "; #it{p}_{T} (GeV/#it{c}); N_{cluster}", {kTProfile, {ptAxis}}); + registryMC.add("Tracks/postSel/hPtNclPIDTPC", "; #it{p}_{T} (GeV/#it{c}); N_{PID cluster}", {kTH2F, {ptAxis, clTpcAxis}}); + registryMC.add("Tracks/postSel/pPtNclPIDTPC", "; #it{p}_{T} (GeV/#it{c}); N_{PID cluster}", {kTProfile, {ptAxis}}); + } + registryMC.addClone("Tracks/postSel/", "Tracks/preSel/"); + } + + if (doprocessMcEfficiency) { + registryMC.add({"Events/ResponseGen", ";N_{ch,FV0};1-#rho_{FV0};", {kTHnSparseF, {multAxis, flatAxis}}}); + registryMC.add("Events/h1flatencityFV0MCGen", "", {kTH1F, {flatAxis}}); + registryMC.add("Events/hVtxZRec", "MC Rec vertex z position", kTH1F, {vtxzAxis}); + registryMC.add("Events/hVtxZGenGtOneRec", "MC Rec vertex z position w/ Nrec > 0", kTH1F, {vtxzAxis}); + registryMC.add("Events/hVtxZGen", "Generated vertex z position", kTH1F, {vtxzAxis}); + registryMC.add("Events/hMCEventStatus", "Number of events; Cut; #Evt Passed Cut", {HistType::kTHnSparseF, {{static_cast(MCEventStatus::nEvtMcEventStatus), -0.5, +static_cast(MCEventStatus::nEvtMcEventStatus) - 0.5}, multAxis, flatAxis}}); std::array(MCEventStatus::nEvtMcEventStatus)> labelEvtSel{ "All", - "Remove split vertices"}; - registryMC.get(HIST("Events/hMCEventStatus"))->SetMinimum(0.1); + "Selected", + "Has Mc Coll", + "Split Vtx Removed"}; for (int iBin = 0; iBin < static_cast(MCEventStatus::nEvtMcEventStatus); iBin++) { - registryMC.get(HIST("Events/hMCEventStatus"))->GetXaxis()->SetBinLabel(iBin + 1, labelEvtSel[iBin].data()); + registryMC.get(HIST("Events/hMCEventStatus"))->GetAxis(0)->SetBinLabel(iBin + 1, labelEvtSel[iBin].data()); } + registryMC.add({"Tracks/hPtFakes", "Fake tracks; Gen Nch (|#eta|<0.8); flat; #it{p}_{T} (GeV/#it{c})", {kTHnSparseF, {nChAxis, flatAxis, ptAxis}}}); + registryMC.add("Tracks/hPtRes", "#it{p}_{T} resolution;;(#it{p}_{T}_{rec} - #it{p}_{T}_{gen})/#it{p}_{T}_{gen};", kTH2F, {ptAxis, {100, -1.0, 1.0}}); - registryMC.add({"Events/ResponseGen", ";N_{ch,FV0};1-#rho_{FV0};", {kTHnSparseF, {multAxis, flatAxis}}}); - registryMC.add("Events/h1flatencityFV0MCGen", "", {kTH1F, {flatAxis}}); + for (int i = 0; i < Npart; i++) { + registryMC.add({fmt::format(CpTeffGenPrimRecEvtF.data(), CspeciesAll[i].data()).c_str(), "Gen evt w/ Nrec > 0; mult; flat; #it{p}_{T} (GeV/#it{c})", {kTHnSparseF, {multAxis, flatAxis, ptAxis}}}); + registryMC.add({fmt::format(CpTeffPrimRecEvtF.data(), CspeciesAll[i].data()).c_str(), "Gen evt w/ Nrec > 0 + Evt sel; mult; flat; #it{p}_{T} (GeV/#it{c})", {kTHnSparseF, {multAxis, flatAxis, ptAxis}}}); + registryMC.add({fmt::format(CpTmcClosureRecF.data(), CspeciesAll[i].data()).c_str(), "Gen Nch w/ Nrec > 0 + Evt. sel; Gen Nch (|#eta|<0.8); flat; #it{p}_{T} (GeV/#it{c})", {kTHnSparseF, {nChAxis, flatAxis, ptAxis}}}); + } + LOG(info) << "processMcEfficiency -> Size of the MC histograms:"; + registryMC.print(); + } + + if (doprocessMcSgnEvtLoss) { registryMC.add("Events/hFlatMCGen", "Events/hFlatMCGen", {kTH1F, {flatAxis}}); registryMC.add("Events/hEvtMcGen", "Events/hEvtMcGen", {kTH1F, {{4, 0.f, 4.f}}}); registryMC.get(HIST("Events/hEvtMcGen"))->GetXaxis()->SetBinLabel(1, "all"); @@ -817,10 +814,6 @@ struct FlattenictyPikp { registryMC.add("Events/hNchGen", "Gen Nch; Gen Nch (|#eta|<0.8)", {kTH1F, {nChAxis}}); registryMC.add("Events/hNchGenCent", "Gen cent; mult", {kTH1F, {multAxis}}); registryMC.add("Events/hNchGenVsCent", "Gen Nch vs Cent; mult; Gen Nch (|#eta|<0.8)", {kTH2F, {nChAxis, multAxis}}); - registryMC.add("Events/hVtxZRec", "MC Rec vertex z position", kTH1F, {vtxzAxis}); - registryMC.add("Events/hVtxZGen", "Generated vertex z position", kTH1F, {vtxzAxis}); - registryMC.add("Events/hNchTVX", "Nch in FT0A+FT0C; Nch; status", {kTH2F, {nChAxis, {2, 0, 2}}}); - registryMC.add({"Tracks/hPtFakes", "Fake tracks; Gen Nch (|#eta|<0.8); flat; #it{p}_{T} (GeV/#it{c})", {kTHnSparseF, {nChAxis, flatAxis, ptAxis}}}); // Event loss registryMC.add("Events/hNchVsFlatGenINELgt0", "Gen Nch w/o Evt sel; Gen Nch (|#eta|<0.8); flat", {kTH2F, {nChAxis, flatAxis}}); registryMC.add("Events/hNchVsFlatGenINELgt0wRecEvtSel", "Gen Nch w/ Nrec > 0 + Evt sel; Gen Nch (|#eta|<0.8); flat", {kTH2F, {nChAxis, flatAxis}}); @@ -828,72 +821,28 @@ struct FlattenictyPikp { registryMC.add("Events/hFlatResponse", "Flattenicity response; flat ture; flat measured", {kTH2F, {flatAxis, flatAxis}}); // Event split registryMC.add("Events/hCentVsFlatRecINELgt0", "Gen evt w/o Evt sel; mult; flat", {kTH2F, {multAxis, flatAxis}}); - registryMC.add("Events/hCentVsFlatRecINELgt0wRecEvt", "Gen evt w/ Nrec > 0; mult; flat", {kTH2F, {multAxis, flatAxis}}); registryMC.add("Events/hCentVsFlatRecINELgt0wRecEvtSel", "Gen evt w/ Nrec > 0 + Evt sel; mult; flat", {kTH2F, {multAxis, flatAxis}}); - registryMC.add("Tracks/hPtRes", "#it{p}_{T} resolution;;(#it{p}_{T}_{rec} - #it{p}_{T}_{gen})/#it{p}_{T}_{gen};", kTH2F, {ptAxis, {100, -1.0, 1.0}}); + // Signal loss for (int i = 0; i < Npart; ++i) { - // Signal loss registryMC.add({fmt::format(CpTgenPrimSgnF.data(), CspeciesAll[i].data()).c_str(), "Gen evt w/o Evt sel; Gen Nch (|#eta|<0.8); flat; #it{p}_{T} (GeV/#it{c})", {kTHnSparseF, {nChAxis, flatAxis, ptAxis}}}); registryMC.add({fmt::format(CpTrecCollPrimSgnF.data(), CspeciesAll[i].data()).c_str(), "Gen Nch w/ Nrec > 0; Gen Nch (|#eta|<0.8); flat; #it{p}_{T} (GeV/#it{c})", {kTHnSparseF, {nChAxis, flatAxis, ptAxis}}}); - // Closure test registryMC.add({fmt::format(CpTmcClosureGenPrimF.data(), CspeciesAll[i].data()).c_str(), "Gen evt w/o Evt sel; Gen Nch (|#eta|<0.8); flat; #it{p}_{T} (GeV/#it{c})", {kTHnSparseF, {nChAxis, flatAxis, ptAxis}}}); - registryMC.add({fmt::format(CpTmcClosureRecF.data(), CspeciesAll[i].data()).c_str(), "Gen Nch w/ Nrec > 0 + Evt. sel; Gen Nch (|#eta|<0.8); flat; #it{p}_{T} (GeV/#it{c})", {kTHnSparseF, {nChAxis, flatAxis, ptAxis}}}); - } - if (defOpt.fillNclVsPhiCutQaHist) { - registryMC.add("Tracks/postSel/hPtPhi", "; #it{p}_{T} (GeV/#it{c}); fmod(#varphi,#pi/9)", {kTH2F, {ptAxis, phiAxisMod}}); - registryMC.add("Tracks/postSel/hPtPhiNclTPC", "; #eta; #it{p}_{T} (GeV/#it{c}); fmod(#varphi,#pi/9); N_{cluster}", {kTHnSparseF, {etaAxis, ptAxis, phiAxisMod, clTpcAxis}}); - registryMC.add("Tracks/postSel/hPtPhiNclPIDTPC", "; #eta; #it{p}_{T} (GeV/#it{c}); fmod(#varphi,#pi/9); N_{PID cluster}", {kTHnSparseF, {etaAxis, ptAxis, phiAxisMod, clTpcAxis}}); - registryMC.add("Tracks/postSel/hPtNclTPC", "; #it{p}_{T} (GeV/#it{c}); N_{cluster}", {kTH2F, {ptAxis, clTpcAxis}}); - registryMC.add("Tracks/postSel/pPtNclTPC", "; #it{p}_{T} (GeV/#it{c}); N_{cluster}", {kTProfile, {ptAxis}}); - registryMC.add("Tracks/postSel/hPtNclPIDTPC", "; #it{p}_{T} (GeV/#it{c}); N_{PID cluster}", {kTH2F, {ptAxis, clTpcAxis}}); - registryMC.add("Tracks/postSel/pPtNclPIDTPC", "; #it{p}_{T} (GeV/#it{c}); N_{PID cluster}", {kTProfile, {ptAxis}}); - } - registryMC.addClone("Tracks/postSel/", "Tracks/preSel/"); - - for (int i = 0; i < NpartChrg; i++) { - const std::string strID = Form("/%s/%s", (i < Npart) ? "pos" : "neg", pID[i % Npart]); - hPtEffRecGoodCollPrimSelTrk[i] = registryMC.add("Tracks/hPtEffRecGoodCollPrimSelTrk" + strID, " ; mult; flat; #it{p}_{T} (GeV/#it{c})", kTHnSparseF, {multAxis, flatAxis, ptAxis}); - hPtGenRecEvt[i] = registryMC.add("Tracks/hPtGenRecEvt" + strID, "Gen evt w/ Nrec > 0; mult; flat; #it{p}_{T} (GeV/#it{c})", kTHnSparseF, {multAxis, flatAxis, ptAxis}); - hPtGenPrimRecEvt[i] = registryMC.add("Tracks/hPtGenPrimRecEvt" + strID, "Gen evt w/ Nrec > 0 (primary); mult; flat; #it{p}_{T} (GeV/#it{c})", kTHnSparseF, {multAxis, flatAxis, ptAxis}); - hPtGenRecEvtGtZero[i] = registryMC.add("Tracks/hPtGenRecEvtGtZero" + strID, "Gen evt w/ Nrec > 0; mult; flat; #it{p}_{T} (GeV/#it{c})", kTHnSparseF, {multAxis, flatAxis, ptAxis}); - hPtGenPrimRecEvtGtZero[i] = registryMC.add("Tracks/hPtGenPrimRecEvtGtZero" + strID, "Gen evt w/ Nrec > 0 (primary); mult; flat; #it{p}_{T} (GeV/#it{c})", kTHnSparseF, {multAxis, flatAxis, ptAxis}); - hPtEffGenPrim[i] = registryMC.add("Tracks/hPtEffGenPrim" + strID, "Gen evt w/o rec Evt; mult; flat; #it{p}_{T} (GeV/#it{c})", kTHnSparseF, {multAxis, flatAxis, ptAxis}); - hPtEffGenWeak[i] = registryMC.add("Tracks/hPtEffGenWeak" + strID, "Gen evt w/o rec Evt; mult; flat; #it{p}_{T} (GeV/#it{c})", kTHnSparseF, {multAxis, flatAxis, ptAxis}); - hPtEffGenMat[i] = registryMC.add("Tracks/hPtEffGenMat" + strID, "Gen evt w/o rec Evt; mult; flat; #it{p}_{T} (GeV/#it{c})", kTHnSparseF, {multAxis, flatAxis, ptAxis}); - hPtEffGenPrimEvtSelGen[i] = registryMC.add("Tracks/hPtEffGenPrimEvtSelGen" + strID, "Gen evt w/o rec Evt; mult; flat; #it{p}_{T} (GeV/#it{c})", kTHnSparseF, {multAxis, flatAxis, ptAxis}); - hPtEffGenWeakEvtSelGen[i] = registryMC.add("Tracks/hPtEffGenWeakEvtSelGen" + strID, "Gen evt w/o rec Evt; mult; flat; #it{p}_{T} (GeV/#it{c})", kTHnSparseF, {multAxis, flatAxis, ptAxis}); - hPtEffGenMatEvtSelGen[i] = registryMC.add("Tracks/hPtEffGenMatEvtSelGen" + strID, "Gen evt w/o rec Evt; mult; flat; #it{p}_{T} (GeV/#it{c})", kTHnSparseF, {multAxis, flatAxis, ptAxis}); } + LOG(info) << "processMcSgnEvtLoss -> Size of the MC histograms:"; + registryMC.print(); + } + if (doprocessMcRecTrack) { for (int i = 0; i < Npart; i++) { - registryMC.add({fmt::format(CpTeffGenPrimRecEvtF.data(), CspeciesAll[i].data()).c_str(), "Gen evt w/ Nrec > 0; mult; flat; #it{p}_{T} (GeV/#it{c})", {kTHnSparseF, {multAxis, flatAxis, ptAxis}}}); - registryMC.add({fmt::format(CpTeffPrimRecEvtF.data(), CspeciesAll[i].data()).c_str(), "Gen evt w/ Nrec > 0 + Evt sel; mult; flat; #it{p}_{T} (GeV/#it{c})", {kTHnSparseF, {multAxis, flatAxis, ptAxis}}}); + registryMC.add({fmt::format(CdEdxMcRecPrimF.data(), CspeciesAll[i].data()).c_str(), "; #eta; mult; flat; #it{p} (GeV/#it{c}); dEdx", {kTHnSparseF, {etaAxis, multAxis, flatAxis, pAxis, dEdxAxis}}}); registryMC.add({fmt::format(CdEdxMcRecPrimSelF.data(), CspeciesAll[i].data()).c_str(), "; #eta; mult; flat; #it{p} (GeV/#it{c}); dEdx", {kTHnSparseF, {etaAxis, multAxis, flatAxis, pAxis, dEdxAxis}}}); registryMC.add({fmt::format(CEtaVsPtVsPMcRecPrimSelF.data(), CspeciesAll[i].data()).c_str(), "; #eta; #it{p}_{T} (GeV/#it{c}); #it{p} (GeV/#it{c})", {kTHnSparseF, {etaAxis, ptAxis, pAxis}}}); - if (defOpt.fillMCRecDCA) { - registryMC.add({fmt::format(CdEdxMcRecPrimF.data(), CspeciesAll[i].data()).c_str(), "; #eta; mult; flat; #it{p} (GeV/#it{c}); dEdx", {kTHnSparseF, {etaAxis, multAxis, flatAxis, pAxis, dEdxAxis}}}); - registryMC.add({fmt::format(CpTvsDCAxyAllF.data(), CspeciesAll[i].data()).c_str(), "; mult; flat; #it{p}_{T} (GeV/#it{c}); DCA_{xy} (cm)", {kTHnSparseF, {multAxis, flatAxis, ptAxis, dcaXYAxis}}}); - registryMC.add({fmt::format(CpTvsDCAxyPrimAllF.data(), CspeciesAll[i].data()).c_str(), "; mult; flat; #it{p}_{T} (GeV/#it{c}); DCA_{xy} (cm)", {kTHnSparseF, {multAxis, flatAxis, ptAxis, dcaXYAxis}}}); - registryMC.add({fmt::format(CpTvsDCAxyWeakAllF.data(), CspeciesAll[i].data()).c_str(), "; mult; flat; #it{p}_{T} (GeV/#it{c}); DCA_{xy} (cm)", {kTHnSparseF, {multAxis, flatAxis, ptAxis, dcaXYAxis}}}); - registryMC.add({fmt::format(CpTvsDCAxyMatAllF.data(), CspeciesAll[i].data()).c_str(), "; mult; flat; #it{p}_{T} (GeV/#it{c}); DCA_{xy} (cm)", {kTHnSparseF, {multAxis, flatAxis, ptAxis, dcaXYAxis}}}); - } - } - - if (defOpt.fillMCRecCheck) { - for (int i = 0; i < NpartChrg; i++) { - const std::string strID = Form("/%s/%s", (i < Npart) ? "pos" : "neg", pID[i % Npart]); - hPtEffRecGoodCollPrim[i] = registryMC.add("Tracks/hPtEffRecGoodCollPrim" + strID, "; mult; flat; #it{p}_{T} (GeV/#it{c})", kTHnSparseF, {multAxis, flatAxis, ptAxis}); - hPtEffRecGoodCollMat[i] = registryMC.add("Tracks/hPtEffRecGoodCollMat" + strID, "; mult; flat; #it{p}_{T} (GeV/#it{c})", kTHnSparseF, {multAxis, flatAxis, ptAxis}); - hPtEffRecGoodCollWeak[i] = registryMC.add("Tracks/hPtEffRecGoodCollWeak" + strID, "; mult; flat; #it{p}_{T} (GeV/#it{c})", kTHnSparseF, {multAxis, flatAxis, ptAxis}); - hDCAxyRecBadCollPrim[i] = registryMC.add("Tracks/hDCAxyRecBadCollPrim" + strID, "; #it{p}_{T} (GeV/#it{c}); DCA_{xy} (cm)", kTH2F, {ptAxis, dcaXYAxis}); - hDCAxyRecBadCollWeak[i] = registryMC.add("Tracks/hDCAxyRecBadCollWeak" + strID, "; #it{p}_{T} (GeV/#it{c}); DCA_{xy} (cm)", kTH2F, {ptAxis, dcaXYAxis}); - hDCAxyRecBadCollMat[i] = registryMC.add("Tracks/hDCAxyRecBadCollMat" + strID, "; #it{p}_{T} (GeV/#it{c}); DCA_{xy} (cm)", kTH2F, {ptAxis, dcaXYAxis}); - hPtVsDCAxyRecGoodCollPrim[i] = registryMC.add("Tracks/hPtVsDCAxyRecGoodCollPrim" + strID, "; #it{p}_{T} (GeV/#it{c}); DCA_{xy} (cm)", kTH2F, {ptAxis, dcaXYAxis}); - hPtVsDCAxyRecGoodCollWeak[i] = registryMC.add("Tracks/hPtVsDCAxyRecGoodCollWeak" + strID, "; #it{p}_{T} (GeV/#it{c}); DCA_{xy} (cm)", kTH2F, {ptAxis, dcaXYAxis}); - hPtVsDCAxyRecGoodCollMat[i] = registryMC.add("Tracks/hPtVsDCAxyRecGoodCollMat" + strID, "; #it{p}_{T} (GeV/#it{c}); DCA_{xy} (cm)", kTH2F, {ptAxis, dcaXYAxis}); - } - } - LOG(info) << "Size of the MC histograms:"; + registryMC.add({fmt::format(CpTvsDCAxyAllF.data(), CspeciesAll[i].data()).c_str(), "; mult; flat; #it{p}_{T} (GeV/#it{c}); DCA_{xy} (cm)", {kTHnSparseF, {multAxis, flatAxis, ptAxis, dcaXYAxis}}}); + registryMC.add({fmt::format(CpTvsDCAxyPrimAllF.data(), CspeciesAll[i].data()).c_str(), "; mult; flat; #it{p}_{T} (GeV/#it{c}); DCA_{xy} (cm)", {kTHnSparseF, {multAxis, flatAxis, ptAxis, dcaXYAxis}}}); + registryMC.add({fmt::format(CpTvsDCAxyWeakAllF.data(), CspeciesAll[i].data()).c_str(), "; mult; flat; #it{p}_{T} (GeV/#it{c}); DCA_{xy} (cm)", {kTHnSparseF, {multAxis, flatAxis, ptAxis, dcaXYAxis}}}); + registryMC.add({fmt::format(CpTvsDCAxyMatAllF.data(), CspeciesAll[i].data()).c_str(), "; mult; flat; #it{p}_{T} (GeV/#it{c}); DCA_{xy} (cm)", {kTHnSparseF, {multAxis, flatAxis, ptAxis, dcaXYAxis}}}); + } + LOG(info) << "processMcRecTrack -> Size of the MC histograms:"; registryMC.print(); } } @@ -1000,15 +949,6 @@ struct FlattenictyPikp { return fCalibDeDxFunc; } - template - bool isPID(const P& mcParticle) - { - static_assert(pidSgn == CnullInt || pidSgn == ConeInt); - static_assert(id > CnullInt && id < Npart); - constexpr int Cidx = id + pidSgn * Npart; - return mcParticle.pdgCode() == pIdSgn[Cidx]; - } - template bool selTOFPi(T const& track) { @@ -1357,36 +1297,6 @@ struct FlattenictyPikp { return std::abs(charge) >= CminCharge; } - template - bool isInelGt0wTVX(P const& particles) - { - int nChrgMc = 0; - int nChrgFT0A = 0; - int nChrgFT0C = 0; - for (auto const& particle : particles) { - if (!isChrgParticle(particle.pdgCode())) { - continue; - } - if (!particle.isPhysicalPrimary()) { - continue; - } - // trigger TVX - if (particle.eta() > CminAccFT0A && particle.eta() < CmaxAccFT0A) { - nChrgFT0A++; - } - if (particle.eta() > CminAccFT0C && particle.eta() < CmaxAccFT0C) { - nChrgFT0C++; - } - nChrgMc++; - } - if (nChrgFT0A == CnullInt || nChrgFT0C == CnullInt) { - registryMC.fill(HIST("Events/hNchTVX"), nChrgMc, 0.5); - return false; - } - registryMC.fill(HIST("Events/hNchTVX"), nChrgMc, 1.5); - return true; - } - template int countTracks(T const& tracks, C const& collision, aod::BCsWithTimestamps const& /*bcs*/, float mult) { @@ -2260,128 +2170,130 @@ struct FlattenictyPikp { return flatFV0; } - template - void fillMCRecTrack(MyLabeledPIDTracks::iterator const& track, const float mult, const float flat) + Preslice perCollTrk = aod::track::collisionId; + + void processMcEfficiency(CollsMCExtraMult::iterator const& mcCollision, + soa::SmallGroups const& collisions, + aod::BCsWithTimestamps const& /*bcs*/, + aod::FV0As const& /*fv0s*/, + aod::McParticles const& particles, + MyLabeledPIDTracks const& tracks) { - static_assert(pidSgn == CnullInt || pidSgn == ConeInt); - constexpr int ChistIdx = id + pidSgn * Npart; - // LOG(debug) << "fillMCRecTrack for pidSgn '" << pidSgn << "' and id '" << static_cast(id) << " with index " << ChistIdx; - const aod::McParticles::iterator& mcParticle = track.mcParticle(); - const CollsGen::iterator& collision = track.collision_as>(); - // const CollsGen::iterator& collision = track.collision_as(); - - if (!isChrgParticle(mcParticle.pdgCode())) { - return; - } - if (std::abs(mcParticle.eta()) > trkSelOpt.trkEtaMax) { - return; + LOGP(debug, "MC col {} has {} reco cols", mcCollision.globalIndex(), collisions.size()); + auto multMC = -999.; + if (evtSelOpt.useMultMCmidrap) { + multMC = mcCollision.multMCNParticlesEta08(); + } else { + multMC = getMultMC(mcCollision); } - if (mcParticle.pt() < trkSelOpt.trkPtMin) { - return; + const float flatMC = fillFlatMC(particles); + + auto multRec = -999.; + auto flatRec = -999.; + bool gtOneRec = false; + for (const auto& collision : collisions) { + if (!isGoodEvent(collision)) { + continue; + } + if (evtSelOpt.removeSplitVertex && collision.globalIndex() != mcCollision.bestCollisionIndex()) { + continue; + } + gtOneRec = true; + multRec = getMult(collision); + flatRec = fillFlat(collision); } - if (!isPID(mcParticle)) { - return; + registryMC.fill(HIST("Events/hVtxZGen"), mcCollision.posZ()); + if (gtOneRec) { + registryMC.fill(HIST("Events/hVtxZGenGtOneRec"), mcCollision.posZ()); } - if (mcParticle.isPhysicalPrimary()) { - if (isDCAxyCut(track)) { - hPtEffRecGoodCollPrimSelTrk[ChistIdx]->Fill(mult, flat, mcParticle.pt()); + /* Tracking eff. den */ + for (const auto& particle : particles) { + if (!isChrgParticle(particle.pdgCode())) { + continue; } - } - if (defOpt.fillMCRecCheck) { // for checking purposes only: use gen Nch, gen Flat - if ((collision.has_mcCollision() && (mcParticle.mcCollisionId() != collision.mcCollisionId())) || !collision.has_mcCollision()) { - if (!mcParticle.isPhysicalPrimary()) { - if (mcParticle.getProcess() == CprocessIdWeak) { - hDCAxyRecBadCollWeak[ChistIdx]->Fill(track.pt(), track.dcaXY()); - } else { - hDCAxyRecBadCollMat[ChistIdx]->Fill(track.pt(), track.dcaXY()); - } - } else { - hDCAxyRecBadCollPrim[ChistIdx]->Fill(track.pt(), track.dcaXY()); - } + if (!particle.isPhysicalPrimary()) { + continue; } - if (collision.has_mcCollision() && (mcParticle.mcCollisionId() == collision.mcCollisionId())) { - if (!mcParticle.isPhysicalPrimary()) { - if (mcParticle.getProcess() == CprocessIdWeak) { - hPtEffRecGoodCollWeak[ChistIdx]->Fill(mult, flat, track.pt()); - hPtVsDCAxyRecGoodCollWeak[ChistIdx]->Fill(track.pt(), track.dcaXY()); - } else { - hPtEffRecGoodCollMat[ChistIdx]->Fill(mult, flat, track.pt()); - hPtVsDCAxyRecGoodCollMat[ChistIdx]->Fill(track.pt(), track.dcaXY()); + if (std::abs(particle.eta()) > trkSelOpt.trkEtaMax) { + continue; + } + if (gtOneRec) { + static_for<0, 4>([&](auto i) { + constexpr int Cidx = i.value; + if (std::fabs(particle.pdgCode()) == pDGs[Cidx]) { + registryMC.fill(HIST(Cprefix) + HIST(CspeciesAll[Cidx]) + HIST(CpTeffGenPrimRecEvt), multRec, flatRec, particle.pt()); } - } else { - hPtEffRecGoodCollPrim[ChistIdx]->Fill(mult, flat, track.pt()); - hPtVsDCAxyRecGoodCollPrim[ChistIdx]->Fill(track.pt(), track.dcaXY()); - } + }); } } - } - - template - void fillMCGenRecEvt(aod::McParticles::iterator const& mcParticle, const float mult, const float flat) - { - static_assert(pidSgn == CnullInt || pidSgn == ConeInt); - constexpr int ChistIdx = id + pidSgn * Npart; - - if (!isPID(mcParticle)) { - return; - } - if constexpr (isGtZeroColl) { - hPtGenRecEvtGtZero[ChistIdx]->Fill(mult, flat, mcParticle.pt()); - if (mcParticle.isPhysicalPrimary()) { - hPtGenPrimRecEvtGtZero[ChistIdx]->Fill(mult, flat, mcParticle.pt()); + /* Tracking eff. num */ + for (const auto& collision : collisions) { + if (trkSelOpt.rejectTrkAtTPCSector || defOpt.applyCalibGainFromCCDB || defOpt.applyCalibVtxFromCCDB) { + auto bc = collision.bc_as(); + int currentRun = bc.runNumber(); + if (runNumber != currentRun) { + initCCDB(bc); + runNumber = currentRun; + } } - } else { - hPtGenRecEvt[ChistIdx]->Fill(mult, flat, mcParticle.pt()); - if (mcParticle.isPhysicalPrimary()) { - hPtGenPrimRecEvt[ChistIdx]->Fill(mult, flat, mcParticle.pt()); + registryMC.fill(HIST("Events/hMCEventStatus"), static_cast(MCEventStatus::kMcEvtAll), multRec, flatRec); + if (!isGoodEvent(collision)) { + continue; } - } - } - - template - void fillMCGen(aod::McParticles::iterator const& mcParticle, const float mult, const float flat) - { - static_assert(pidSgn == CnullInt || pidSgn == ConeInt); - constexpr int ChistIdx = id + pidSgn * Npart; - - if (!isPID(mcParticle)) { - return; - } + registryMC.fill(HIST("Events/hMCEventStatus"), static_cast(MCEventStatus::kEvtSelected), multRec, flatRec); + if (!collision.has_mcCollision()) { + continue; + } + registryMC.fill(HIST("Events/hMCEventStatus"), static_cast(MCEventStatus::kEvtHasMcColl), multRec, flatRec); + if (evtSelOpt.removeSplitVertex && collision.globalIndex() != mcCollision.bestCollisionIndex()) { + continue; + } + registryMC.fill(HIST("Events/hMCEventStatus"), static_cast(MCEventStatus::kEvtSplitVtxRemoved), multRec, flatRec); + registryMC.fill(HIST("Events/hVtxZRec"), collision.posZ()); - if constexpr (evtSel) { - if (!mcParticle.isPhysicalPrimary()) { - if (mcParticle.getProcess() == CprocessIdWeak) { - hPtEffGenWeakEvtSelGen[ChistIdx]->Fill(mult, flat, mcParticle.pt()); - } else { - hPtEffGenMatEvtSelGen[ChistIdx]->Fill(mult, flat, mcParticle.pt()); + const auto& groupedTrks = tracks.sliceBy(perCollTrk, collision.globalIndex()); + int nTrk = 0; + for (const auto& track : groupedTrks) { + if (!isGoodTrack(track, magField)) { + continue; } - } else { - hPtEffGenPrimEvtSelGen[ChistIdx]->Fill(mult, flat, mcParticle.pt()); - } - } else { - if (!mcParticle.isPhysicalPrimary()) { - if (mcParticle.getProcess() == CprocessIdWeak) { - hPtEffGenWeak[ChistIdx]->Fill(mult, flat, mcParticle.pt()); - } else { - hPtEffGenMat[ChistIdx]->Fill(mult, flat, mcParticle.pt()); + if (!track.has_mcParticle()) { + registryMC.fill(HIST("Tracks/hPtFakes"), multRec, flatRec, track.pt()); + continue; } - } else { - hPtEffGenPrim[ChistIdx]->Fill(mult, flat, mcParticle.pt()); + const auto& particle = track.mcParticle_as(); + if (collision.mcCollisionId() != particle.mcCollisionId()) { + continue; + } + if (!isChrgParticle(particle.pdgCode())) { + continue; + } + if (!particle.isPhysicalPrimary()) { + continue; + } + static_for<0, 4>([&](auto i) { + constexpr int Cidx = i.value; + if (std::fabs(particle.pdgCode()) == pDGs[Cidx]) { + registryMC.fill(HIST(Cprefix) + HIST(CspeciesAll[Cidx]) + HIST(CpTeffPrimRecEvt), multRec, flatRec, track.pt()); + registryMC.fill(HIST(Cprefix) + HIST(CspeciesAll[Cidx]) + HIST(CpTmcClosureRec), multMC, flatMC, track.pt()); // closure + } + }); + nTrk++; + registryMC.fill(HIST("Tracks/hPtRes"), particle.pt(), (track.pt() - particle.pt()) / particle.pt()); } + registryQC.fill(HIST("Events/hNchVsCent"), nTrk, multRec); } } + PROCESS_SWITCH(FlattenictyPikp, processMcEfficiency, "process tracking efficiency", false); - Preslice perCollTrk = aod::track::collisionId; - - void processMC(CollsMCExtraMult::iterator const& mcCollision, - soa::SmallGroups const& collisions, - aod::BCsWithTimestamps const& /*bcs*/, - aod::FV0As const& /*fv0s*/, - aod::McParticles const& particles, - MyLabeledPIDTracks const& tracks) + void processMcSgnEvtLoss(CollsMCExtraMult::iterator const& mcCollision, + soa::SmallGroups const& collisions, + aod::BCsWithTimestamps const& /*bcs*/, + aod::FV0As const& /*fv0s*/, + aod::McParticles const& particles) { LOGP(debug, "MC col {} has {} reco cols", mcCollision.globalIndex(), collisions.size()); - auto multMC = -1.; + auto multMC = -999; if (evtSelOpt.useMultMCmidrap) { multMC = mcCollision.multMCNParticlesEta08(); } else { @@ -2405,103 +2317,129 @@ struct FlattenictyPikp { } } */ - registryMC.fill(HIST("Events/hMCEventStatus"), static_cast(MCEventStatus::kMcEvtAll)); registryMC.fill(HIST("Events/hNchGen"), multMC); registryMC.fill(HIST("Events/hNchGenCent"), getGenCent(mcCollision)); - - const float flatMC = fillFlatMC(particles); + const float flatMC = fillFlatMC(particles); registryMC.fill(HIST("Events/hFlatMCGen"), flatMC); - // Loop on rec collisions - // Obtain here: Numerator of tracking efficiency; Secondary contamination correction - for (const auto& collision : collisions) { - if (trkSelOpt.rejectTrkAtTPCSector || defOpt.applyCalibGainFromCCDB || defOpt.applyCalibVtxFromCCDB) { - auto bc = collision.bc_as(); - int currentRun = bc.runNumber(); - if (runNumber != currentRun) { - initCCDB(bc); - runNumber = currentRun; - } + /* Signal/event loss den */ + registryMC.fill(HIST("Events/hEvtMcGen"), 0.5); + if (evtSelOpt.zVtxCutMC && std::abs(mcCollision.posZ()) > evtSelOpt.cutVtxZ) { + return; + } + registryMC.fill(HIST("Events/hEvtMcGen"), 1.5); + if (evtSelOpt.useINELCutMC) { + if (!o2::pwglf::isINELgt0mc(particles, pdg)) { + return; } - registryMC.fill(HIST("Events/hCentVsFlatRecINELgt0"), getMult(collision), fillFlat(collision)); // Evt split den - if (evtSelOpt.removeSplitVertex && collision.globalIndex() != mcCollision.bestCollisionIndex()) { + } + registryMC.fill(HIST("Events/hEvtMcGen"), 2.5); + if (evtSelOpt.useInelgt0wTVX && !(mcCollision.multMCFT0C() <= 0 || mcCollision.multMCFT0A() <= 0)) { + return; + } + registryMC.fill(HIST("Events/hEvtMcGen"), 3.5); + registryMC.fill(HIST("Events/hNchVsFlatGenINELgt0"), multMC, flatMC); // Evt loss den + + for (const auto& particle : particles) { + if (!isChrgParticle(particle.pdgCode())) { continue; } - registryMC.fill(HIST("Events/hMCEventStatus"), static_cast(MCEventStatus::kMcEvtSplit)); - registryMC.fill(HIST("Events/hCentVsFlatRecINELgt0wRecEvt"), getMult(collision), fillFlat(collision)); // Evt split num, w/ Nrec > 0 - - for (const auto& particle : particles) { - if (!isChrgParticle(particle.pdgCode())) { - continue; - } - if (!particle.isPhysicalPrimary()) { - continue; - } - if (std::abs(particle.eta()) > trkSelOpt.trkEtaMax) { - continue; - } - if (particle.pt() < trkSelOpt.trkPtMin) { - continue; - } - static_for<0, 1>([&](auto pidSgn) { - fillMCGenRecEvt(particle, multMC, flatMC); - fillMCGenRecEvt(particle, multMC, flatMC); - fillMCGenRecEvt(particle, multMC, flatMC); - }); + if (!particle.isPhysicalPrimary()) { + continue; + } + if (std::abs(particle.eta()) > trkSelOpt.trkEtaMax) { + continue; } + static_for<0, 4>([&](auto i) { + constexpr int Cidx = i.value; + if (std::fabs(particle.pdgCode()) == pDGs[Cidx]) { + registryMC.fill(HIST(Cprefix) + HIST(CspeciesAll[Cidx]) + HIST(CpTgenPrimSgn), multMC, flatMC, particle.pt()); // Sgn loss den + registryMC.fill(HIST(Cprefix) + HIST(CspeciesAll[Cidx]) + HIST(CpTmcClosureGenPrim), multMC, flatMC, particle.pt()); // closure + } + }); + } + + /* Signal/event loss num */ + auto multRec = -999.; + auto flatRec = -999.; + bool gtOneRec = false; + for (const auto& collision : collisions) { if (!isGoodEvent(collision)) { continue; } - const float multRecGt1 = getMult(collision); - const float flatRec = fillFlat(collision); + if (evtSelOpt.removeSplitVertex && collision.globalIndex() != mcCollision.bestCollisionIndex()) { + continue; + } + gtOneRec = true; + multRec = getMult(collision); + flatRec = fillFlat(collision); + } - registryMC.fill(HIST("Events/hVtxZRec"), collision.posZ()); - registryMC.fill(HIST("Events/hCentVsFlatRecINELgt0wRecEvtSel"), multRecGt1, flatRec); // Evt split num, w/ Nrec > 0 + Evt. sel - registryMC.fill(HIST("Events/hNchGenVsCent"), multMC, multRecGt1); - registryMC.fill(HIST("Events/hNchVsFlatGenINELgt0wRecEvtSel"), multMC, flatMC); // Evt loss num, w/ Nrec > 0 + Evt. sel + registryMC.fill(HIST("Events/hCentVsFlatRecINELgt0"), multRec, flatRec); // Evt split den + + if (gtOneRec) { + registryMC.fill(HIST("Events/hCentVsFlatRecINELgt0wRecEvtSel"), multRec, flatRec); // Evt split num + registryMC.fill(HIST("Events/hNchGenVsCent"), multMC, multRec); + registryMC.fill(HIST("Events/hNchVsFlatGenINELgt0wRecEvtSel"), multMC, flatMC); // Evt loss num registryMC.fill(HIST("Events/hFlatResponse"), flatMC, flatRec); + } - // Obtain here: Denominator of tracking efficiency; Numerator event and signal loss - for (const auto& particle : particles) { - if (!isChrgParticle(particle.pdgCode())) { - continue; - } - if (!particle.isPhysicalPrimary()) { - continue; - } - if (std::abs(particle.eta()) > trkSelOpt.trkEtaMax) { - continue; - } - if (particle.pt() < trkSelOpt.trkPtMin) { - continue; - } - static_for<0, 1>([&](auto pidSgn) { - fillMCGenRecEvt(particle, multMC, flatMC); - fillMCGenRecEvt(particle, multMC, flatMC); - fillMCGenRecEvt(particle, multMC, flatMC); - }); + for (const auto& particle : particles) { + if (!isChrgParticle(particle.pdgCode())) { + continue; + } + if (!particle.isPhysicalPrimary()) { + continue; + } + if (std::abs(particle.eta()) > trkSelOpt.trkEtaMax) { + continue; + } + if (gtOneRec) { static_for<0, 4>([&](auto i) { constexpr int Cidx = i.value; if (std::fabs(particle.pdgCode()) == pDGs[Cidx]) { - registryMC.fill(HIST(Cprefix) + HIST(CspeciesAll[Cidx]) + HIST(CpTrecCollPrimSgn), multMC, flatMC, particle.pt()); // Sgn loss num - registryMC.fill(HIST(Cprefix) + HIST(CspeciesAll[Cidx]) + HIST(CpTeffGenPrimRecEvt), multRecGt1, flatRec, particle.pt()); // Tracking eff. den + registryMC.fill(HIST(Cprefix) + HIST(CspeciesAll[Cidx]) + HIST(CpTrecCollPrimSgn), multMC, flatMC, particle.pt()); // Sgn loss num } }); } + } + } + PROCESS_SWITCH(FlattenictyPikp, processMcSgnEvtLoss, "process signal/event loss", false); + + void processMcRecTrack(CollsMCExtraMult::iterator const& mcCollision, + soa::SmallGroups const& collisions, + aod::BCsWithTimestamps const& /*bcs*/, + aod::FV0As const& /*fv0s*/, + aod::McParticles const& /*particles*/, + MyLabeledPIDTracks const& tracks) + { + for (const auto& collision : collisions) { + if (trkSelOpt.rejectTrkAtTPCSector || defOpt.applyCalibGainFromCCDB || defOpt.applyCalibVtxFromCCDB) { + auto bc = collision.bc_as(); + int currentRun = bc.runNumber(); + if (runNumber != currentRun) { + initCCDB(bc); + runNumber = currentRun; + } + } + if (evtSelOpt.removeSplitVertex && collision.globalIndex() != mcCollision.bestCollisionIndex()) { + continue; + } + if (!isGoodEvent(collision)) { + continue; + } + const float multRec = getMult(collision); + const float flatRec = fillFlat(collision); - // Rec tracks; track selection w/ DCA open (for secondaries), w/ DCA close (for efficiency) - // Obtain here: DCAxy for sec contamination, MC closure const auto& groupedTrks = tracks.sliceBy(perCollTrk, collision.globalIndex()); - int nTrk = 0; for (const auto& track : groupedTrks) { if (!track.has_collision()) { continue; } - if (!isGoodTrack(track, magField)) { + if (!isGoodTrack(track, magField)) { continue; } if (!track.has_mcParticle()) { - registryMC.fill(HIST("Tracks/hPtFakes"), multMC, flatMC, track.pt()); continue; } auto particle = track.mcParticle_as(); @@ -2514,118 +2452,37 @@ struct FlattenictyPikp { if (std::abs(particle.eta()) > trkSelOpt.trkEtaMax) { continue; } - if (particle.pt() < trkSelOpt.trkPtMin) { - continue; - } - static_for<0, 1>([&](auto pidSgn) { - fillMCRecTrack(track, multMC, flatMC); - fillMCRecTrack(track, multMC, flatMC); - fillMCRecTrack(track, multMC, flatMC); - }); - if (defOpt.fillMCRecDCA) { - static_for<0, 4>([&](auto i) { - constexpr int Cidx = i.value; - if (std::fabs(particle.pdgCode()) == pDGs[Cidx]) { - if (!particle.isPhysicalPrimary()) { - if (particle.getProcess() == CprocessIdWeak) { - registryMC.fill(HIST(Cprefix) + HIST(CspeciesAll[Cidx]) + HIST(CpTvsDCAxyWeakAll), multRecGt1, flatRec, track.pt(), track.dcaXY()); - } else { - registryMC.fill(HIST(Cprefix) + HIST(CspeciesAll[Cidx]) + HIST(CpTvsDCAxyMatAll), multRecGt1, flatRec, track.pt(), track.dcaXY()); - } + static_for<0, 4>([&](auto i) { + constexpr int Cidx = i.value; + if (std::fabs(particle.pdgCode()) == pDGs[Cidx]) { + if (!particle.isPhysicalPrimary()) { + if (particle.getProcess() == CprocessIdWeak) { + registryMC.fill(HIST(Cprefix) + HIST(CspeciesAll[Cidx]) + HIST(CpTvsDCAxyWeakAll), multRec, flatRec, track.pt(), track.dcaXY()); } else { - registryMC.fill(HIST(Cprefix) + HIST(CspeciesAll[Cidx]) + HIST(CpTvsDCAxyPrimAll), multRecGt1, flatRec, track.pt(), track.dcaXY()); - registryMC.fill(HIST(Cprefix) + HIST(CspeciesAll[Cidx]) + HIST(CdEdxMcRecPrim), track.eta(), multRecGt1, flatRec, track.p(), track.tpcSignal()); + registryMC.fill(HIST(Cprefix) + HIST(CspeciesAll[Cidx]) + HIST(CpTvsDCAxyMatAll), multRec, flatRec, track.pt(), track.dcaXY()); } - registryMC.fill(HIST(Cprefix) + HIST(CspeciesAll[Cidx]) + HIST(CpTvsDCAxyAll), multRecGt1, flatRec, track.pt(), track.dcaXY()); + } else { + registryMC.fill(HIST(Cprefix) + HIST(CspeciesAll[Cidx]) + HIST(CpTvsDCAxyPrimAll), multRec, flatRec, track.pt(), track.dcaXY()); + registryMC.fill(HIST(Cprefix) + HIST(CspeciesAll[Cidx]) + HIST(CdEdxMcRecPrim), track.eta(), multRec, flatRec, track.p(), track.tpcSignal()); } - }); - } - if (isGoodTrack(track, magField)) { + registryMC.fill(HIST(Cprefix) + HIST(CspeciesAll[Cidx]) + HIST(CpTvsDCAxyAll), multRec, flatRec, track.pt(), track.dcaXY()); + } + }); + if (isGoodTrack(track, magField)) { static_for<0, 4>([&](auto i) { constexpr int Cidx = i.value; if (std::fabs(particle.pdgCode()) == pDGs[Cidx]) { if (particle.isPhysicalPrimary()) { - registryMC.fill(HIST(Cprefix) + HIST(CspeciesAll[Cidx]) + HIST(CdEdxMcRecPrimSel), track.eta(), multRecGt1, flatRec, track.p(), track.tpcSignal()); + registryMC.fill(HIST(Cprefix) + HIST(CspeciesAll[Cidx]) + HIST(CdEdxMcRecPrimSel), track.eta(), multRec, flatRec, track.p(), track.tpcSignal()); registryMC.fill(HIST(Cprefix) + HIST(CspeciesAll[Cidx]) + HIST(CEtaVsPtVsPMcRecPrimSel), track.eta(), track.pt(), track.p()); - registryMC.fill(HIST(Cprefix) + HIST(CspeciesAll[Cidx]) + HIST(CpTeffPrimRecEvt), multRecGt1, flatRec, track.pt()); // Tracking eff. num - registryMC.fill(HIST(Cprefix) + HIST(CspeciesAll[Cidx]) + HIST(CpTmcClosureRec), multMC, flatMC, track.pt()); // closure } } }); - nTrk++; } - registryMC.fill(HIST("Tracks/hPtRes"), particle.pt(), (track.pt() - particle.pt()) / particle.pt()); - } - registryQC.fill(HIST("Events/hNchVsCent"), nTrk, multRecGt1); - } - - // Loop on generated particles (no requirement on availaability of reconstructed collision; no event selection) - // - for (const auto& particle : particles) { - if (!isChrgParticle(particle.pdgCode())) { - continue; } - if (std::abs(particle.eta()) > trkSelOpt.trkEtaMax) { - continue; - } - if (particle.pt() < trkSelOpt.trkPtMin) { - continue; - } - static_for<0, 1>([&](auto pidSgn) { - fillMCGen(particle, multMC, flatMC); - fillMCGen(particle, multMC, flatMC); - fillMCGen(particle, multMC, flatMC); - }); - } - - // Obtain here: Denominator of signal loss and event loss; MC closure - // - registryMC.fill(HIST("Events/hEvtMcGen"), 0.5); - if (evtSelOpt.zVtxCutMC && std::abs(mcCollision.posZ()) > evtSelOpt.cutVtxZ) { - return; - } - registryMC.fill(HIST("Events/hVtxZGen"), mcCollision.posZ()); - registryMC.fill(HIST("Events/hEvtMcGen"), 1.5); - if (evtSelOpt.useINELCutMC) { - if (!o2::pwglf::isINELgt0mc(particles, pdg)) { - return; - } - } - registryMC.fill(HIST("Events/hEvtMcGen"), 2.5); - if (evtSelOpt.useInelgt0wTVX && !isInelGt0wTVX(particles)) { // TVX trigger: FT0A + FT0C acceptance - return; - } - registryMC.fill(HIST("Events/hEvtMcGen"), 3.5); - registryMC.fill(HIST("Events/hNchVsFlatGenINELgt0"), multMC, flatMC); // Evt loss den - - for (const auto& particle : particles) { - if (!isChrgParticle(particle.pdgCode())) { - continue; - } - if (std::abs(particle.eta()) > trkSelOpt.trkEtaMax) { - continue; - } - if (particle.pt() < trkSelOpt.trkPtMin) { - continue; - } - static_for<0, 1>([&](auto pidSgn) { - fillMCGen(particle, multMC, flatMC); - fillMCGen(particle, multMC, flatMC); - fillMCGen(particle, multMC, flatMC); - }); - static_for<0, 4>([&](auto i) { - constexpr int Cidx = i.value; - // LOG(debug) << "fillMCGen for pidSgn '" << pidSgn << "' and id '" << static_cast(id) << " with index " << ChistIdx; - if (particle.isPhysicalPrimary()) { - if (std::fabs(particle.pdgCode()) == pDGs[Cidx]) { - registryMC.fill(HIST(Cprefix) + HIST(CspeciesAll[Cidx]) + HIST(CpTgenPrimSgn), multMC, flatMC, particle.pt()); // Sgn loss den - registryMC.fill(HIST(Cprefix) + HIST(CspeciesAll[Cidx]) + HIST(CpTmcClosureGenPrim), multMC, flatMC, particle.pt()); // closure - } - } - }); } } - PROCESS_SWITCH(FlattenictyPikp, processMC, "process MC", false); + PROCESS_SWITCH(FlattenictyPikp, processMcRecTrack, "process MC rec tracks", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)