From 1963f8cc6e08119cb2b536f321596c1176cc0c8c Mon Sep 17 00:00:00 2001 From: Emil Gorm Nielsen Date: Wed, 12 Aug 2026 13:23:28 +0200 Subject: [PATCH 1/4] Add Nch MC response matrix --- .../GenericFramework/Tasks/flowGfwNonflow.cxx | 64 +++++++++++++++++-- 1 file changed, 60 insertions(+), 4 deletions(-) diff --git a/PWGCF/GenericFramework/Tasks/flowGfwNonflow.cxx b/PWGCF/GenericFramework/Tasks/flowGfwNonflow.cxx index 949724ac0ea..ebbee2546b9 100644 --- a/PWGCF/GenericFramework/Tasks/flowGfwNonflow.cxx +++ b/PWGCF/GenericFramework/Tasks/flowGfwNonflow.cxx @@ -78,7 +78,7 @@ struct FlowGfwNonflow { Configurable cfgMpar{"cfgMpar", 4, "Highest order of pt-pt correlations"}; Configurable cfgCentEstimator{"cfgCentEstimator", 0, "0:FT0C; 1:FT0CVariant1; 2:FT0M; 3:FT0A, 4:NTPV, 5:NGlobal, 6:MFT"}; Configurable cfgUseNch{"cfgUseNch", false, "Do correlations as function of Nch"}; - Configurable cfgUseNchCorrection{"cfgUseNchCorrection", 1, "Use correction for Nch; 0: Use size of tracks table, 1: Use efficiency-corrected Nch values, 2: Use uncorrected Nch values"}; + Configurable cfgUseNchCorrection{"cfgUseNchCorrection", 1, "Nch used on the x-axis; 0: tracks table size, 1: efficiency-corrected, 2: accepted reconstructed, 3: response-matrix corrected"}; Configurable cfgRunByRun{"cfgRunByRun", false, "Use run-by-run NUA"}; Configurable cfgFillQA{"cfgFillQA", false, "Fill QA histograms"}; Configurable cfgUseCentralMoments{"cfgUseCentralMoments", true, "Use central moments in vn-pt calculations"}; @@ -86,6 +86,7 @@ struct FlowGfwNonflow { struct : ConfigurableGroup { Configurable cfgEfficiencyPath{"cfgEfficiencyPath", "", "CCDB path to efficiency object"}; Configurable cfgUse2DEfficiency{"cfgUse2DEfficiency", false, "Toggle the use of 2D (pt, centrality) efficiency versus centrality integrated efficiency"}; + Configurable cfgNchResponsePath{"cfgNchResponsePath", "", "CCDB path to TH2 response matrix (reconstructed Nch on x, generated Nch on y)"}; Configurable cfgAcceptancePath{"cfgAcceptancePath", "", "CCDB path to acceptance object"}; } cfgCorrections; struct : ConfigurableGroup { @@ -181,6 +182,7 @@ struct FlowGfwNonflow { struct Config { TH1* mEfficiency = nullptr; + TH2* mNchResponse = nullptr; std::vector mAcceptance; bool correctionsLoaded = false; } correctionsConfig; @@ -379,6 +381,9 @@ struct FlowGfwNonflow { registry.add("eventQA/before/occ_mult_cent", "; occupancy; N_{ch}; centrality (%)", {HistType::kTH3D, {occAxis, nchAxis, centAxis}}); } } + if (doprocessMCReco) { + registry.add("MCReco/Nch_reco_gen", "; N_{ch}^{reco}; N_{ch}^{gen}", {HistType::kTH2D, {nchAxis, nchAxis}}); + } registry.add("eventQA/before/centrality", "; centrality (%); Counts", {HistType::kTH1D, {centAxis}}); registry.add("eventQA/before/multiplicity", "; N_{ch}; Counts", {HistType::kTH1D, {nchAxis}}); registry.addClone("eventQA/before/", "eventQA/after/"); @@ -639,9 +644,43 @@ struct FlowGfwNonflow { } LOGF(info, "Loaded efficiency histogram from %s", cfgCorrections.cfgEfficiencyPath.value.c_str()); } + if (!cfgCorrections.cfgNchResponsePath.value.empty()) { + correctionsConfig.mNchResponse = ccdb->getForTimeStamp(cfgCorrections.cfgNchResponsePath, timestamp); + if (correctionsConfig.mNchResponse == nullptr) { + LOGF(fatal, "Could not load Nch response matrix from %s", cfgCorrections.cfgNchResponsePath.value.c_str()); + } + LOGF(info, "Loaded Nch response matrix from %s", cfgCorrections.cfgNchResponsePath.value.c_str()); + } else if (cfgUseNchCorrection == 3) { + LOGF(fatal, "cfgUseNchCorrection=3 requires cfgNchResponsePath"); + } correctionsConfig.correctionsLoaded = true; } + float getResponseCorrectedNch(const unsigned int multReconstructed) const + { + if (!correctionsConfig.mNchResponse) { + return multReconstructed; + } + const auto* response = correctionsConfig.mNchResponse; + const int recoBin = response->GetXaxis()->FindFixBin(multReconstructed); + if (recoBin < 1 || recoBin > response->GetNbinsX()) { + LOGF(warn, "Reconstructed Nch %u is outside the response matrix; using the uncorrected value", multReconstructed); + return reconstructedNch; + } + double sumWeights = 0.; + double sumGeneratedNch = 0.; + for (int genBin = 1; genBin <= response->GetNbinsY(); ++genBin) { + const double weight = response->GetBinContent(recoBin, genBin); + sumWeights += weight; + sumGeneratedNch += weight * response->GetYaxis()->GetBinCenter(genBin); + } + if (sumWeights <= 0.) { + LOGF(warn, "Response matrix has no entries for reconstructed Nch %u; using the uncorrected value", multReconstructed); + return multReconstructed; + } + return sumGeneratedNch / sumWeights; + } + template double getAcceptance(const TTrack& track, const double& vtxz) { // 0 ref, 1 ch, 2 pi, 3 ka, 4 pr @@ -1014,7 +1053,7 @@ struct FlowGfwNonflow { }; template - void processCollision(const TCollision& collision, const TTracks& tracks, const float& centrality, const float& field) + void processCollision(const TCollision& collision, const TTracks& tracks, const float& centrality, const float& field, const int generatedNch = -1) { if (tracks.size() < 1) { return; @@ -1037,6 +1076,11 @@ struct FlowGfwNonflow { for (const auto& track : tracks) { processTrack(track, vtxz, field, centrality, acceptedTracks); } + if constexpr (dt == Reco) { + if (generatedNch >= 0) { + registry.fill(HIST("MCReco/Nch_reco_gen"), acceptedTracks.totaluncorr, generatedNch); + } + } if (dt != Gen && cfgFillQA) { registry.fill(HIST("trackQA/after/Nch_corrected"), acceptedTracks.total); registry.fill(HIST("trackQA/after/Nch_uncorrected"), acceptedTracks.totaluncorr); @@ -1053,6 +1097,9 @@ struct FlowGfwNonflow { case 2: multiplicity = acceptedTracks.totaluncorr; break; + case 3: + multiplicity = (dt == Gen) ? acceptedTracks.totaluncorr : getResponseCorrectedNch(acceptedTracks.totaluncorr); + break; default: multiplicity = tracks.size(); break; @@ -1155,8 +1202,10 @@ struct FlowGfwNonflow { using GFWCollisions = soa::Filtered>; using GFWMCCollisions = soa::Join; + using FilteredGFWMCCollisions = soa::Filtered; using GFWTracks = soa::Filtered>; using GFWMCTracks = soa::Filtered>; + Preslice particlesPerMcCollision = aod::mcparticle::mcCollisionId; SliceCache cache; Partition posTracks = aod::track::signed1Pt > 0.0f; @@ -1221,7 +1270,7 @@ struct FlowGfwNonflow { } PROCESS_SWITCH(FlowGfwNonflow, processData, "Process analysis for non-derived data", true); - void processMCReco(GFWCollisions::iterator const& collision, aod::BCsWithTimestamps const&, GFWMCTracks const& tracks, aod::McParticles const&) + void processMCReco(FilteredGFWMCCollisions::iterator const& collision, aod::BCsWithTimestamps const&, GFWMCTracks const& tracks, aod::McParticles const& particles) { auto bc = collision.bc_as(); int run = bc.runNumber(); @@ -1265,9 +1314,16 @@ struct FlowGfwNonflow { if (cfgFillQA) { fillEventQA(collision, tracks); } + unsigned int generatedNch = 0; + const auto particlesThisCollision = particles.sliceBy(particlesPerMcCollision, collision.mcCollisionId()); + for (const auto& particle : particlesThisCollision) { + if (particle.isPhysicalPrimary() && particle.eta() > cfgKinematics.cfgEtaNch->first && particle.eta() < cfgKinematics.cfgEtaNch->second && particle.pt() > gfwMemberCache.ptlow && particle.pt() < gfwMemberCache.ptup) { + ++generatedNch; + } + } loadCorrections(bc); auto field = (cfgEventSelection.cfgMagField == DefaultMagneticFieldCut) ? getMagneticField(bc.timestamp()) : static_cast(cfgEventSelection.cfgMagField); - processCollision(collision, tracks, centrality, field); + processCollision(collision, tracks, centrality, field, generatedNch); } PROCESS_SWITCH(FlowGfwNonflow, processMCReco, "Process analysis for MC reconstructed events", false); From 0c2a24a6b40aed848411ac7e3070011828e1ea17 Mon Sep 17 00:00:00 2001 From: Emil Gorm Nielsen Date: Mon, 17 Aug 2026 09:49:47 +0200 Subject: [PATCH 2/4] Add guard conditions on histogram filling --- .../Tasks/flowGenericFramework.cxx | 34 +++++++++++++------ .../GenericFramework/Tasks/flowGfwNonflow.cxx | 18 ++++++---- 2 files changed, 36 insertions(+), 16 deletions(-) diff --git a/PWGCF/GenericFramework/Tasks/flowGenericFramework.cxx b/PWGCF/GenericFramework/Tasks/flowGenericFramework.cxx index 77ca0b782c5..710e892d87a 100644 --- a/PWGCF/GenericFramework/Tasks/flowGenericFramework.cxx +++ b/PWGCF/GenericFramework/Tasks/flowGenericFramework.cxx @@ -686,8 +686,6 @@ struct FlowGenericFramework { registryQA.add("trackQA/after/Nch_uncorrected", "; N_{ch}; Counts", {HistType::kTH1D, {nchAxis}}); registryQA.add("trackQA/after/etaNch", "; #eta; Counts", {HistType::kTH1D, {etaAxis}}); registryQA.add("trackQA/after/etaPtPt", "; #eta; Counts", {HistType::kTH1D, {etaAxis}}); - registryQA.add("trackQA/after/etaV02", "; #eta; Counts", {HistType::kTH1D, {etaAxis}}); - registryQA.add("trackQA/after/etaV0", "; #eta; Counts", {HistType::kTH1D, {etaAxis}}); if (!cfgFill.cfgFillRunByRunQA) { if (cfgUsePID) { registryQA.add("phi_eta_vtxz_ref", "", {HistType::kTH3D, {phiAxis, etaAxis, vtxAxis}}); @@ -732,6 +730,10 @@ struct FlowGenericFramework { AxisSpec axisLambdaMass = {resoSwitchVals[MassBins][Lambda], resoCutVals[MassMin][Lambda], resoCutVals[MassMax][Lambda]}; AxisSpec yAxis = {100, -1, 1}; // QA histograms for V0s + if (cfgFill.cfgFillV0QA && (resoSwitchVals[UseParticle][K0] != 0 || resoSwitchVals[UseParticle][Lambda] != 0)) { + registryQA.add("trackQA/after/etaV02", "; #eta; Counts", {HistType::kTH1D, {etaAxis}}); + registryQA.add("trackQA/after/etaV0", "; #eta; Counts", {HistType::kTH1D, {etaAxis}}); + } if (resoSwitchVals[UseParticle][K0] != 0) { if (cfgFill.cfgFillV0QA) { registryQA.add("K0/PiPlusTPC_K0", "", {HistType::kTH2D, {{ptAxis, axisNsigmaTPC}}}); @@ -2373,7 +2375,9 @@ struct FlowGenericFramework { registryQA.fill(HIST("K0/hK0Count"), FillDaughterTrackSelected); selection.selected = true; selection.isK0 = true; - registryQA.fill(HIST("K0/hK0AP"), v0.alpha(), v0.qtarm()); + if (cfgFill.cfgFillV0QA) { + registryQA.fill(HIST("K0/hK0AP"), v0.alpha(), v0.qtarm()); + } return selection; } @@ -2481,7 +2485,7 @@ struct FlowGenericFramework { if (!selectionV0Daughter(postrack, Protons) || !selectionV0Daughter(negtrack, Pions)) { return selection; } - if (fillSelectionQA) { + if (fillSelectionQA && cfgFill.cfgFillV0QA) { registryQA.fill(HIST("Lambda/hLambdaAP"), v0.alpha(), v0.qtarm()); } } @@ -2489,7 +2493,7 @@ struct FlowGenericFramework { if (!selectionV0Daughter(postrack, Pions) || !selectionV0Daughter(negtrack, Protons)) { return selection; } - if (fillSelectionQA) { + if (fillSelectionQA && cfgFill.cfgFillV0QA) { registryQA.fill(HIST("Lambda/hAntiLambdaAP"), v0.alpha(), v0.qtarm()); } } @@ -2827,11 +2831,15 @@ struct FlowGenericFramework { if (cfgEventSelection.cfgDoOccupancySel) { int occupancy = collision.trackOccupancyInTimeRange(); - registryQA.fill(HIST("eventQA/before/occ_mult_cent"), occupancy, tracks.size(), centrality); + if (cfgFill.cfgFillQA) { + registryQA.fill(HIST("eventQA/before/occ_mult_cent"), occupancy, tracks.size(), centrality); + } if (occupancy < 0 || occupancy > cfgEventSelection.cfgOccupancySelection) { return; } - registryQA.fill(HIST("eventQA/after/occ_mult_cent"), occupancy, tracks.size(), centrality); + if (cfgFill.cfgFillQA) { + registryQA.fill(HIST("eventQA/after/occ_mult_cent"), occupancy, tracks.size(), centrality); + } } registryQA.fill(HIST("eventQA/eventSel"), 2.5); if (cfgFill.cfgFillRunByRunQA) { @@ -2874,7 +2882,9 @@ struct FlowGenericFramework { void processOnTheFly(soa::Filtered::iterator const& mcCollision, aod::McParticles const& mcParticles, aod::V0Datas const& v0s) { int run = 0; - registryQA.fill(HIST("MCGen/impactParameter"), mcCollision.impactParameter(), mcParticles.size()); + if (cfgFill.cfgFillQA) { + registryQA.fill(HIST("MCGen/impactParameter"), mcCollision.impactParameter(), mcParticles.size()); + } processCollision(mcCollision, mcParticles, v0s, mcCollision.impactParameter(), -999, run); } PROCESS_SWITCH(FlowGenericFramework, processOnTheFly, "Process analysis for MC on-the-fly generated events", false); @@ -3269,7 +3279,9 @@ struct FlowGenericFramework { continue; } fillGeneratedEfficiencyTrack(particle, selectedCentrality); - fillGeneratedLambdaFeeddownXi(particle, selectedCentrality); + if (cfgFill.cfgFillV0QA) { + fillGeneratedLambdaFeeddownXi(particle, selectedCentrality); + } if (isGeneratedEfficiencyV0(particle, PDG_t::kK0Short, K0) && resoSwitchVals[UseParticle][K0] != 0) { fillGeneratedEfficiencyV0(particle, EfficiencyK0, selectedCentrality); } @@ -3294,7 +3306,9 @@ struct FlowGenericFramework { if (v0.collisionId() != bestCollisionIndex) { continue; } - fillLambdaFeeddownReco(v0, collision, tracks, selectedCentrality); + if (cfgFill.cfgFillV0QA) { + fillLambdaFeeddownReco(v0, collision, tracks, selectedCentrality); + } fillEfficiencyRecoV0(v0, collision, tracks, selectedCentrality); } break; diff --git a/PWGCF/GenericFramework/Tasks/flowGfwNonflow.cxx b/PWGCF/GenericFramework/Tasks/flowGfwNonflow.cxx index ebbee2546b9..b904f9e938f 100644 --- a/PWGCF/GenericFramework/Tasks/flowGfwNonflow.cxx +++ b/PWGCF/GenericFramework/Tasks/flowGfwNonflow.cxx @@ -242,6 +242,12 @@ struct FlowGfwNonflow { ProtonID, SpeciesCount }; + enum NchSelector { + TableSize, + Corrected, + Uncorrected, + ResponseMatrixCorrected + }; // Generic Framework GFW* fGFW = new GFW(); @@ -650,7 +656,7 @@ struct FlowGfwNonflow { LOGF(fatal, "Could not load Nch response matrix from %s", cfgCorrections.cfgNchResponsePath.value.c_str()); } LOGF(info, "Loaded Nch response matrix from %s", cfgCorrections.cfgNchResponsePath.value.c_str()); - } else if (cfgUseNchCorrection == 3) { + } else if (cfgUseNchCorrection == NchSelector::ResponseMatrixCorrected) { LOGF(fatal, "cfgUseNchCorrection=3 requires cfgNchResponsePath"); } correctionsConfig.correctionsLoaded = true; @@ -665,7 +671,7 @@ struct FlowGfwNonflow { const int recoBin = response->GetXaxis()->FindFixBin(multReconstructed); if (recoBin < 1 || recoBin > response->GetNbinsX()) { LOGF(warn, "Reconstructed Nch %u is outside the response matrix; using the uncorrected value", multReconstructed); - return reconstructedNch; + return multReconstructed; } double sumWeights = 0.; double sumGeneratedNch = 0.; @@ -1088,16 +1094,16 @@ struct FlowGfwNonflow { float multiplicity = 0.f; switch (cfgUseNchCorrection) { - case 0: + case NchSelector::TableSize: multiplicity = tracks.size(); break; - case 1: + case NchSelector::Corrected: multiplicity = acceptedTracks.total; break; - case 2: + case NchSelector::Uncorrected: multiplicity = acceptedTracks.totaluncorr; break; - case 3: + case NchSelector::ResponseMatrixCorrected: multiplicity = (dt == Gen) ? acceptedTracks.totaluncorr : getResponseCorrectedNch(acceptedTracks.totaluncorr); break; default: From 397df7cd22c84a2508795b86631e02e894220c38 Mon Sep 17 00:00:00 2001 From: Emil Gorm Nielsen Date: Mon, 17 Aug 2026 12:46:45 +0200 Subject: [PATCH 3/4] codecheck --- .../GenericFramework/Tasks/flowGfwNonflow.cxx | 31 ++++++++++--------- 1 file changed, 16 insertions(+), 15 deletions(-) diff --git a/PWGCF/GenericFramework/Tasks/flowGfwNonflow.cxx b/PWGCF/GenericFramework/Tasks/flowGfwNonflow.cxx index b904f9e938f..288492d4347 100644 --- a/PWGCF/GenericFramework/Tasks/flowGfwNonflow.cxx +++ b/PWGCF/GenericFramework/Tasks/flowGfwNonflow.cxx @@ -44,6 +44,7 @@ #include #include +#include #include #include #include @@ -722,22 +723,22 @@ struct FlowGfwNonflow { return -1.; } return 1. / eff; - } else { - auto* effHist = dynamic_cast(correctionsConfig.mEfficiency); - if (!effHist) { - LOGF(error, "Efficiency object at %s is not a TH1D", cfgCorrections.cfgEfficiencyPath.value.c_str()); - return -1.; - } - bin = effHist->FindBin(track.pt()); - if (!bin) { - return -1.; - } - const double eff = effHist->GetBinContent(bin); - if (!std::isfinite(eff) || eff <= 0.) { - return -1.; - } - return 1. / eff; } + + auto* effHist = dynamic_cast(correctionsConfig.mEfficiency); + if (!effHist) { + LOGF(error, "Efficiency object at %s is not a TH1D", cfgCorrections.cfgEfficiencyPath.value.c_str()); + return -1.; + } + bin = effHist->FindBin(track.pt()); + if (!bin) { + return -1.; + } + const double eff = effHist->GetBinContent(bin); + if (!std::isfinite(eff) || eff <= 0.) { + return -1.; + } + return 1. / eff; } template From 176caaee34d2c944ee917e174c765368d83904c2 Mon Sep 17 00:00:00 2001 From: Emil Gorm Nielsen Date: Mon, 17 Aug 2026 22:27:38 +0200 Subject: [PATCH 4/4] force rebuild --- PWGCF/GenericFramework/Tasks/flowGfwNonflow.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PWGCF/GenericFramework/Tasks/flowGfwNonflow.cxx b/PWGCF/GenericFramework/Tasks/flowGfwNonflow.cxx index 288492d4347..cdccf5bd2fd 100644 --- a/PWGCF/GenericFramework/Tasks/flowGfwNonflow.cxx +++ b/PWGCF/GenericFramework/Tasks/flowGfwNonflow.cxx @@ -79,7 +79,7 @@ struct FlowGfwNonflow { Configurable cfgMpar{"cfgMpar", 4, "Highest order of pt-pt correlations"}; Configurable cfgCentEstimator{"cfgCentEstimator", 0, "0:FT0C; 1:FT0CVariant1; 2:FT0M; 3:FT0A, 4:NTPV, 5:NGlobal, 6:MFT"}; Configurable cfgUseNch{"cfgUseNch", false, "Do correlations as function of Nch"}; - Configurable cfgUseNchCorrection{"cfgUseNchCorrection", 1, "Nch used on the x-axis; 0: tracks table size, 1: efficiency-corrected, 2: accepted reconstructed, 3: response-matrix corrected"}; + Configurable cfgUseNchCorrection{"cfgUseNchCorrection", 1, "Nch used on the x-axis; 0: tracks table size, 1: efficiency-corrected, 2: accepted reconstructed, 3: Reco.-gen. response-matrix corrected"}; Configurable cfgRunByRun{"cfgRunByRun", false, "Use run-by-run NUA"}; Configurable cfgFillQA{"cfgFillQA", false, "Fill QA histograms"}; Configurable cfgUseCentralMoments{"cfgUseCentralMoments", true, "Use central moments in vn-pt calculations"};