diff --git a/PWGCF/EbyEFluctuations/Tasks/meanPtFlucId.cxx b/PWGCF/EbyEFluctuations/Tasks/meanPtFlucId.cxx index 3605e90661d..9af2293e4e3 100644 --- a/PWGCF/EbyEFluctuations/Tasks/meanPtFlucId.cxx +++ b/PWGCF/EbyEFluctuations/Tasks/meanPtFlucId.cxx @@ -61,51 +61,71 @@ struct MeanPtFlucId { Configurable nPBins{"nPBins", 300, ""}; Configurable nPhiBins{"nPhiBins", 100, ""}; Configurable nPtBins{"nPtBins", 50, "Number of pT bins"}; - Configurable cfgCutPtMax{"cfgCutPtMax", 2.0, "maximum pT"}; - Configurable cfgCutPtMin{"cfgCutPtMin", 0.2, "minimum pT"}; - Configurable cfgCutEta{"cfgCutEta", 0.8, "Eta cut"}; - Configurable cfgDcaPtDepOffset{"cfgDcaPtDepOffset", 0.004, "Offset of the pT-dependent DCA cut (cm)"}; - Configurable cfgDcaPtDepScale{"cfgDcaPtDepScale", 0.013, "1/pT coefficient of the pT-dependent DCA cut (cm GeV/c)"}; - Configurable cfgCutPosZ{"cfgCutPosZ", 7.0, "cut for vertex Z"}; - Configurable cfgPosZ{"cfgPosZ", true, "Position Z"}; - Configurable cfgCutNSig2{"cfgCutNSig2", 2.0, "nSigma cut: 2"}; - Configurable cfgCutNSig3{"cfgCutNSig3", 3.0, "nSigma cut: 3"}; - Configurable cfgCutNSig5{"cfgCutNSig5", 5.0, "nSigma cut: 5"}; - Configurable cfgCutNSigTpcTof{"cfgCutNSigTpcTof", 3.0, "nSigma TPC cut when TPC + TOF "}; - Configurable cfgSelCutNSigTpcPi{"cfgSelCutNSigTpcPi", 3.0, "nSigma TPC cut for pion selection"}; - Configurable cfgSelCutNSigTpcKa{"cfgSelCutNSigTpcKa", 3.0, "nSigma TPC cut for kaon selection"}; - Configurable cfgSelCutNSigTpcPr{"cfgSelCutNSigTpcPr", 3.0, "nSigma TPC cut for proton selection"}; - Configurable cfgSelCutNSigTofPi{"cfgSelCutNSigTofPi", 3.0, "nSigma TOF cut for pion selection"}; - Configurable cfgSelCutNSigTofKa{"cfgSelCutNSigTofKa", 3.0, "nSigma TOF cut for kaon selection"}; - Configurable cfgSelCutNSigTofPr{"cfgSelCutNSigTofPr", 3.0, "nSigma TOF cut for proton selection"}; - // Signed TOF nSigma offsets: corrected = raw - shift. Zero preserves the original PID. - Configurable cfgTofShiftPi{"cfgTofShiftPi", 0.f, "Pion TOF nSigma offset to subtract (not a time shift in ps)"}; - Configurable cfgTofShiftKa{"cfgTofShiftKa", 0.f, "Kaon TOF nSigma offset to subtract (not a time shift in ps)"}; - Configurable cfgTofShiftPr{"cfgTofShiftPr", 0.f, "Proton TOF nSigma offset to subtract (not a time shift in ps)"}; - Configurable cfgCutPiPtMin{"cfgCutPiPtMin", 0.2, "Minimum pion p_{T} cut"}; - Configurable cfgCutKaPtMin{"cfgCutKaPtMin", 0.3, "Minimum kaon p_{T} cut"}; - Configurable cfgCutPrPtMin{"cfgCutPrPtMin", 0.5, "Minimum proton p_{T} cut"}; - Configurable cfgCutPiThrsldP{"cfgCutPiThrsldP", 0.7, "Momentum threshold for requiring pion TOF PID"}; - Configurable cfgCutKaThrsldP{"cfgCutKaThrsldP", 0.55, "Momentum threshold for requiring kaon TOF PID"}; - Configurable cfgCutPrThrsldP{"cfgCutPrThrsldP", 0.8, "Momentum threshold for requiring proton TOF PID"}; - Configurable cfgP0U{"cfgP0U", 6.36269, "p_{0} for upper cut for N_{TPC} vs N{sim}"}; - Configurable cfgP1U{"cfgP1U", 2.22108, "p_{1} for upper cut for N_{TPC} vs N{sim}"}; - Configurable cfgP2U{"cfgP2U", -0.0118451, "p_{2} for upper cut for N_{TPC} vs N{sim}"}; - Configurable cfgP0L{"cfgP0L", 23.8213, "p_{0} for lower cut for N_{TPC} vs N{sim}"}; - Configurable cfgP1L{"cfgP1L", -0.402328, "p_{1} for lower cut for N_{TPC} vs N{sim}"}; - Configurable cfgP2L{"cfgP2L", 0.0128423, "p_{2} for lower cut for N_{TPC} vs N{sim}"}; + struct : ConfigurableGroup { + Configurable cfgCutPtMax{"cfgCutPtMax", 2.0, "maximum pT"}; + Configurable cfgCutPtMin{"cfgCutPtMin", 0.2, "minimum pT"}; + Configurable cfgCutEta{"cfgCutEta", 0.8, "Eta cut"}; + Configurable cfgCutITSChi2NCl{"cfgCutITSChi2NCl", 36.0, "ITS Chi2/NCl cut"}; + Configurable cfgCutTPCChi2NCl{"cfgCutTPCChi2NCl", 4.0, "TPC Chi2/NCl cut"}; + Configurable cfgCutTPCNCls{"cfgCutTPCNCls", 70.0, "TPC NCls cut"}; + Configurable cfgCutDCAxy{"cfgCutDCAxy", 0.1, "DCAxy cut"}; + Configurable cfgCutDCAz{"cfgCutDCAz", 0.2, "DCAz cut"}; + Configurable cfgDcaPtDepOffset{"cfgDcaPtDepOffset", 0.004, "Offset of the pT-dependent DCA cut (cm)"}; + Configurable cfgDcaPtDepScale{"cfgDcaPtDepScale", 0.013, "1/pT coefficient of the pT-dependent DCA cut (cm GeV/c)"}; + Configurable cfgCutNCrossedRowsTPC{"cfgCutNCrossedRowsTPC", 70, "Minimum number of crossed TPC rows"}; + Configurable cfgGlobalTrack{"cfgGlobalTrack", true, "Require the standard global-track selection (includes its DCA cuts)"}; + Configurable cfgGlobalTrackWoDCA{"cfgGlobalTrackWoDCA", false, "Require global-track quality without built-in DCA cuts; disable cfgGlobalTrack to use only this selection"}; + Configurable cfgDcaXYPtDep{"cfgDcaXYPtDep", true, "Apply the pT-dependent absolute DCAxy cut; takes precedence over cfgDcaXY"}; + Configurable cfgNTPCCrossedRows{"cfgNTPCCrossedRows", false, "Apply the minimum number of crossed TPC rows"}; + Configurable cfgTPCChi2NCl{"cfgTPCChi2NCl", false, "Chi2 per N_{TPC} Clusters Cut "}; + Configurable cfgITSChi2NCl{"cfgITSChi2NCl", false, "Chi2 per N_{ITS} Clusters Cut "}; + Configurable cfgDcaXY{"cfgDcaXY", false, "DCA XY Cut "}; + Configurable cfgDcaZ{"cfgDcaZ", false, "DCA Z Cut "}; + Configurable cfgDcaZPtDep{"cfgDcaZPtDep", false, "DCA Z pT dependent Cut "}; + } trackCuts; + struct : ConfigurableGroup { + Configurable cfgCutPosZ{"cfgCutPosZ", 7.0, "cut for vertex Z"}; + Configurable cfgP0U{"cfgP0U", 6.36269, "p_{0} for upper cut for N_{TPC} vs N{sim}"}; + Configurable cfgP1U{"cfgP1U", 2.22108, "p_{1} for upper cut for N_{TPC} vs N{sim}"}; + Configurable cfgP2U{"cfgP2U", -0.0118451, "p_{2} for upper cut for N_{TPC} vs N{sim}"}; + Configurable cfgP0L{"cfgP0L", 23.8213, "p_{0} for lower cut for N_{TPC} vs N{sim}"}; + Configurable cfgP1L{"cfgP1L", -0.402328, "p_{1} for lower cut for N_{TPC} vs N{sim}"}; + Configurable cfgP2L{"cfgP2L", 0.0128423, "p_{2} for lower cut for N_{TPC} vs N{sim}"}; + Configurable cfgSel8{"cfgSel8", true, "Sel8 trigger"}; + Configurable cfgPosZ{"cfgPosZ", true, "Position Z"}; + Configurable cfgCentFT0C{"cfgCentFT0C", false, "Use FT0C for centrality instead of FT0M"}; + Configurable cfgCentFV0A{"cfgCentFV0A", false, "Use FV0A for centrality instead of FT0M"}; + Configurable cfgNoSameBunchPileup{"cfgNoSameBunchPileup", true, "kNoSameBunchPileup"}; + Configurable cfgIsVertexITSTPC{"cfgIsVertexITSTPC", true, "kIsVertexITSTPC"}; + Configurable cfgLightIonCuts{"cfgLightIonCuts", false, "Light Ion Cuts"}; + Configurable cfgIsGoodITSLayersAll{"cfgIsGoodITSLayersAll", true, "kIsGoodITSLayersAll"}; + Configurable cfgIsGoodZvtxFT0vsPV{"cfgIsGoodZvtxFT0vsPV", true, "kIsGoodZvtxFT0vsPV"}; + Configurable cfgNtpcEventCut{"cfgNtpcEventCut", false, "N_{TPC} Event Cut "}; + } eventCuts; + struct : ConfigurableGroup { + Configurable cfgCutNSig3{"cfgCutNSig3", 3.0, "nSigma cut: 3"}; + Configurable cfgCutNSig5{"cfgCutNSig5", 5.0, "nSigma cut: 5"}; + Configurable cfgCutNSigTpcTof{"cfgCutNSigTpcTof", 3.0, "nSigma TPC cut when TPC + TOF "}; + Configurable cfgSelCutNSigTpcPi{"cfgSelCutNSigTpcPi", 3.0, "nSigma TPC cut for pion selection"}; + Configurable cfgSelCutNSigTpcKa{"cfgSelCutNSigTpcKa", 3.0, "nSigma TPC cut for kaon selection"}; + Configurable cfgSelCutNSigTpcPr{"cfgSelCutNSigTpcPr", 3.0, "nSigma TPC cut for proton selection"}; + Configurable cfgSelCutNSigTofPi{"cfgSelCutNSigTofPi", 3.0, "nSigma TOF cut for pion selection"}; + Configurable cfgSelCutNSigTofKa{"cfgSelCutNSigTofKa", 3.0, "nSigma TOF cut for kaon selection"}; + Configurable cfgSelCutNSigTofPr{"cfgSelCutNSigTofPr", 3.0, "nSigma TOF cut for proton selection"}; + Configurable cfgTofShiftPi{"cfgTofShiftPi", 0.f, "Pion TOF nSigma offset to subtract (not a time shift in ps)"}; + Configurable cfgTofShiftKa{"cfgTofShiftKa", 0.f, "Kaon TOF nSigma offset to subtract (not a time shift in ps)"}; + Configurable cfgTofShiftPr{"cfgTofShiftPr", 0.f, "Proton TOF nSigma offset to subtract (not a time shift in ps)"}; + Configurable cfgCutPiPtMin{"cfgCutPiPtMin", 0.2, "Minimum pion p_{T} cut"}; + Configurable cfgCutKaPtMin{"cfgCutKaPtMin", 0.3, "Minimum kaon p_{T} cut"}; + Configurable cfgCutPrPtMin{"cfgCutPrPtMin", 0.5, "Minimum proton p_{T} cut"}; + Configurable cfgCutPiThrsldP{"cfgCutPiThrsldP", 0.7, "Momentum threshold for requiring pion TOF PID"}; + Configurable cfgCutKaThrsldP{"cfgCutKaThrsldP", 0.55, "Momentum threshold for requiring kaon TOF PID"}; + Configurable cfgCutPrThrsldP{"cfgCutPrThrsldP", 0.8, "Momentum threshold for requiring proton TOF PID"}; + Configurable cfgRejTrk{"cfgRejTrk", true, "Rejected Tracks"}; + } pidConfig; Configurable cfgP0corr{"cfgP0corr", 0.0001, "p_{0} for corrected N_{TPC} "}; Configurable cfgP1corr{"cfgP1corr", 0.0001, "p_{1} for corrected N_{TPC} "}; Configurable cfgP2corr{"cfgP2corr", 0.0001, "p_{2} for corrected N_{TPC} "}; - Configurable cfgSel8{"cfgSel8", true, "Sel8 trigger"}; - Configurable cfgCentFT0C{"cfgCentFT0C", false, "Use FT0C for centrality instead of FT0M"}; - Configurable cfgNoSameBunchPileup{"cfgNoSameBunchPileup", true, "kNoSameBunchPileup"}; - Configurable cfgIsVertexITSTPC{"cfgIsVertexITSTPC", true, "kIsVertexITSTPC"}; - Configurable cfgLightIonCuts{"cfgLightIonCuts", false, "Light Ion Cuts"}; - Configurable cfgIsGoodITSLayersAll{"cfgIsGoodITSLayersAll", true, "kIsGoodITSLayersAll"}; - Configurable cfgIsGoodZvtxFT0vsPV{"cfgIsGoodZvtxFT0vsPV", true, "kIsGoodZvtxFT0vsPV"}; - Configurable cfgRejTrk{"cfgRejTrk", true, "Rejected Tracks"}; - Configurable cfgNtpcEventCut{"cfgNtpcEventCut", true, "N_{TPC} Event Cut for reco"}; Configurable cfgApplyEfficiencyCorrection{"cfgApplyEfficiencyCorrection", false, "Apply measured track-by-track efficiency weights"}; Configurable cfgUseCombinedPurityEfficiencyWeight{"cfgUseCombinedPurityEfficiencyWeight", false, "Use the legacy Ngen/Nselected weight instead of the efficiency-only 1/epsilon weight"}; Configurable cfgEfficiencyCCDBUrl{"cfgEfficiencyCCDBUrl", "http://ccdb-test.cern.ch:8080", "Efficiency CCDB URL"}; @@ -122,6 +142,16 @@ struct MeanPtFlucId { Configurable> ptBins{"ptBins", {0.20, 0.25, 0.30, 0.35, 0.40, 0.45, 0.50, 0.55, 0.60, 0.65, 0.70, 0.75, 0.80, 0.85, 0.90, 0.95, 1.00, 1.05, 1.10, 1.15, 1.20, 1.25, 1.30, 1.35, 1.40, 1.45, 1.50, 1.55, 1.60, 1.65, 1.70, 1.75, 1.80, 1.85, 1.90, 1.95, 2.00}, "p_{T} bins"}; Configurable> etaBins{"etaBins", {-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}, "#eta bins"}; + struct : ConfigurableGroup { + Configurable cfgGenCutPosZ{"cfgGenCutPosZ", false, "Apply a generated collision vertex-Z cut in generator-only mode"}; + Configurable cfgGenMaxPosZ{"cfgGenMaxPosZ", 7.f, "Maximum absolute generated vertex Z (cm)"}; + Configurable cfgGenCentralityCalibration{"cfgGenCentralityCalibration", false, "Collect forward multiplicity only; derive percentile lookup before analysis"}; + Configurable cfgGenForwardMultMax{"cfgGenForwardMultMax", 10000, "Largest integer forward multiplicity in calibration histogram"}; + Configurable> cfgGenCentralityLookup{"cfgGenCentralityLookup", {}, "Percentile indexed by integer generated forward multiplicity (fixed reference sample)"}; + Configurable cfgGenCentralityEstimator{"cfgGenCentralityEstimator", "", "Estimator of the supplied lookup: FT0M, FT0C or FV0A"}; + Configurable cfgGenUseEventWeight{"cfgGenUseEventWeight", false, "Use positive mcCollision.weight() consistently for generator calibration, events and matrices"}; + } genConfig; + Service pdg{}; Service ccdb{}; @@ -172,10 +202,14 @@ struct MeanPtFlucId { Gen_Charged, Gen_Pion, Gen_Kaon, - Gen_Proton + Gen_Proton, + Model_Charged, + Model_Pion, + Model_Kaon, + Model_Proton }; - static constexpr std::array Directories{ + static constexpr std::array Directories{ "QA/after/", "QA/Pion/", "QA/Kaon/", @@ -187,7 +221,11 @@ struct MeanPtFlucId { "Gen/Charged/", "Gen/Pion/", "Gen/Kaon/", - "Gen/Proton/"}; + "Gen/Proton/", + "Model/Charged/", + "Model/Pion/", + "Model/Kaon/", + "Model/Proton/"}; static constexpr auto const& dire = Directories; enum CorrectedMode { @@ -245,20 +283,20 @@ struct MeanPtFlucId { void init(InitContext const&) { - if (!std::isfinite(cfgTofShiftPi.value) || !std::isfinite(cfgTofShiftKa.value) || !std::isfinite(cfgTofShiftPr.value)) { + if (!std::isfinite(pidConfig.cfgTofShiftPi.value) || !std::isfinite(pidConfig.cfgTofShiftKa.value) || !std::isfinite(pidConfig.cfgTofShiftPr.value)) { LOGF(fatal, "TOF nSigma shifts must be finite: pion=%g, kaon=%g, proton=%g", - cfgTofShiftPi.value, cfgTofShiftKa.value, cfgTofShiftPr.value); + pidConfig.cfgTofShiftPi.value, pidConfig.cfgTofShiftKa.value, pidConfig.cfgTofShiftPr.value); } - if (nPtBins <= 0 || cfgCutPtMin >= cfgCutPtMax || cfgCutPiPtMin >= cfgCutPtMax || - cfgCutKaPtMin >= cfgCutPtMax || cfgCutPrPtMin >= cfgCutPtMax) { + if (nPtBins <= 0 || trackCuts.cfgCutPtMin >= trackCuts.cfgCutPtMax || pidConfig.cfgCutPiPtMin >= trackCuts.cfgCutPtMax || + pidConfig.cfgCutKaPtMin >= trackCuts.cfgCutPtMax || pidConfig.cfgCutPrPtMin >= trackCuts.cfgCutPtMax) { LOGF(fatal, "Invalid pT matrix binning: nPtBins=%d, max=%g, minima=(%g,%g,%g,%g)", - static_cast(nPtBins), static_cast(cfgCutPtMax), - static_cast(cfgCutPtMin), static_cast(cfgCutPiPtMin), - static_cast(cfgCutKaPtMin), static_cast(cfgCutPrPtMin)); + static_cast(nPtBins), static_cast(trackCuts.cfgCutPtMax), + static_cast(trackCuts.cfgCutPtMin), static_cast(pidConfig.cfgCutPiPtMin), + static_cast(pidConfig.cfgCutKaPtMin), static_cast(pidConfig.cfgCutPrPtMin)); } - if (cfgApplyEfficiencyCorrection) { + if (cfgApplyEfficiencyCorrection && (doprocessRun3 || doprocessMCRecoSimRun3)) { ccdb->setURL(cfgEfficiencyCCDBUrl.value); ccdb->setCaching(true); auto* efficiencyList = ccdb->getForTimeStamp(cfgEfficiencyCCDBPath.value, cfgEfficiencyTimestamp.value); @@ -282,10 +320,10 @@ struct MeanPtFlucId { auto validateEfficiencyMap = [&](const TH3* map, const char* name, double minPt) { constexpr double AxisTolerance = 1.e-6; constexpr double MaximumCentrality = 100.; - const bool coversEta = map->GetXaxis()->GetXmin() <= -cfgCutEta + AxisTolerance && - map->GetXaxis()->GetXmax() >= cfgCutEta - AxisTolerance; + const bool coversEta = map->GetXaxis()->GetXmin() <= -trackCuts.cfgCutEta + AxisTolerance && + map->GetXaxis()->GetXmax() >= trackCuts.cfgCutEta - AxisTolerance; const bool coversPt = map->GetYaxis()->GetXmin() <= minPt + AxisTolerance && - map->GetYaxis()->GetXmax() >= cfgCutPtMax - AxisTolerance; + map->GetYaxis()->GetXmax() >= trackCuts.cfgCutPtMax - AxisTolerance; const bool coversCentrality = map->GetZaxis()->GetXmin() <= AxisTolerance && map->GetZaxis()->GetXmax() >= MaximumCentrality - AxisTolerance; if (!coversEta || !coversPt || !coversCentrality) { @@ -297,20 +335,20 @@ struct MeanPtFlucId { map->GetXaxis()->GetXmin(), map->GetXaxis()->GetXmax(), map->GetYaxis()->GetXmin(), map->GetYaxis()->GetXmax(), map->GetZaxis()->GetXmin(), map->GetZaxis()->GetXmax(), - -static_cast(cfgCutEta), static_cast(cfgCutEta), - minPt, static_cast(cfgCutPtMax)); + -static_cast(trackCuts.cfgCutEta), static_cast(trackCuts.cfgCutEta), + minPt, static_cast(trackCuts.cfgCutPtMax)); } }; - validateEfficiencyMap(efficiencyWeightCh, chargedMapName, cfgCutPtMin); - validateEfficiencyMap(efficiencyWeightPi, pionMapName, cfgCutPiPtMin); - validateEfficiencyMap(efficiencyWeightKa, kaonMapName, cfgCutKaPtMin); - validateEfficiencyMap(efficiencyWeightPr, protonMapName, cfgCutPrPtMin); + validateEfficiencyMap(efficiencyWeightCh, chargedMapName, trackCuts.cfgCutPtMin); + validateEfficiencyMap(efficiencyWeightPi, pionMapName, pidConfig.cfgCutPiPtMin); + validateEfficiencyMap(efficiencyWeightKa, kaonMapName, pidConfig.cfgCutKaPtMin); + validateEfficiencyMap(efficiencyWeightPr, protonMapName, pidConfig.cfgCutPrPtMin); } - const AxisSpec axisPtCh{nPtBins, cfgCutPtMin, cfgCutPtMax, "p_{T} (GeV/c)"}; - const AxisSpec axisPtPi{nPtBins, cfgCutPiPtMin, cfgCutPtMax, "p_{T} (GeV/c)"}; - const AxisSpec axisPtKa{nPtBins, cfgCutKaPtMin, cfgCutPtMax, "p_{T} (GeV/c)"}; - const AxisSpec axisPtPr{nPtBins, cfgCutPrPtMin, cfgCutPtMax, "p_{T} (GeV/c)"}; + const AxisSpec axisPtCh{nPtBins, trackCuts.cfgCutPtMin, trackCuts.cfgCutPtMax, "p_{T} (GeV/c)"}; + const AxisSpec axisPtPi{nPtBins, pidConfig.cfgCutPiPtMin, trackCuts.cfgCutPtMax, "p_{T} (GeV/c)"}; + const AxisSpec axisPtKa{nPtBins, pidConfig.cfgCutKaPtMin, trackCuts.cfgCutPtMax, "p_{T} (GeV/c)"}; + const AxisSpec axisPtPr{nPtBins, pidConfig.cfgCutPrPtMin, trackCuts.cfgCutPtMax, "p_{T} (GeV/c)"}; const AxisSpec axisCol{3, 1, 4, ""}; const AxisSpec axisTrack{5, 1, 6, ""}; @@ -325,8 +363,10 @@ struct MeanPtFlucId { const AxisSpec axisMultCorr{multCorrBins, "N_{Corr} "}; const AxisSpec axisMultMC{multMCBins, "N_{Gen} "}; const AxisSpec axisMultFT0M{multFT0MBins, "N_{FT0M}"}; - const AxisSpec axisCentFT0M{101, 0, 101, "FT0M (%)"}; - const AxisSpec axisCent{centBins, "FT0M (%)"}; + // Match the estimator precedence used by getConfiguredCentrality(). + const std::string centralityLabel = eventCuts.cfgCentFT0C ? "FT0C (%)" : (eventCuts.cfgCentFV0A ? "FV0A (%)" : "FT0M (%)"); + const AxisSpec axisCentFT0M{101, 0, 101, centralityLabel}; + const AxisSpec axisCent{centBins, centralityLabel}; const AxisSpec axisVtxZ{80, -20., 20., "V_{Z} (cm)"}; const AxisSpec axisDCAz{dcaZBins, "DCA_{Z} (cm)"}; const AxisSpec axisDCAxy{dcaXYBins, "DCA_{XY} (cm)"}; @@ -353,7 +393,7 @@ struct MeanPtFlucId { hist.add("QA/before/h_VtxZ", "V_{Z}", kTH1D, {axisVtxZ}); hist.add("QA/before/h_NTPC", "N_{TPC}", kTH1D, {axisMultTPC}); hist.add("QA/before/h_NFT0M", "FT0M Multiplicity", kTH1D, {axisMultFT0M}); - hist.add("QA/before/h_CentM", "FT0M (%)", kTH1D, {axisCentFT0M}); + hist.add("QA/before/h_CentM", centralityLabel.c_str(), kTH1D, {axisCentFT0M}); hist.add("QA/before/h2_TPCSignal", "TPC Signal", tpcSignalHist); hist.add("QA/before/h2_TOFSignal", "TOF Signal", tofSignalHist); @@ -374,10 +414,10 @@ struct MeanPtFlucId { hist.add("QA/after/h_TPCChi2perCluster", "TPC #Chi^{2}/Cluster", kTH1D, {axisChi2}); hist.add("QA/after/h_ITSChi2perCluster", "ITS #Chi^{2}/Cluster", kTH1D, {axisChi2}); hist.add("QA/after/h_crossedTPC", "Crossed TPC", kTH1D, {axisCrossedTPC}); - hist.add("QA/after/h2_NTPC_CentM", "N_{TPC} vs FT0M(%)", kTH2D, {{axisCentFT0M}, {axisMultTPC}}); + hist.add("QA/after/h2_NTPC_CentM", ("N_{TPC} vs " + centralityLabel).c_str(), kTH2D, {{axisCentFT0M}, {axisMultTPC}}); hist.add("QA/after/h2_NTPC_NFT0M", "N_{TPC} vs N_{FT0M}", kTH2D, {{axisMultFT0M}, {axisMultTPC}}); hist.add("QA/after/p_NTPC_NFT0M", "N_{TPC} vs N_{FT0M} (Profile)", kTProfile, {axisMultFT0M}); - hist.add("QA/after/p_NTPC_CentM", "N_{TPC} vs FT0M(%) (Profile)", kTProfile, {axisCentFT0M}); + hist.add("QA/after/p_NTPC_CentM", ("N_{TPC} vs " + centralityLabel + " (Profile)").c_str(), kTProfile, {axisCentFT0M}); hist.add("QA/after/h_DCAxy_primary", "DCA_{XY} (Primary)", kTH1D, {axisDCAxy}); hist.add("QA/after/h_DCAz_primary", "DCA_{Z} (Primary)", kTH1D, {axisDCAz}); hist.add("QA/after/h_DCAxy_secondary", "DCA_{XY} (Secondary)", kTH1D, {axisDCAxy}); @@ -387,7 +427,7 @@ struct MeanPtFlucId { hist.add("QA/after/h2_NTPC_NSim", "Truth vs Reco Multiplicty TPC", kTH2D, {{axisMultTPC}, {axisMultMC}}); hist.add("QA/after/p_NTPC_NSim", "Truth vs Reco Multiplicty TPC", kTProfile, {axisMultTPC}); hist.add("QA/after/h2_Ncorr_NSim", "Truth vs Corrected Multiplicty TPC", kTH2D, {{axisMultCorr}, {axisMultMC}}); - hist.add("QA/after/h2_Ncorr_CentM", "N_{corr} vs FT0M(%)", kTH2D, {{axisCentFT0M}, {axisMultCorr}}); + hist.add("QA/after/h2_Ncorr_CentM", ("N_{corr} vs " + centralityLabel).c_str(), kTH2D, {{axisCentFT0M}, {axisMultCorr}}); }; addEventQAHistograms(); @@ -427,7 +467,7 @@ struct MeanPtFlucId { hist.add("QA/Charged/h2_Pt_Eta", "p_{T} vs #eta ", kTH2D, {{axisEta}, {axisPtCh}}); hist.add("QA/Charged/h3_Pt_Eta_centFT0M", "p_{T} vs #eta in centrality classes", kTH3D, {{axisEta}, {axisPtCh}, {axisCent}}); hist.add("QA/Charged/h3_Pt_EtaMC_centFT0M", "p_{T} vs #eta in centrality classes", kTH3D, {{axisEta}, {axisPtCh}, {axisCent}}); - hist.add("QA/Charged/h3_Pt_EtaTruth_primary_centFT0M", "Selected truth primary tracks in reconstructed coordinates;#eta;p_{T};FT0M (%)", kTH3D, {{axisEta}, {axisPtCh}, {axisCent}}); + hist.add("QA/Charged/h3_Pt_EtaTruth_primary_centFT0M", ("Selected truth primary tracks in reconstructed coordinates;#eta;p_{T};" + centralityLabel).c_str(), kTH3D, {{axisEta}, {axisPtCh}, {axisCent}}); hist.add("QA/Charged/h2_Pt_EtaMC", "p_{T} vs #eta (MC)", kTH2D, {{axisEta}, {axisPtCh}}); }; addTrackQAHistograms(); @@ -438,7 +478,7 @@ struct MeanPtFlucId { hist.add("QA/Pion/h2_Pt_Eta", "p_{T} vs #eta ", kTH2D, {{axisEta}, {axisPtPi}}); hist.add("QA/Pion/h3_Pt_Eta_centFT0M", "p_{T} vs #eta in centrality classes", kTH3D, {{axisEta}, {axisPtPi}, {axisCent}}); hist.add("QA/Pion/h3_Pt_EtaMC_centFT0M", "p_{T} vs #eta in centrality classes", kTH3D, {{axisEta}, {axisPtPi}, {axisCent}}); - hist.add("QA/Pion/h3_Pt_EtaTruth_primary_centFT0M", "Selected truth primary tracks in reconstructed coordinates;#eta;p_{T};FT0M (%)", kTH3D, {{axisEta}, {axisPtPi}, {axisCent}}); + hist.add("QA/Pion/h3_Pt_EtaTruth_primary_centFT0M", ("Selected truth primary tracks in reconstructed coordinates;#eta;p_{T};" + centralityLabel).c_str(), kTH3D, {{axisEta}, {axisPtPi}, {axisCent}}); hist.add("QA/Pion/h2_Pt_EtaMC", "p_{T} vs #eta (MC)", kTH2D, {{axisEta}, {axisPtPi}}); hist.add("QA/Pion/h_PtPos", "p_{T} (positive) ", kTH1D, {axisPtPi}); hist.add("QA/Pion/h_PtNeg", "p_{T} (negative) ", kTH1D, {axisPtPi}); @@ -455,7 +495,7 @@ struct MeanPtFlucId { hist.add("QA/Kaon/h2_Pt_Eta", "p_{T} vs #eta ", kTH2D, {{axisEta}, {axisPtKa}}); hist.add("QA/Kaon/h3_Pt_Eta_centFT0M", "p_{T} vs #eta in centrality classes", kTH3D, {{axisEta}, {axisPtKa}, {axisCent}}); hist.add("QA/Kaon/h3_Pt_EtaMC_centFT0M", "p_{T} vs #eta in centrality classes", kTH3D, {{axisEta}, {axisPtKa}, {axisCent}}); - hist.add("QA/Kaon/h3_Pt_EtaTruth_primary_centFT0M", "Selected truth primary tracks in reconstructed coordinates;#eta;p_{T};FT0M (%)", kTH3D, {{axisEta}, {axisPtKa}, {axisCent}}); + hist.add("QA/Kaon/h3_Pt_EtaTruth_primary_centFT0M", ("Selected truth primary tracks in reconstructed coordinates;#eta;p_{T};" + centralityLabel).c_str(), kTH3D, {{axisEta}, {axisPtKa}, {axisCent}}); hist.add("QA/Kaon/h2_Pt_EtaMC", "p_{T} vs #eta (MC)", kTH2D, {{axisEta}, {axisPtKa}}); hist.add("QA/Kaon/h_PtPos", "p_{T} (positive) ", kTH1D, {axisPtKa}); hist.add("QA/Kaon/h_PtNeg", "p_{T} (negative) ", kTH1D, {axisPtKa}); @@ -472,7 +512,7 @@ struct MeanPtFlucId { hist.add("QA/Proton/h2_Pt_Eta", "p_{T} vs #eta ", kTH2D, {{axisEta}, {axisPtPr}}); hist.add("QA/Proton/h3_Pt_Eta_centFT0M", "p_{T} vs #eta in centrality classes", kTH3D, {{axisEta}, {axisPtPr}, {axisCent}}); hist.add("QA/Proton/h3_Pt_EtaMC_centFT0M", "p_{T} vs #eta in centrality classes", kTH3D, {{axisEta}, {axisPtPr}, {axisCent}}); - hist.add("QA/Proton/h3_Pt_EtaTruth_primary_centFT0M", "Selected truth primary tracks in reconstructed coordinates;#eta;p_{T};FT0M (%)", kTH3D, {{axisEta}, {axisPtPr}, {axisCent}}); + hist.add("QA/Proton/h3_Pt_EtaTruth_primary_centFT0M", ("Selected truth primary tracks in reconstructed coordinates;#eta;p_{T};" + centralityLabel).c_str(), kTH3D, {{axisEta}, {axisPtPr}, {axisCent}}); hist.add("QA/Proton/h2_Pt_EtaMC", "p_{T} vs #eta (MC)", kTH2D, {{axisEta}, {axisPtPr}}); hist.add("QA/Proton/h_PtPos", "p_{T} (positive) ", kTH1D, {axisPtPr}); hist.add("QA/Proton/h_PtNeg", "p_{T} (negative) ", kTH1D, {axisPtPr}); @@ -567,6 +607,53 @@ struct MeanPtFlucId { }; addGeneratedHistograms(); + if (doprocessMCGenOnly) { + if (genConfig.cfgGenCutPosZ && (!std::isfinite(genConfig.cfgGenMaxPosZ.value) || genConfig.cfgGenMaxPosZ <= 0.f)) { + LOGF(fatal, "cfgGenMaxPosZ must be finite and positive"); + } + const std::string estimator = eventCuts.cfgCentFT0C ? "FT0C" : (eventCuts.cfgCentFV0A ? "FV0A" : "FT0M"); + constexpr size_t MinimumCentralityLookupEntries = 2; + constexpr double MinimumCentralityPercentile = 0.; + constexpr double MaximumCentralityPercentile = 100.; + if (genConfig.cfgGenForwardMultMax < 1) { + LOGF(fatal, "cfgGenForwardMultMax must be positive"); + } + if (!genConfig.cfgGenCentralityCalibration) { + const auto& lookup = genConfig.cfgGenCentralityLookup.value; + if (genConfig.cfgGenCentralityEstimator.value != estimator || lookup.size() < MinimumCentralityLookupEntries) { + LOGF(fatal, "Supply a centrality lookup for %s, or enable cfgGenCentralityCalibration for the first pass", estimator.c_str()); + } + for (size_t i = 0; i < lookup.size(); ++i) { + if (!std::isfinite(lookup[i]) || lookup[i] < MinimumCentralityPercentile || lookup[i] > MaximumCentralityPercentile || (i != 0 && lookup[i] > lookup[i - 1])) { + LOGF(fatal, "Generator centrality lookup must be finite, within [0,100] and non-increasing"); + } + } + } + const AxisSpec axisModelCent{centBins, "Generated " + estimator + " percentile (%)"}; + const AxisSpec axisForward{genConfig.cfgGenForwardMultMax.value + 1, -0.5, genConfig.cfgGenForwardMultMax.value + 0.5, "Generated " + estimator + " multiplicity"}; + hist.add("Model/h_ForwardMult", "Reference forward multiplicity after generator event cuts", kTH1D, {axisForward}); + hist.add("Model/h_Centrality", "Events used for matrices (including zero selected particles)", kTH1D, {axisModelCent}); + hist.add("Model/h_CentralityStatus", "Centrality lookup;Status;Events", kTH1D, {{2, 0.5, 2.5}}); + hist.get(HIST("Model/h_CentralityStatus"))->GetXaxis()->SetBinLabel(1, "Calibrated"); + hist.get(HIST("Model/h_CentralityStatus"))->GetXaxis()->SetBinLabel(2, "Outside lookup"); + const AxisSpec axisModelMult{multMCBins, "N_{ch}^{gen} (physical primaries in |#eta| acceptance, all p_{T})"}; + hist.add("Model/h_EventCount", "Generated events;Selection;Events", kTH1D, {{2, 0.5, 2.5}}); + hist.get(HIST("Model/h_EventCount"))->GetXaxis()->SetBinLabel(1, "All generated"); + hist.get(HIST("Model/h_EventCount"))->GetXaxis()->SetBinLabel(2, "Accepted"); + hist.add("Model/h_VtxZ", "Generated vertex Z", kTH1D, {axisVtxZ}); + hist.add("Model/h_Nch", "Accepted events (including zero selected particles)", kTH1D, {axisModelMult}); + auto addModelSpecies = [&](const std::string& directory, const AxisSpec& axisPtId) { + hist.add(directory + "h_Pt", "Generated physical primaries", kTH1D, {axisPtId}); + hist.add(directory + "h_Eta", "Generated physical primaries", kTH1D, {axisEta}); + hist.add(directory + "hN1Matrix", "Generated single-particle density", kTH2D, {axisModelCent, axisPtId}); + hist.add(directory + "hN2Matrix", "Generated ordered-pair density", kTH3D, {axisModelCent, axisPtId, axisPtId}); + }; + addModelSpecies("Model/Charged/", axisPtCh); + addModelSpecies("Model/Pion/", axisPtPi); + addModelSpecies("Model/Kaon/", axisPtKa); + addModelSpecies("Model/Proton/", axisPtPr); + } + hist.add("QA/h_collisions_info", "Collisions info", kTH1D, {axisCol}); hist.add("Gen/h_collisions_info", "Collisions info", kTH1D, {axisCol}); hist.add("Gen/h_collision_recgen", "Number of Collisions ", kTH1D, {axisCol}); @@ -586,7 +673,7 @@ struct MeanPtFlucId { hist.get(HIST("QA/after/h_counts_evSelCuts"))->GetXaxis()->SetBinLabel(SelCollisionLabels::kSelColSel8, "kSelColSel8"); hist.get(HIST("QA/after/h_counts_evSelCuts"))->GetXaxis()->SetBinLabel(SelCollisionLabels::kSelColNoSameBunchPileup, "kSelColNoSameBunchPileup"); hist.get(HIST("QA/after/h_counts_evSelCuts"))->GetXaxis()->SetBinLabel(SelCollisionLabels::kSelColIsVertexITSTPC, "kSelColIsVertexITSTPC"); - if (cfgLightIonCuts) { + if (eventCuts.cfgLightIonCuts) { hist.get(HIST("QA/after/h_counts_evSelCuts"))->GetXaxis()->SetBinLabel(SelCollisionLabels::kSelColIsGoodITSLayersAll, "kSelColIsGoodITSLayersAll"); hist.get(HIST("QA/after/h_counts_evSelCuts"))->GetXaxis()->SetBinLabel(SelCollisionLabels::kSelColIsGoodZvtxFT0vsPV, "kSelColIsGoodZvtxFT0vsPV"); } @@ -596,54 +683,62 @@ struct MeanPtFlucId { int nTPC = 0, nFT0M = 0; int nSim = 0; + template + float getConfiguredCentrality(T const& col) const + { + if (eventCuts.cfgCentFT0C.value) { + return col.centFT0C(); + } + if (eventCuts.cfgCentFV0A.value) { + return col.centFV0A(); + } + return col.centFT0M(); + } + // Event selection cuts: template bool selRun3Col(T const& col) { - if (cfgCentFT0C) { - centFT0M = col.centFT0C(); - } else { - centFT0M = col.centFT0M(); - } + centFT0M = getConfiguredCentrality(col); nTPC = col.multNTracksHasTPC(); nFT0M = col.multFT0M(); hist.fill(HIST("QA/after/h_counts_evSelCuts"), kBeforeSelCol); - if (cfgPosZ) { - if (std::abs(col.posZ()) > cfgCutPosZ) { + if (eventCuts.cfgPosZ) { + if (std::abs(col.posZ()) > eventCuts.cfgCutPosZ) { return false; } hist.fill(HIST("QA/after/h_counts_evSelCuts"), kSelColPosZ); } - if (cfgSel8) { + if (eventCuts.cfgSel8) { if (!col.sel8()) { return false; } hist.fill(HIST("QA/after/h_counts_evSelCuts"), kSelColSel8); } - if (cfgNoSameBunchPileup) { + if (eventCuts.cfgNoSameBunchPileup) { if (!col.selection_bit(o2::aod::evsel::kNoSameBunchPileup)) { return false; } hist.fill(HIST("QA/after/h_counts_evSelCuts"), kSelColNoSameBunchPileup); } - if (cfgIsVertexITSTPC) { + if (eventCuts.cfgIsVertexITSTPC) { if (!col.selection_bit(o2::aod::evsel::kIsVertexITSTPC)) { return false; } hist.fill(HIST("QA/after/h_counts_evSelCuts"), kSelColIsVertexITSTPC); } - if (cfgLightIonCuts) { - if (cfgIsGoodITSLayersAll) { + if (eventCuts.cfgLightIonCuts) { + if (eventCuts.cfgIsGoodITSLayersAll) { if (!col.selection_bit(o2::aod::evsel::kIsGoodITSLayersAll)) { return false; } hist.fill(HIST("QA/after/h_counts_evSelCuts"), kSelColIsGoodITSLayersAll); } - if (cfgIsGoodZvtxFT0vsPV) { + if (eventCuts.cfgIsGoodZvtxFT0vsPV) { if (!col.selection_bit(o2::aod::evsel::kIsGoodZvtxFT0vsPV)) { return false; } @@ -658,15 +753,19 @@ struct MeanPtFlucId { template bool selTrack(T const& track) { - if (!track.isGlobalTrack()) { + if (trackCuts.cfgGlobalTrack && !track.isGlobalTrack()) { return false; } - if (track.pt() < cfgCutPtMin) { + if (trackCuts.cfgGlobalTrackWoDCA && !track.isGlobalTrackWoDCA()) { return false; } - if (track.pt() >= cfgCutPtMax) { + if (track.pt() < trackCuts.cfgCutPtMin) { + return false; + } + + if (track.pt() >= trackCuts.cfgCutPtMax) { return false; } @@ -674,16 +773,41 @@ struct MeanPtFlucId { return false; } - const float maxDca = cfgDcaPtDepOffset + cfgDcaPtDepScale / track.pt(); - // if (std::abs(track.dcaXY()) >= maxDca) { - // return false; - // } + const float maxDca = trackCuts.cfgDcaPtDepOffset + trackCuts.cfgDcaPtDepScale / track.pt(); + + if (trackCuts.cfgDcaXYPtDep) { + if (std::abs(track.dcaXY()) >= maxDca) { + return false; + } + } else if (trackCuts.cfgDcaXY) { + if (std::abs(track.dcaXY()) >= trackCuts.cfgCutDCAxy) { + return false; + } + } + + if (trackCuts.cfgDcaZPtDep) { + if (std::abs(track.dcaZ()) >= maxDca) { + return false; + } + } else if (trackCuts.cfgDcaZ) { + if (std::abs(track.dcaZ()) >= trackCuts.cfgCutDCAz) { + return false; + } + } + + if (std::abs(track.eta()) >= trackCuts.cfgCutEta) { + return false; + } + + if (trackCuts.cfgTPCChi2NCl && track.tpcChi2NCl() > trackCuts.cfgCutTPCChi2NCl) { + return false; + } - if (std::abs(track.dcaZ()) >= maxDca) { + if (trackCuts.cfgITSChi2NCl && track.itsChi2NCl() > trackCuts.cfgCutITSChi2NCl) { return false; } - if (std::abs(track.eta()) >= cfgCutEta) { + if (trackCuts.cfgNTPCCrossedRows && track.tpcNClsCrossedRows() < trackCuts.cfgCutNCrossedRowsTPC) { return false; } @@ -694,11 +818,11 @@ struct MeanPtFlucId { template bool rejectTracks(T const& track) { - return ((track.tpcNSigmaEl()) > -cfgCutNSig3 && - (track.tpcNSigmaEl()) < cfgCutNSig5) && - (std::fabs(track.tpcNSigmaPi()) > cfgCutNSig3 && - std::fabs(track.tpcNSigmaKa()) > cfgCutNSig3 && - std::fabs(track.tpcNSigmaPr()) > cfgCutNSig3); + return ((track.tpcNSigmaEl()) > -pidConfig.cfgCutNSig3 && + (track.tpcNSigmaEl()) < pidConfig.cfgCutNSig5) && + (std::fabs(track.tpcNSigmaPi()) > pidConfig.cfgCutNSig3 && + std::fabs(track.tpcNSigmaKa()) > pidConfig.cfgCutNSig3 && + std::fabs(track.tpcNSigmaPr()) > pidConfig.cfgCutNSig3); } // Share the corrected hypotheses between PID selection and selected-particle QA. @@ -708,9 +832,9 @@ struct MeanPtFlucId { { std::array values = {-999.f, track.tofNSigmaPi(), track.tofNSigmaKa(), track.tofNSigmaPr()}; if (track.hasTOF()) { - values[static_cast(PIDType::kPions)] -= cfgTofShiftPi.value; - values[static_cast(PIDType::kKaons)] -= cfgTofShiftKa.value; - values[static_cast(PIDType::kProtons)] -= cfgTofShiftPr.value; + values[static_cast(PIDType::kPions)] -= pidConfig.cfgTofShiftPi.value; + values[static_cast(PIDType::kKaons)] -= pidConfig.cfgTofShiftKa.value; + values[static_cast(PIDType::kProtons)] -= pidConfig.cfgTofShiftPr.value; } return values; } @@ -730,16 +854,16 @@ struct MeanPtFlucId { switch (species) { case PIDType::kPions: - tpcSelectionCut = cfgSelCutNSigTpcPi; - tofSelectionCut = cfgSelCutNSigTofPi; + tpcSelectionCut = pidConfig.cfgSelCutNSigTpcPi; + tofSelectionCut = pidConfig.cfgSelCutNSigTofPi; break; case PIDType::kKaons: - tpcSelectionCut = cfgSelCutNSigTpcKa; - tofSelectionCut = cfgSelCutNSigTofKa; + tpcSelectionCut = pidConfig.cfgSelCutNSigTpcKa; + tofSelectionCut = pidConfig.cfgSelCutNSigTofKa; break; case PIDType::kProtons: - tpcSelectionCut = cfgSelCutNSigTpcPr; - tofSelectionCut = cfgSelCutNSigTofPr; + tpcSelectionCut = pidConfig.cfgSelCutNSigTpcPr; + tofSelectionCut = pidConfig.cfgSelCutNSigTofPr; break; default: return false; @@ -750,7 +874,7 @@ struct MeanPtFlucId { const bool passesTOF = selectedTofNSigma < tofSelectionCut && std::abs(tofNSigma[rejSpecies1]) > selectedTofNSigma && std::abs(tofNSigma[rejSpecies2]) > selectedTofNSigma; - const bool passesTPC = std::abs(tpcNSigma[selSpecies]) < cfgCutNSigTpcTof; + const bool passesTPC = std::abs(tpcNSigma[selSpecies]) < pidConfig.cfgCutNSigTpcTof; return passesTPC && passesTOF; } @@ -780,7 +904,7 @@ struct MeanPtFlucId { hist.fill(HIST("QA/before/h_VtxZ"), col.posZ()); hist.fill(HIST("QA/before/h_Counts"), 2); hist.fill(HIST("QA/before/h_NTPC"), col.multNTracksHasTPC()); - hist.fill(HIST("QA/before/h_CentM"), col.centFT0M()); + hist.fill(HIST("QA/before/h_CentM"), getConfiguredCentrality(col)); hist.fill(HIST("QA/before/h_NFT0M"), col.multFT0M()); } @@ -943,10 +1067,10 @@ struct MeanPtFlucId { } template - void fillMatrices(const std::vector& tracks, double ptMin) + void fillMatrices(const std::vector& tracks, double ptMin, double eventCoordinate, double eventWeight = 1.) { const int nBins = nPtBins; - const double binWidth = (cfgCutPtMax - ptMin) / nBins; + const double binWidth = (trackCuts.cfgCutPtMax - ptMin) / nBins; matrixNInBin.assign(nBins, 0); matrixOccupiedBins.clear(); matrixOccupiedBins.reserve(std::min(static_cast(nBins), tracks.size())); @@ -964,7 +1088,7 @@ struct MeanPtFlucId { for (const int& ptBin : matrixOccupiedBins) { const double ptCenter = ptMin + (ptBin + 0.5) * binWidth; - hist.fill(HIST(dire[Mode]) + HIST("hN1Matrix"), centFT0M, ptCenter, static_cast(matrixNInBin[ptBin])); + hist.fill(HIST(dire[Mode]) + HIST("hN1Matrix"), eventCoordinate, ptCenter, eventWeight * static_cast(matrixNInBin[ptBin])); } for (size_t i = 0; i < matrixOccupiedBins.size(); ++i) { @@ -972,19 +1096,25 @@ struct MeanPtFlucId { const double ptCenter1 = ptMin + (ptBin1 + 0.5) * binWidth; const double nPairDiagonal = static_cast(matrixNInBin[ptBin1]) * (matrixNInBin[ptBin1] - 1); if (nPairDiagonal > 0.0) { - hist.fill(HIST(dire[Mode]) + HIST("hN2Matrix"), centFT0M, ptCenter1, ptCenter1, nPairDiagonal); + hist.fill(HIST(dire[Mode]) + HIST("hN2Matrix"), eventCoordinate, ptCenter1, ptCenter1, eventWeight * nPairDiagonal); } for (size_t j = i + 1; j < matrixOccupiedBins.size(); ++j) { const int ptBin2 = matrixOccupiedBins[j]; const double ptCenter2 = ptMin + (ptBin2 + 0.5) * binWidth; const double nPair = static_cast(matrixNInBin[ptBin1]) * matrixNInBin[ptBin2]; - hist.fill(HIST(dire[Mode]) + HIST("hN2Matrix"), centFT0M, ptCenter1, ptCenter2, nPair); - hist.fill(HIST(dire[Mode]) + HIST("hN2Matrix"), centFT0M, ptCenter2, ptCenter1, nPair); + hist.fill(HIST(dire[Mode]) + HIST("hN2Matrix"), eventCoordinate, ptCenter1, ptCenter2, eventWeight * nPair); + hist.fill(HIST(dire[Mode]) + HIST("hN2Matrix"), eventCoordinate, ptCenter2, ptCenter1, eventWeight * nPair); } } } + template + void fillMatrices(const std::vector& tracks, double ptMin) + { + fillMatrices(tracks, ptMin, centFT0M); + } + bool getEfficiencyWeight(PIDType species, double pt, double eta, double& weight) const { TH3* weightMap = nullptr; @@ -1023,7 +1153,7 @@ struct MeanPtFlucId { void fillCorrectedMatrices(const std::vector& tracks, double ptMin) { const int nBins = nPtBins; - const double binWidth = (cfgCutPtMax - ptMin) / nBins; + const double binWidth = (trackCuts.cfgCutPtMax - ptMin) / nBins; matrixWInBin.assign(nBins, 0.); matrixW2InBin.assign(nBins, 0.); matrixOccupiedBins.clear(); @@ -1189,14 +1319,14 @@ struct MeanPtFlucId { fillBeforePIDQAHistos(track); - if (cfgRejTrk && rejectTracks(track)) { + if (pidConfig.cfgRejTrk && rejectTracks(track)) { continue; } - auto selIDPion = identifyParticle(track, PIDType::kPions, PIDType::kKaons, PIDType::kProtons, p, cfgCutPiThrsldP, tofNSigma); - auto selIDKaon = identifyParticle(track, PIDType::kKaons, PIDType::kPions, PIDType::kProtons, p, cfgCutKaThrsldP, tofNSigma); - auto selIDProton = identifyParticle(track, PIDType::kProtons, PIDType::kPions, PIDType::kKaons, p, cfgCutPrThrsldP, tofNSigma); + auto selIDPion = identifyParticle(track, PIDType::kPions, PIDType::kKaons, PIDType::kProtons, p, pidConfig.cfgCutPiThrsldP, tofNSigma); + auto selIDKaon = identifyParticle(track, PIDType::kKaons, PIDType::kPions, PIDType::kProtons, p, pidConfig.cfgCutKaThrsldP, tofNSigma); + auto selIDProton = identifyParticle(track, PIDType::kProtons, PIDType::kPions, PIDType::kKaons, p, pidConfig.cfgCutPrThrsldP, tofNSigma); - if (selIDPion && track.pt() >= cfgCutPiPtMin) { + if (selIDPion && track.pt() >= pidConfig.cfgCutPiPtMin) { if (cfgApplyEfficiencyCorrection) { double weight = 0.; if (!getEfficiencyWeight(PIDType::kPions, pt, eta, weight)) { @@ -1206,7 +1336,7 @@ struct MeanPtFlucId { } fillIdParticleQAHistos(track, ptMC, etaMC, rapPi, nSigmaTPCPi, nSigmaTOFPi, pionTracks); } - if (selIDKaon && track.pt() >= cfgCutKaPtMin) { + if (selIDKaon && track.pt() >= pidConfig.cfgCutKaPtMin) { if (cfgApplyEfficiencyCorrection) { double weight = 0.; if (!getEfficiencyWeight(PIDType::kKaons, pt, eta, weight)) { @@ -1216,7 +1346,7 @@ struct MeanPtFlucId { } fillIdParticleQAHistos(track, ptMC, etaMC, rapKa, nSigmaTPCKa, nSigmaTOFKa, kaonTracks); } - if (selIDProton && track.pt() >= cfgCutPrPtMin) { + if (selIDProton && track.pt() >= pidConfig.cfgCutPrPtMin) { if (cfgApplyEfficiencyCorrection) { double weight = 0.; if (!getEfficiencyWeight(PIDType::kProtons, pt, eta, weight)) { @@ -1230,7 +1360,7 @@ struct MeanPtFlucId { if constexpr (RecoFlag) { auto mc = track.template mcParticle_as(); int pid = mc.pdgCode(); - if (selIDPion && track.pt() >= cfgCutPiPtMin) { + if (selIDPion && track.pt() >= pidConfig.cfgCutPiPtMin) { if (std::abs(pid) == kPiPlus) { hist.fill(HIST("QA/Pion/h_PtTruth"), track.pt()); fillPtMCHist(false, pt, eta, rapPi, phi, pid, kPiPlus, kPiMinus); @@ -1250,7 +1380,7 @@ struct MeanPtFlucId { } } } - if (selIDKaon && track.pt() >= cfgCutKaPtMin) { + if (selIDKaon && track.pt() >= pidConfig.cfgCutKaPtMin) { if (std::abs(pid) == kKPlus) { hist.fill(HIST("QA/Kaon/h_PtTruth"), track.pt()); fillPtMCHist(false, pt, eta, rapKa, phi, pid, kKPlus, kKMinus); @@ -1270,7 +1400,7 @@ struct MeanPtFlucId { } } } - if (selIDProton && track.pt() >= cfgCutPrPtMin) { + if (selIDProton && track.pt() >= pidConfig.cfgCutPrPtMin) { if (std::abs(pid) == kProton) { hist.fill(HIST("QA/Proton/h_PtTruth"), track.pt()); fillPtMCHist(false, pt, eta, rapPr, phi, pid, kProton, kProtonBar); @@ -1292,22 +1422,22 @@ struct MeanPtFlucId { } } } - fillMatrices(chargedTracks, cfgCutPtMin); - fillMatrices(pionTracks, cfgCutPiPtMin); - fillMatrices(kaonTracks, cfgCutKaPtMin); - fillMatrices(protonTracks, cfgCutPrPtMin); + fillMatrices(chargedTracks, trackCuts.cfgCutPtMin); + fillMatrices(pionTracks, pidConfig.cfgCutPiPtMin); + fillMatrices(kaonTracks, pidConfig.cfgCutKaPtMin); + fillMatrices(protonTracks, pidConfig.cfgCutPrPtMin); if (cfgApplyEfficiencyCorrection) { - fillCorrectedMatrices(chargedTracksCorrected, cfgCutPtMin); - fillCorrectedMatrices(pionTracksCorrected, cfgCutPiPtMin); - fillCorrectedMatrices(kaonTracksCorrected, cfgCutKaPtMin); - fillCorrectedMatrices(protonTracksCorrected, cfgCutPrPtMin); + fillCorrectedMatrices(chargedTracksCorrected, trackCuts.cfgCutPtMin); + fillCorrectedMatrices(pionTracksCorrected, pidConfig.cfgCutPiPtMin); + fillCorrectedMatrices(kaonTracksCorrected, pidConfig.cfgCutKaPtMin); + fillCorrectedMatrices(protonTracksCorrected, pidConfig.cfgCutPrPtMin); if constexpr (RecoFlag) { if (!cfgUseCombinedPurityEfficiencyWeight) { - fillCorrectedMatrices(chargedTracksTruthMatchedCorrected, cfgCutPtMin); - fillCorrectedMatrices(pionTracksTruthMatchedCorrected, cfgCutPiPtMin); - fillCorrectedMatrices(kaonTracksTruthMatchedCorrected, cfgCutKaPtMin); - fillCorrectedMatrices(protonTracksTruthMatchedCorrected, cfgCutPrPtMin); + fillCorrectedMatrices(chargedTracksTruthMatchedCorrected, trackCuts.cfgCutPtMin); + fillCorrectedMatrices(pionTracksTruthMatchedCorrected, pidConfig.cfgCutPiPtMin); + fillCorrectedMatrices(kaonTracksTruthMatchedCorrected, pidConfig.cfgCutKaPtMin); + fillCorrectedMatrices(protonTracksTruthMatchedCorrected, pidConfig.cfgCutPrPtMin); } } } @@ -1318,7 +1448,7 @@ struct MeanPtFlucId { { int multiplicity = 0; for (auto const& mcPart : mcParticles) { - if (!mcPart.isPhysicalPrimary() || std::abs(mcPart.eta()) >= cfgCutEta) { + if (!mcPart.isPhysicalPrimary() || std::abs(mcPart.eta()) >= trackCuts.cfgCutEta) { continue; } auto* particle = pdg->GetParticle(mcPart.pdgCode()); @@ -1331,11 +1461,11 @@ struct MeanPtFlucId { bool passesNtpcEventCut() { - if (!cfgNtpcEventCut) { + if (!eventCuts.cfgNtpcEventCut) { return true; } - const float lower = cfgP2L * nSim * nSim + cfgP1L * nSim - cfgP0L; - const float upper = cfgP2U * nSim * nSim + cfgP1U * nSim + cfgP0U; + const float lower = eventCuts.cfgP2L * nSim * nSim + eventCuts.cfgP1L * nSim - eventCuts.cfgP0L; + const float upper = eventCuts.cfgP2U * nSim * nSim + eventCuts.cfgP1U * nSim + eventCuts.cfgP0U; return nTPC >= lower && nTPC <= upper; } @@ -1368,11 +1498,11 @@ struct MeanPtFlucId { eta = mcPart.eta(); phi = mcPart.phi(); - if (std::abs(eta) < cfgCutEta) { + if (std::abs(eta) < trackCuts.cfgCutEta) { nSim++; } - if (pt >= cfgCutPtMin && pt < cfgCutPtMax && std::abs(eta) < cfgCutEta) { + if (pt >= trackCuts.cfgCutPtMin && pt < trackCuts.cfgCutPtMax && std::abs(eta) < trackCuts.cfgCutEta) { chargedTracksGen.push_back({pt}); hist.fill(HIST("Gen/Charged/h_PtTruth"), pt); hist.fill(HIST("Gen/Charged/h_PtTruth_e"), pt); @@ -1385,25 +1515,25 @@ struct MeanPtFlucId { rap = mcPart.y(); auto pid = mcPart.pdgCode(); - if (std::abs(pid) == kPiPlus && pt >= cfgCutPiPtMin) { + if (std::abs(pid) == kPiPlus && pt >= pidConfig.cfgCutPiPtMin) { pionTracksGen.push_back({pt}); fillPtMCHist(true, pt, eta, rap, phi, pid, kPiPlus, kPiMinus); } - if (std::abs(pid) == kKPlus && pt >= cfgCutKaPtMin) { + if (std::abs(pid) == kKPlus && pt >= pidConfig.cfgCutKaPtMin) { kaonTracksGen.push_back({pt}); fillPtMCHist(true, pt, eta, rap, phi, pid, kKPlus, kKMinus); } - if (std::abs(pid) == kProton && pt >= cfgCutPrPtMin) { + if (std::abs(pid) == kProton && pt >= pidConfig.cfgCutPrPtMin) { protonTracksGen.push_back({pt}); fillPtMCHist(true, pt, eta, rap, phi, pid, kProton, kProtonBar); } } } - fillMatrices(chargedTracksGen, cfgCutPtMin); - fillMatrices(pionTracksGen, cfgCutPiPtMin); - fillMatrices(kaonTracksGen, cfgCutKaPtMin); - fillMatrices(protonTracksGen, cfgCutPrPtMin); + fillMatrices(chargedTracksGen, trackCuts.cfgCutPtMin); + fillMatrices(pionTracksGen, pidConfig.cfgCutPiPtMin); + fillMatrices(kaonTracksGen, pidConfig.cfgCutKaPtMin); + fillMatrices(protonTracksGen, pidConfig.cfgCutPrPtMin); hist.fill(HIST("Gen/h_Counts"), 2); hist.fill(HIST("Gen/h_VtxZ"), mcCol.posZ()); @@ -1451,14 +1581,14 @@ struct MeanPtFlucId { aod::pidTOFFullPi, aod::pidTPCFullPi, aod::pidTOFFullPr, aod::pidTPCFullPr, aod::pidTOFFullKa, aod::pidTPCFullKa, aod::pidTOFFullEl, aod::pidTPCFullEl, aod::pidTOFbeta, aod::pidTOFmass>; - using MyRun3Collisions = soa::Join; - using MyRun3MCCollisions = soa::Join; + using MyRun3Collisions = soa::Join; + using MyRun3MCCollisions = soa::Join; using MyMCTracks = soa::Join; SliceCache cache; - Preslice perCollision = aod::track::collisionId; + PresliceOptional perCollision = aod::track::collisionId; void processRun3(MyRun3Collisions::iterator const& col, MyAllTracks const& tracks) { @@ -1472,6 +1602,115 @@ struct MeanPtFlucId { analyzeMC(mcCol, cols, tracks, mcParts); } PROCESS_SWITCH(MeanPtFlucId, processMCRecoSimRun3, "process MC Reconstructed & Truth Run-3", true); + + template + int countGenForwardMultiplicity(M const& particles) + { + constexpr double MinimumChargedParticleCharge = 1.e-3; + constexpr double Ft0aEtaMin = 3.5; + constexpr double Ft0aEtaMax = 4.9; + constexpr double Ft0cEtaMin = -3.3; + constexpr double Ft0cEtaMax = -2.1; + constexpr double Fv0aEtaMin = 2.2; + constexpr double Fv0aEtaMax = 5.0; + int ft0a = 0, ft0c = 0, fv0a = 0; + // Match Common/Tools/Multiplicity/MultModule.h::collisionProcessMonteCarlo. + // No analysis pT/eta cuts: these are forward charged physical primaries. + for (auto const& particle : particles) { + if (!particle.isPhysicalPrimary()) { + continue; + } + auto* definition = pdg->GetParticle(particle.pdgCode()); + if (!definition || std::abs(definition->Charge()) < MinimumChargedParticleCharge) { + continue; + } + const double eta = particle.eta(); + ft0a += (Ft0aEtaMin < eta && eta < Ft0aEtaMax); + ft0c += (Ft0cEtaMin < eta && eta < Ft0cEtaMax); + fv0a += (Fv0aEtaMin < eta && eta < Fv0aEtaMax); + } + return eventCuts.cfgCentFT0C ? ft0c : (eventCuts.cfgCentFV0A ? fv0a : ft0a + ft0c); + } + + double getGenCentrality(int forwardMultiplicity) const + { + const auto& lookup = genConfig.cfgGenCentralityLookup.value; + if (forwardMultiplicity < 0 || static_cast(forwardMultiplicity) >= lookup.size()) { + return -1.; // Do not extrapolate beyond the reference distribution. + } + return std::min(lookup[forwardMultiplicity], std::nextafter(100., 0.)); + } + + // Only MC tables are inputs: no reconstructed event, track, centrality or PID selection. + // The framework groups mcParts by the current mcCol through mcCollisionId. + void processMCGenOnly(aod::McCollisions::iterator const& mcCol, aod::McParticles const& mcParts) + { + const double eventWeight = genConfig.cfgGenUseEventWeight ? mcCol.weight() : 1.; + if (!std::isfinite(eventWeight) || eventWeight <= 0.) { + LOGF(fatal, "Generator percentile analysis requires finite positive event weights"); + } + hist.fill(HIST("Model/h_EventCount"), 1, eventWeight); + if (genConfig.cfgGenCutPosZ && std::abs(mcCol.posZ()) > genConfig.cfgGenMaxPosZ) { + return; + } + hist.fill(HIST("Model/h_EventCount"), 2, eventWeight); + hist.fill(HIST("Model/h_VtxZ"), mcCol.posZ(), eventWeight); + const int forwardMultiplicity = countGenForwardMultiplicity(mcParts); + hist.fill(HIST("Model/h_ForwardMult"), forwardMultiplicity, eventWeight); + if (genConfig.cfgGenCentralityCalibration) { + return; + } + const double centrality = getGenCentrality(forwardMultiplicity); + if (centrality < 0.) { + hist.fill(HIST("Model/h_CentralityStatus"), 2, eventWeight); + return; + } + hist.fill(HIST("Model/h_CentralityStatus"), 1, eventWeight); + hist.fill(HIST("Model/h_Centrality"), centrality, eventWeight); + const int multiplicity = countSimMultiplicity(mcParts); + hist.fill(HIST("Model/h_Nch"), multiplicity, eventWeight); + + std::array, 4> modelTracks; + for (auto const& particle : mcParts) { + if (!particle.isPhysicalPrimary()) { + continue; + } + auto* definition = pdg->GetParticle(particle.pdgCode()); + if (!definition || definition->Charge() == 0.) { + continue; + } + const float pt = particle.pt(); + const float eta = particle.eta(); + if (!std::isfinite(pt) || !std::isfinite(eta) || + pt < trackCuts.cfgCutPtMin || pt >= trackCuts.cfgCutPtMax || std::abs(eta) >= trackCuts.cfgCutEta) { + continue; + } + modelTracks[0].push_back({pt}); + hist.fill(HIST("Model/Charged/h_Pt"), pt, eventWeight); + hist.fill(HIST("Model/Charged/h_Eta"), eta, eventWeight); + const int code = std::abs(particle.pdgCode()); + if (code == kPiPlus && pt >= pidConfig.cfgCutPiPtMin) { + modelTracks[1].push_back({pt}); + hist.fill(HIST("Model/Pion/h_Pt"), pt, eventWeight); + hist.fill(HIST("Model/Pion/h_Eta"), eta, eventWeight); + } + if (code == kKPlus && pt >= pidConfig.cfgCutKaPtMin) { + modelTracks[2].push_back({pt}); + hist.fill(HIST("Model/Kaon/h_Pt"), pt, eventWeight); + hist.fill(HIST("Model/Kaon/h_Eta"), eta, eventWeight); + } + if (code == kProton && pt >= pidConfig.cfgCutPrPtMin) { + modelTracks[3].push_back({pt}); + hist.fill(HIST("Model/Proton/h_Pt"), pt, eventWeight); + hist.fill(HIST("Model/Proton/h_Eta"), eta, eventWeight); + } + } + fillMatrices(modelTracks[0], trackCuts.cfgCutPtMin, centrality, eventWeight); + fillMatrices(modelTracks[1], pidConfig.cfgCutPiPtMin, centrality, eventWeight); + fillMatrices(modelTracks[2], pidConfig.cfgCutKaPtMin, centrality, eventWeight); + fillMatrices(modelTracks[3], pidConfig.cfgCutPrPtMin, centrality, eventWeight); + } + PROCESS_SWITCH(MeanPtFlucId, processMCGenOnly, "Independent generator-only model comparison", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& context)