From 01dc96838c867478600b8837056efec322df39ea Mon Sep 17 00:00:00 2001 From: abmodak <67369858+abmodak@users.noreply.github.com> Date: Mon, 6 Oct 2025 20:14:03 +0200 Subject: [PATCH 1/5] Add multiplicity classifier --- .../Tasks/longrangeCorrelation.cxx | 181 ++++++++++++++---- 1 file changed, 141 insertions(+), 40 deletions(-) diff --git a/PWGCF/TwoParticleCorrelations/Tasks/longrangeCorrelation.cxx b/PWGCF/TwoParticleCorrelations/Tasks/longrangeCorrelation.cxx index a1eb8e9a3eb..bc5949ba294 100644 --- a/PWGCF/TwoParticleCorrelations/Tasks/longrangeCorrelation.cxx +++ b/PWGCF/TwoParticleCorrelations/Tasks/longrangeCorrelation.cxx @@ -93,6 +93,17 @@ static constexpr std::string_view kEvntType[] = {"SE/", "ME/"}; auto static constexpr kMinFt0cCell = 96; AxisSpec axisEvent{10, 0.5, 9.5, "#Event", "EventAxis"}; +namespace o2::aod +{ +namespace longrangemultclass +{ +DECLARE_SOA_COLUMN(Multiplicity, multiplicity, float); //! Centrality/multiplicity value +} // namespace longrangemultclass +DECLARE_SOA_TABLE(LRMultTables, "AOD", "LRMULTTABLE", longrangemultclass::Multiplicity); //! Transient multiplicity table + +using LRMultTable = LRMultTables::iterator; +} // namespace o2::aod + struct LongrangeCorrelation { struct : ConfigurableGroup { @@ -139,7 +150,7 @@ struct LongrangeCorrelation { ConfigurableAxis amplitudeFt0a{"amplitudeFt0a", {5000, 0, 10000}, "FT0A amplitude"}; ConfigurableAxis channelFt0aAxis{"channelFt0aAxis", {96, 0.0, 96.0}, "FT0A channel"}; - using CollTable = soa::Join; + using CollTable = soa::Join; using TrksTable = soa::Filtered>; using MftTrkTable = soa::Filtered; Preslice perColGlobal = aod::track::collisionId; @@ -325,11 +336,11 @@ struct LongrangeCorrelation { } template - void fillCorrFt0aGlobal(TTarget target, TTriggers const& triggers, TFT0s const& ft0, bool mixing, float vz) + void fillCorrFt0aGlobal(TTarget target, TTriggers const& triggers, TFT0s const& ft0, bool mixing, float vz, float multiplicity) { int fSampleIndex = gRandom->Uniform(0, cfgSampleSize); if (!mixing) - histos.fill(HIST("Ft0aGlobal/SE/hMult_used"), triggers.size()); + histos.fill(HIST("Ft0aGlobal/SE/hMult_used"), multiplicity); for (auto const& triggerTrack : triggers) { if (!mixing) histos.fill(HIST("Ft0aGlobal/SE/Trig_hist"), fSampleIndex, vz, triggerTrack.pt()); @@ -367,11 +378,11 @@ struct LongrangeCorrelation { } // fillCorrFt0aGlobal template - void fillCorrFt0cGlobal(TTarget target, TTriggers const& triggers, TFT0s const& ft0, bool mixing, float vz) + void fillCorrFt0cGlobal(TTarget target, TTriggers const& triggers, TFT0s const& ft0, bool mixing, float vz, float multiplicity) { int fSampleIndex = gRandom->Uniform(0, cfgSampleSize); if (!mixing) - histos.fill(HIST("Ft0cGlobal/SE/hMult_used"), triggers.size()); + histos.fill(HIST("Ft0cGlobal/SE/hMult_used"), multiplicity); for (auto const& triggerTrack : triggers) { if (!mixing) histos.fill(HIST("Ft0cGlobal/SE/Trig_hist"), fSampleIndex, vz, triggerTrack.pt()); @@ -409,11 +420,11 @@ struct LongrangeCorrelation { } // fillCorrFt0cGlobal template - void fillCorrMftGlobal(TTarget target, TTriggers const& triggers, TMFTs const& mft, bool mixing, float vz) + void fillCorrMftGlobal(TTarget target, TTriggers const& triggers, TMFTs const& mft, bool mixing, float vz, float multiplicity) { int fSampleIndex = gRandom->Uniform(0, cfgSampleSize); if (!mixing) - histos.fill(HIST("MftGlobal/SE/hMult_used"), triggers.size()); + histos.fill(HIST("MftGlobal/SE/hMult_used"), multiplicity); for (auto const& triggerTrack : triggers) { if (!mixing) histos.fill(HIST("MftGlobal/SE/Trig_hist"), fSampleIndex, vz, triggerTrack.pt()); @@ -444,12 +455,12 @@ struct LongrangeCorrelation { } // trigger tracks } // fillCorrMftGlobal - template - void fillCorrFt0aMft(TTarget target, TTracks const& tracks, TTriggers const& triggers, TFT0s const& ft0, bool mixing, float vz) + template + void fillCorrFt0aMft(TTarget target, TTriggers const& triggers, TFT0s const& ft0, bool mixing, float vz, float multiplicity) { int fSampleIndex = gRandom->Uniform(0, cfgSampleSize); if (!mixing) - histos.fill(HIST("Ft0aMft/SE/hMult_used"), tracks.size()); + histos.fill(HIST("Ft0aMft/SE/hMult_used"), multiplicity); for (auto const& triggerTrack : triggers) { if (!isMftTrackSelected(triggerTrack)) { continue; @@ -493,12 +504,12 @@ struct LongrangeCorrelation { } // trigger tracks } // fillCorrFt0aMft - template - void fillCorrFt0aFt0c(TTarget target, TTriggers const& triggers, TFT0As const& ft0a, TFT0Cs const& ft0c, bool mixing, float vz) + template + void fillCorrFt0aFt0c(TTarget target, TFT0As const& ft0a, TFT0Cs const& ft0c, bool mixing, float vz, float multiplicity) { int fSampleIndex = gRandom->Uniform(0, cfgSampleSize); if (!mixing) - histos.fill(HIST("Ft0aFt0c/SE/hMult_used"), triggers.size()); + histos.fill(HIST("Ft0aFt0c/SE/hMult_used"), multiplicity); for (std::size_t iChA = 0; iChA < ft0a.channelA().size(); iChA++) { if (!mixing) @@ -581,10 +592,11 @@ struct LongrangeCorrelation { if (col.has_foundFT0()) { fillYield(tracks); const auto& ft0 = col.foundFT0(); - if (tracks.size() < cfgMinMult || tracks.size() >= cfgMaxMult) { + auto multiplicity = col.multiplicity(); + if (multiplicity < cfgMinMult || multiplicity >= cfgMaxMult) { return; } - fillCorrFt0aGlobal(sameFt0aGlobal, tracks, ft0, false, col.posZ()); + fillCorrFt0aGlobal(sameFt0aGlobal, tracks, ft0, false, col.posZ(), multiplicity); } } // same event @@ -596,10 +608,11 @@ struct LongrangeCorrelation { if (col.has_foundFT0()) { fillYield(tracks); const auto& ft0 = col.foundFT0(); - if (tracks.size() < cfgMinMult || tracks.size() >= cfgMaxMult) { + auto multiplicity = col.multiplicity(); + if (multiplicity < cfgMinMult || multiplicity >= cfgMaxMult) { return; } - fillCorrFt0cGlobal(sameFt0cGlobal, tracks, ft0, false, col.posZ()); + fillCorrFt0cGlobal(sameFt0cGlobal, tracks, ft0, false, col.posZ(), multiplicity); } } // same event @@ -609,13 +622,14 @@ struct LongrangeCorrelation { return; } fillYield(tracks); - if (tracks.size() < cfgMinMult || tracks.size() >= cfgMaxMult) { + auto multiplicity = col.multiplicity(); + if (multiplicity < cfgMinMult || multiplicity >= cfgMaxMult) { return; } - fillCorrMftGlobal(sameMftGlobal, tracks, mfttracks, false, col.posZ()); + fillCorrMftGlobal(sameMftGlobal, tracks, mfttracks, false, col.posZ(), multiplicity); } // same event - void processFt0aMftSE(CollTable::iterator const& col, aod::FT0s const&, TrksTable const& tracks, MftTrkTable const& mfttracks) + void processFt0aMftSE(CollTable::iterator const& col, aod::FT0s const&, MftTrkTable const& mfttracks) { if (!isEventSelected(col)) { return; @@ -623,25 +637,27 @@ struct LongrangeCorrelation { if (col.has_foundFT0()) { fillYield(mfttracks); const auto& ft0 = col.foundFT0(); - if (tracks.size() < cfgMinMult || tracks.size() >= cfgMaxMult) { + auto multiplicity = col.multiplicity(); + if (multiplicity < cfgMinMult || multiplicity >= cfgMaxMult) { return; } - fillCorrFt0aMft(sameFt0aMft, tracks, mfttracks, ft0, false, col.posZ()); + fillCorrFt0aMft(sameFt0aMft, mfttracks, ft0, false, col.posZ(), multiplicity); } } // same event - void processFt0aFt0cSE(CollTable::iterator const& col, aod::FT0s const&, TrksTable const& tracks) + void processFt0aFt0cSE(CollTable::iterator const& col, aod::FT0s const&) { if (!isEventSelected(col)) { return; } if (col.has_foundFT0()) { - histos.fill(HIST("Ft0aFt0c/SE/hMult"), tracks.size()); + auto multiplicity = col.multiplicity(); + histos.fill(HIST("Ft0aFt0c/SE/hMult"), multiplicity); const auto& ft0 = col.foundFT0(); - if (tracks.size() < cfgMinMult || tracks.size() >= cfgMaxMult) { + if (multiplicity < cfgMinMult || multiplicity >= cfgMaxMult) { return; } - fillCorrFt0aFt0c(sameFt0aFt0c, tracks, ft0, ft0, false, col.posZ()); + fillCorrFt0aFt0c(sameFt0aFt0c, ft0, ft0, false, col.posZ(), multiplicity); } } // same event @@ -665,10 +681,11 @@ struct LongrangeCorrelation { auto slicedTriggerTracks = tracks.sliceBy(perColGlobal, col1.globalIndex()); fillYield(slicedTriggerTracks); const auto& ft0 = col2.foundFT0(); - if (slicedTriggerTracks.size() < cfgMinMult || slicedTriggerTracks.size() >= cfgMaxMult) { + auto multiplicity = col1.multiplicity(); + if (multiplicity < cfgMinMult || multiplicity >= cfgMaxMult) { continue; } - fillCorrFt0aGlobal(mixedFt0aGlobal, slicedTriggerTracks, ft0, true, col1.posZ()); + fillCorrFt0aGlobal(mixedFt0aGlobal, slicedTriggerTracks, ft0, true, col1.posZ(), multiplicity); } } } // mixed event @@ -693,10 +710,11 @@ struct LongrangeCorrelation { auto slicedTriggerTracks = tracks.sliceBy(perColGlobal, col1.globalIndex()); fillYield(slicedTriggerTracks); const auto& ft0 = col2.foundFT0(); - if (slicedTriggerTracks.size() < cfgMinMult || slicedTriggerTracks.size() >= cfgMaxMult) { + auto multiplicity = col1.multiplicity(); + if (multiplicity < cfgMinMult || multiplicity >= cfgMaxMult) { continue; } - fillCorrFt0cGlobal(mixedFt0cGlobal, slicedTriggerTracks, ft0, true, col1.posZ()); + fillCorrFt0cGlobal(mixedFt0cGlobal, slicedTriggerTracks, ft0, true, col1.posZ(), multiplicity); } } } // mixed event @@ -716,10 +734,11 @@ struct LongrangeCorrelation { if (!isEventSelected(col1) || !isEventSelected(col2)) { continue; } - if ((tracks1.size() < cfgMinMult || tracks1.size() >= cfgMaxMult)) { + auto multiplicity = col1.multiplicity(); + if (multiplicity < cfgMinMult || multiplicity >= cfgMaxMult) { continue; } - fillCorrMftGlobal(mixedMftGlobal, tracks1, tracks2, true, col1.posZ()); + fillCorrMftGlobal(mixedMftGlobal, tracks1, tracks2, true, col1.posZ(), multiplicity); } } // mixed event @@ -740,14 +759,14 @@ struct LongrangeCorrelation { continue; } if (col1.has_foundFT0() && col2.has_foundFT0()) { - auto slicedGlobalTracks = tracks.sliceBy(perColGlobal, col1.globalIndex()); auto slicedTriggerMftTracks = mfttracks.sliceBy(perColMft, col1.globalIndex()); fillYield(slicedTriggerMftTracks); const auto& ft0 = col2.foundFT0(); - if (slicedGlobalTracks.size() < cfgMinMult || slicedGlobalTracks.size() >= cfgMaxMult) { + auto multiplicity = col1.multiplicity(); + if (multiplicity < cfgMinMult || multiplicity >= cfgMaxMult) { continue; } - fillCorrFt0aMft(mixedFt0aMft, slicedGlobalTracks, slicedTriggerMftTracks, ft0, true, col1.posZ()); + fillCorrFt0aMft(mixedFt0aMft, slicedTriggerMftTracks, ft0, true, col1.posZ(), multiplicity); } } } // mixed event @@ -769,14 +788,14 @@ struct LongrangeCorrelation { continue; } if (col1.has_foundFT0() && col2.has_foundFT0()) { - auto slicedTriggerTracks = tracks.sliceBy(perColGlobal, col1.globalIndex()); - histos.fill(HIST("Ft0aFt0c/ME/hMult"), slicedTriggerTracks.size()); + auto multiplicity = col1.multiplicity(); + histos.fill(HIST("Ft0aFt0c/ME/hMult"), multiplicity); const auto& ft0a = col1.foundFT0(); const auto& ft0c = col2.foundFT0(); - if (slicedTriggerTracks.size() < cfgMinMult || slicedTriggerTracks.size() >= cfgMaxMult) { + if (multiplicity < cfgMinMult || multiplicity >= cfgMaxMult) { continue; } - fillCorrFt0aFt0c(mixedFt0aFt0c, slicedTriggerTracks, ft0a, ft0c, true, col1.posZ()); + fillCorrFt0aFt0c(mixedFt0aFt0c, ft0a, ft0c, true, col1.posZ(), multiplicity); } } } // mixed event @@ -794,7 +813,89 @@ struct LongrangeCorrelation { PROCESS_SWITCH(LongrangeCorrelation, processFt0aFt0cME, "mixed event FT0a vs FT0c", false); }; +struct MultiplicityClassifier { + Produces multvalue; + Configurable cfgEtaCut{"cfgEtaCut", 0.8f, "Eta range to consider"}; + Configurable dcaZ{"dcaZ", 0.2f, "Custom DCA Z cut (ignored if negative)"}; + Configurable cfgPtCutMin{"cfgPtCutMin", 0.2f, "minimum accepted track pT"}; + Configurable cfgPtCutMax{"cfgPtCutMax", 3.0f, "maximum accepted track pT"}; + HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; + + Filter fTrackSelectionITS = ncheckbit(aod::track::v001::detectorMap, (uint8_t)o2::aod::track::ITS) && + ncheckbit(aod::track::trackCutFlag, TrackSelectionIts); + Filter fTrackSelectionTPC = ifnode(ncheckbit(aod::track::v001::detectorMap, (uint8_t)o2::aod::track::TPC), + ncheckbit(aod::track::trackCutFlag, TrackSelectionTpc), true); + Filter fTrackSelectionDCA = ifnode(dcaZ.node() > 0.f, nabs(aod::track::dcaZ) <= dcaZ && ncheckbit(aod::track::trackCutFlag, TrackSelectionDcaxyOnly), + ncheckbit(aod::track::trackCutFlag, TrackSelectionDca)); + Filter fTracksEta = nabs(aod::track::eta) < cfgEtaCut; + Filter fTracksPt = (aod::track::pt > cfgPtCutMin) && (aod::track::pt < cfgPtCutMax); + + void init(InitContext const&) + { + int enabledFunctions = 0; + if (doprocessTracks) { + histos.add("htrackPt", "htrackPt", {HistType::kTH1F, {{10, 0., 10.}}}); + histos.add("htrackMult", "htrackMult", {HistType::kTH1F, {{500, 0., 500.}}}); + enabledFunctions++; + } + if (doprocessFT0C) { + histos.add("hCentFt0c", "hCentFt0c", {HistType::kTH1F, {{100, 0., 100.}}}); + enabledFunctions++; + } + if (doprocessFV0A) { + histos.add("hCentFv0a", "hCentFv0a", {HistType::kTH1F, {{100, 0., 100.}}}); + enabledFunctions++; + } + if (doprocessFT0M) { + histos.add("hCentFt0m", "hCentFt0m", {HistType::kTH1F, {{100, 0., 100.}}}); + enabledFunctions++; + } + if (enabledFunctions != 1) { + LOGP(fatal, "{} multiplicity classifier enabled but we need exactly 1.", enabledFunctions); + } + } + + void processTracks(aod::Collision const&, soa::Filtered> const& tracks) + { + multvalue(tracks.size()); + histos.fill(HIST("htrackMult"), tracks.size()); + for (auto const& iTrk : tracks) + histos.fill(HIST("htrackPt"), iTrk.pt()); + } + + void processFT0C(aod::CentFT0Cs const& centralities) + { + for (auto const& c : centralities) { + multvalue(c.centFT0C()); + histos.fill(HIST("hCentFt0c"), c.centFT0C()); + } + } + + void processFV0A(aod::CentFV0As const& centralities) + { + for (auto const& c : centralities) { + multvalue(c.centFV0A()); + histos.fill(HIST("hCentFv0a"), c.centFV0A()); + } + } + + void processFT0M(aod::CentFT0Ms const& centralities) + { + for (auto const& c : centralities) { + multvalue(c.centFT0M()); + histos.fill(HIST("hCentFt0m"), c.centFT0M()); + } + } + + PROCESS_SWITCH(MultiplicityClassifier, processTracks, "Select track count as multiplicity", false); + PROCESS_SWITCH(MultiplicityClassifier, processFT0C, "Select FT0C centrality as multiplicity", false); + PROCESS_SWITCH(MultiplicityClassifier, processFV0A, "Select FV0A centrality as multiplicity", false); + PROCESS_SWITCH(MultiplicityClassifier, processFT0M, "Select FT0M centrality as multiplicity", false); +}; + WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { - return WorkflowSpec{adaptAnalysisTask(cfgc)}; + return WorkflowSpec{ + adaptAnalysisTask(cfgc), + adaptAnalysisTask(cfgc)}; } From bc288807bda996e77103d92652c5a987d0000c49 Mon Sep 17 00:00:00 2001 From: abmodak <67369858+abmodak@users.noreply.github.com> Date: Mon, 6 Oct 2025 20:37:41 +0200 Subject: [PATCH 2/5] Add efficiency estimation --- .../Tasks/longrangeCorrelation.cxx | 85 ++++++++++++++++++- 1 file changed, 84 insertions(+), 1 deletion(-) diff --git a/PWGCF/TwoParticleCorrelations/Tasks/longrangeCorrelation.cxx b/PWGCF/TwoParticleCorrelations/Tasks/longrangeCorrelation.cxx index bc5949ba294..b490b511768 100644 --- a/PWGCF/TwoParticleCorrelations/Tasks/longrangeCorrelation.cxx +++ b/PWGCF/TwoParticleCorrelations/Tasks/longrangeCorrelation.cxx @@ -44,11 +44,13 @@ #include "Framework/HistogramRegistry.h" #include "Framework/StepTHn.h" #include "Framework/runDataProcessing.h" +#include "Framework/O2DatabasePDGPlugin.h" #include "ReconstructionDataFormats/Track.h" #include #include #include +#include #include #include @@ -91,6 +93,7 @@ enum KindOfCorrType { static constexpr std::string_view kCorrType[] = {"Ft0aGlobal/", "Ft0cGlobal/", "MftGlobal/", "Ft0aMft/", "Ft0aFt0c/"}; static constexpr std::string_view kEvntType[] = {"SE/", "ME/"}; auto static constexpr kMinFt0cCell = 96; +auto static constexpr kMinCharge = 3.f; AxisSpec axisEvent{10, 0.5, 9.5, "#Event", "EventAxis"}; namespace o2::aod @@ -113,7 +116,9 @@ struct LongrangeCorrelation { SliceCache cache; Service ccdb; + Service pdg; o2::ccdb::CcdbApi ccdbApi; + o2::ft0::Geometry ft0Det; std::vector* offsetFT0; HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; Configurable cfgVtxCut{"cfgVtxCut", 10.0f, "Vertex Z range to consider"}; @@ -153,9 +158,11 @@ struct LongrangeCorrelation { using CollTable = soa::Join; using TrksTable = soa::Filtered>; using MftTrkTable = soa::Filtered; + using CollTableMC = soa::SmallGroups>; + using TrksTableMC = soa::Filtered>; Preslice perColGlobal = aod::track::collisionId; + Preslice perColMC = aod::track::collisionId; Preslice perColMft = aod::fwdtrack::collisionId; - o2::ft0::Geometry ft0Det; OutputObj sameFt0aGlobal{Form("sameEventFt0aGlobal_%i_%i", static_cast(cfgMinMult), static_cast(cfgMaxMult))}; OutputObj mixedFt0aGlobal{Form("mixedEventFt0aGlobal_%i_%i", static_cast(cfgMinMult), static_cast(cfgMaxMult))}; @@ -256,6 +263,13 @@ struct LongrangeCorrelation { sameFt0aFt0c.setObject(new CorrelationContainer(Form("sameEventFt0aFt0c_%i_%i", static_cast(cfgMinMult), static_cast(cfgMaxMult)), Form("sameEventFt0aFt0c_%i_%i", static_cast(cfgMinMult), static_cast(cfgMaxMult)), corrAxis, effAxis, userAxis)); mixedFt0aFt0c.setObject(new CorrelationContainer(Form("mixedEventFt0aFt0c_%i_%i", static_cast(cfgMinMult), static_cast(cfgMaxMult)), Form("mixedEventFt0aFt0c_%i_%i", static_cast(cfgMinMult), static_cast(cfgMaxMult)), corrAxis, effAxis, userAxis)); } + + if (doprocessEff) { + histos.add("hmcgendndptPrimary", "hmcgendndptPrimary", kTHnSparseD, {axisEtaTrig, axisPtTrigger, axisMultiplicity, axisVtxZ}, false); + histos.add("hmcrecdndptRecoPrimary", "hmcrecdndptRecoPrimary", kTHnSparseD, {axisEtaTrig, axisPtTrigger, axisMultiplicity, axisVtxZ}, false); + histos.add("hmcrecdndptRecoAll", "hmcrecdndptRecoAll", kTHnSparseD, {axisEtaTrig, axisPtTrigger, axisMultiplicity, axisVtxZ}, false); + histos.add("hmcrecdndptFake", "hmcrecdndptFake", kTHnSparseD, {axisEtaTrig, axisPtTrigger, axisMultiplicity, axisVtxZ}, false); + } } Filter fTrackSelectionITS = ncheckbit(aod::track::v001::detectorMap, (uint8_t)o2::aod::track::ITS) && @@ -800,6 +814,74 @@ struct LongrangeCorrelation { } } // mixed event + template + bool isGenTrackSelected(CheckGenTrack const& track) + { + if (!track.isPhysicalPrimary()) { + return false; + } + if (!track.producedByGenerator()) { + return false; + } + auto pdgTrack = pdg->GetParticle(track.pdgCode()); + if (pdgTrack == nullptr) { + return false; + } + if (std::abs(pdgTrack->Charge()) < kMinCharge) { + return false; + } + if (std::abs(track.eta()) >= cfgEtaCut) { + return false; + } + return true; + } + + void processEff(aod::McCollisions::iterator const& mcCollision, CollTableMC const& RecCols, aod::McParticles const& GenParticles, TrksTableMC const& RecTracks) + { + if (std::abs(mcCollision.posZ()) >= cfgVtxCut) { + return; + } + + auto multiplicity = -999.; + auto numcontributors = -999; + for (const auto& RecCol : RecCols) { + if (!isEventSelected(RecCol)) { + continue; + } + if (RecCol.numContrib() <= numcontributors) { + continue; + } else { + numcontributors = RecCol.numContrib(); + } + multiplicity = RecCol.multiplicity(); + } + + for (const auto& particle : GenParticles) { + if (!isGenTrackSelected(particle)) { + continue; + } + histos.fill(HIST("hmcgendndptPrimary"), particle.eta(), particle.pt(), multiplicity, mcCollision.posZ()); + } // track (mcgen) loop + + for (const auto& RecCol : RecCols) { + if (!isEventSelected(RecCol)) { + continue; + } + auto recTracksPart = RecTracks.sliceBy(perColMC, RecCol.globalIndex()); + for (const auto& Rectrack : recTracksPart) { + if (Rectrack.has_mcParticle()) { + auto mcpart = Rectrack.mcParticle(); + histos.fill(HIST("hmcrecdndptRecoAll"), mcpart.eta(), mcpart.pt(), multiplicity, mcCollision.posZ()); + if (mcpart.isPhysicalPrimary()) { + histos.fill(HIST("hmcrecdndptRecoPrimary"), mcpart.eta(), mcpart.pt(), multiplicity, mcCollision.posZ()); + } + } else { + histos.fill(HIST("hmcrecdndptFake"), Rectrack.eta(), Rectrack.pt(), multiplicity, mcCollision.posZ()); + } + } // track (mcrec) loop + } // rec collision + } + PROCESS_SWITCH(LongrangeCorrelation, processEventStat, "event stat", false); PROCESS_SWITCH(LongrangeCorrelation, processFt0aGlobalSE, "same event FT0a vs global", false); PROCESS_SWITCH(LongrangeCorrelation, processFt0aGlobalME, "mixed event FT0a vs global", false); @@ -811,6 +893,7 @@ struct LongrangeCorrelation { PROCESS_SWITCH(LongrangeCorrelation, processFt0aMftME, "mixed event FT0a vs MFT", false); PROCESS_SWITCH(LongrangeCorrelation, processFt0aFt0cSE, "same event FT0a vs FT0c", false); PROCESS_SWITCH(LongrangeCorrelation, processFt0aFt0cME, "mixed event FT0a vs FT0c", false); + PROCESS_SWITCH(LongrangeCorrelation, processEff, "Estimate efficiency", false); }; struct MultiplicityClassifier { From 9a4f3c26a7e06b6effe7a0d55a8c34bc88313f5b Mon Sep 17 00:00:00 2001 From: abmodak <67369858+abmodak@users.noreply.github.com> Date: Wed, 8 Oct 2025 17:38:17 +0200 Subject: [PATCH 3/5] Load and apply efficiency correction --- .../Tasks/longrangeCorrelation.cxx | 196 +++++++++++++++--- 1 file changed, 168 insertions(+), 28 deletions(-) diff --git a/PWGCF/TwoParticleCorrelations/Tasks/longrangeCorrelation.cxx b/PWGCF/TwoParticleCorrelations/Tasks/longrangeCorrelation.cxx index b490b511768..430e70cf74f 100644 --- a/PWGCF/TwoParticleCorrelations/Tasks/longrangeCorrelation.cxx +++ b/PWGCF/TwoParticleCorrelations/Tasks/longrangeCorrelation.cxx @@ -42,9 +42,9 @@ #include "Framework/AnalysisDataModel.h" #include "Framework/AnalysisTask.h" #include "Framework/HistogramRegistry.h" +#include "Framework/O2DatabasePDGPlugin.h" #include "Framework/StepTHn.h" #include "Framework/runDataProcessing.h" -#include "Framework/O2DatabasePDGPlugin.h" #include "ReconstructionDataFormats/Track.h" #include @@ -137,6 +137,9 @@ struct LongrangeCorrelation { Configurable isApplySameBunchPileup{"isApplySameBunchPileup", false, "Enable SameBunchPileup cut"}; Configurable isApplyGoodZvtxFT0vsPV{"isApplyGoodZvtxFT0vsPV", false, "Enable GoodZvtxFT0vsPV cut"}; Configurable isReadoutCenter{"isReadoutCenter", false, "Enable Readout Center"}; + Configurable isUseEffCorr{"isUseEffCorr", false, "Enable efficiency correction"}; + Configurable cfgEffccdbPath{"cfgEffccdbPath", "/alice/data/CCDB/Users/a/abmodak/OO/Efficiency", "Browse track eff object from CCDB"}; + Configurable cfgMultccdbPath{"cfgMultccdbPath", "/alice/data/CCDB/Users/a/abmodak/OO/Multiplicity", "Browse mult efficiency object from CCDB"}; ConfigurableAxis axisDeltaPhi{"axisDeltaPhi", {72, -PIHalf, PIHalf * 3}, "delta phi axis for histograms"}; ConfigurableAxis axisDeltaEta{"axisDeltaEta", {40, -6, -2}, "delta eta axis for histograms"}; ConfigurableAxis axisPtTrigger{"axisPtTrigger", {VARIABLE_WIDTH, 0.5, 1.0, 1.5, 2.0, 3.0, 4.0, 6.0, 10.0}, "pt trigger axis for histograms"}; @@ -175,6 +178,12 @@ struct LongrangeCorrelation { OutputObj sameFt0aFt0c{Form("sameEventFt0aFt0c_%i_%i", static_cast(cfgMinMult), static_cast(cfgMaxMult))}; OutputObj mixedFt0aFt0c{Form("mixedEventFt0aFt0c_%i_%i", static_cast(cfgMinMult), static_cast(cfgMaxMult))}; + // corrections + TH3D* hTrkEff = nullptr; + TH1D* hMultEff = nullptr; + bool fLoadTrkEffCorr = false; + bool fLoadMultEffCorr = false; + template void addHistos() { @@ -349,15 +358,82 @@ struct LongrangeCorrelation { } } + void loadEffCorrection(uint64_t timestamp) + { + if (fLoadTrkEffCorr) { + return; + } + if (cfgEffccdbPath.value.empty() == false) { + hTrkEff = ccdb->getForTimeStamp(cfgEffccdbPath, timestamp); + if (hTrkEff == nullptr) { + LOGF(fatal, "Could not load efficiency histogram for trigger particles from %s", cfgEffccdbPath.value.c_str()); + } + LOGF(info, "Loaded efficiency histogram from %s (%p)", cfgEffccdbPath.value.c_str(), (void*)hTrkEff); + } + fLoadTrkEffCorr = true; + } + + void loadMultCorrection(uint64_t timestamp) + { + if (fLoadMultEffCorr) { + return; + } + if (cfgMultccdbPath.value.empty() == false) { + hMultEff = ccdb->getForTimeStamp(cfgMultccdbPath, timestamp); + if (hMultEff == nullptr) { + LOGF(fatal, "Could not load efficiency histogram for multiplicity from %s", cfgMultccdbPath.value.c_str()); + } + LOGF(info, "Loaded efficiency histogram from %s (%p)", cfgMultccdbPath.value.c_str(), (void*)hMultEff); + } + fLoadMultEffCorr = true; + } + + float getTrkEffCorr(float eta, float pt, float posZ) + { + float eff = 1.; + if (hTrkEff) { + int etaBin = hTrkEff->GetXaxis()->FindBin(eta); + int ptBin = hTrkEff->GetYaxis()->FindBin(pt); + int zBin = hTrkEff->GetZaxis()->FindBin(posZ); + eff = hTrkEff->GetBinContent(etaBin, ptBin, zBin); + } else { + eff = 1.0; + } + if (eff < 0.001) + eff = 1.0; + + return eff; + } + + float getMultEffCorr(float mult) + { + float eff = 1.; + if (hMultEff) { + int multBin = hMultEff->GetXaxis()->FindBin(mult); + eff = hMultEff->GetBinContent(multBin); + } else { + eff = 1.0; + } + if (eff < 0.001) + eff = 1.0; + + return eff; + } + template void fillCorrFt0aGlobal(TTarget target, TTriggers const& triggers, TFT0s const& ft0, bool mixing, float vz, float multiplicity) { int fSampleIndex = gRandom->Uniform(0, cfgSampleSize); if (!mixing) histos.fill(HIST("Ft0aGlobal/SE/hMult_used"), multiplicity); + for (auto const& triggerTrack : triggers) { + float trkeffw = 1.0f; + if (isUseEffCorr) + trkeffw = getTrkEffCorr(triggerTrack.eta(), triggerTrack.pt(), vz); + if (!mixing) - histos.fill(HIST("Ft0aGlobal/SE/Trig_hist"), fSampleIndex, vz, triggerTrack.pt()); + histos.fill(HIST("Ft0aGlobal/SE/Trig_hist"), fSampleIndex, vz, triggerTrack.pt(), trkeffw); for (std::size_t iCh = 0; iCh < ft0.channelA().size(); iCh++) { auto chanelid = ft0.channelA()[iCh]; @@ -383,10 +459,10 @@ struct LongrangeCorrelation { float deltaPhi = RecoDecay::constrainAngle(triggerTrack.phi() - phi, -PIHalf); float deltaEta = triggerTrack.eta() - eta; if (mixing) - histos.fill(HIST("Ft0aGlobal/ME/deltaEta_deltaPhi"), deltaPhi, deltaEta); + histos.fill(HIST("Ft0aGlobal/ME/deltaEta_deltaPhi"), deltaPhi, deltaEta, trkeffw); else - histos.fill(HIST("Ft0aGlobal/SE/deltaEta_deltaPhi"), deltaPhi, deltaEta); - target->getPairHist()->Fill(step, fSampleIndex, vz, triggerTrack.pt(), triggerTrack.pt(), deltaPhi, deltaEta); + histos.fill(HIST("Ft0aGlobal/SE/deltaEta_deltaPhi"), deltaPhi, deltaEta, trkeffw); + target->getPairHist()->Fill(step, fSampleIndex, vz, triggerTrack.pt(), triggerTrack.pt(), deltaPhi, deltaEta, trkeffw); } // associated ft0 tracks } // trigger tracks } // fillCorrFt0aGlobal @@ -398,8 +474,12 @@ struct LongrangeCorrelation { if (!mixing) histos.fill(HIST("Ft0cGlobal/SE/hMult_used"), multiplicity); for (auto const& triggerTrack : triggers) { + float trkeffw = 1.0f; + if (isUseEffCorr) + trkeffw = getTrkEffCorr(triggerTrack.eta(), triggerTrack.pt(), vz); + if (!mixing) - histos.fill(HIST("Ft0cGlobal/SE/Trig_hist"), fSampleIndex, vz, triggerTrack.pt()); + histos.fill(HIST("Ft0cGlobal/SE/Trig_hist"), fSampleIndex, vz, triggerTrack.pt(), trkeffw); for (std::size_t iCh = 0; iCh < ft0.channelC().size(); iCh++) { auto chanelid = ft0.channelC()[iCh] + 96; @@ -425,10 +505,10 @@ struct LongrangeCorrelation { float deltaPhi = RecoDecay::constrainAngle(triggerTrack.phi() - phi, -PIHalf); float deltaEta = triggerTrack.eta() - eta; if (mixing) - histos.fill(HIST("Ft0cGlobal/ME/deltaEta_deltaPhi"), deltaPhi, deltaEta); + histos.fill(HIST("Ft0cGlobal/ME/deltaEta_deltaPhi"), deltaPhi, deltaEta, trkeffw); else - histos.fill(HIST("Ft0cGlobal/SE/deltaEta_deltaPhi"), deltaPhi, deltaEta); - target->getPairHist()->Fill(step, fSampleIndex, vz, triggerTrack.pt(), triggerTrack.pt(), deltaPhi, deltaEta); + histos.fill(HIST("Ft0cGlobal/SE/deltaEta_deltaPhi"), deltaPhi, deltaEta, trkeffw); + target->getPairHist()->Fill(step, fSampleIndex, vz, triggerTrack.pt(), triggerTrack.pt(), deltaPhi, deltaEta, trkeffw); } // associated ft0 tracks } // trigger tracks } // fillCorrFt0cGlobal @@ -440,8 +520,12 @@ struct LongrangeCorrelation { if (!mixing) histos.fill(HIST("MftGlobal/SE/hMult_used"), multiplicity); for (auto const& triggerTrack : triggers) { + float trkeffw = 1.0f; + if (isUseEffCorr) + trkeffw = getTrkEffCorr(triggerTrack.eta(), triggerTrack.pt(), vz); + if (!mixing) - histos.fill(HIST("MftGlobal/SE/Trig_hist"), fSampleIndex, vz, triggerTrack.pt()); + histos.fill(HIST("MftGlobal/SE/Trig_hist"), fSampleIndex, vz, triggerTrack.pt(), trkeffw); for (auto const& assoTrack : mft) { if (!isMftTrackSelected(assoTrack)) { @@ -461,10 +545,10 @@ struct LongrangeCorrelation { float deltaPhi = RecoDecay::constrainAngle(triggerTrack.phi() - phi, -PIHalf); float deltaEta = triggerTrack.eta() - assoTrack.eta(); if (mixing) - histos.fill(HIST("MftGlobal/ME/deltaEta_deltaPhi"), deltaPhi, deltaEta); + histos.fill(HIST("MftGlobal/ME/deltaEta_deltaPhi"), deltaPhi, deltaEta, trkeffw); else - histos.fill(HIST("MftGlobal/SE/deltaEta_deltaPhi"), deltaPhi, deltaEta); - target->getPairHist()->Fill(step, fSampleIndex, vz, triggerTrack.pt(), assoTrack.pt(), deltaPhi, deltaEta); + histos.fill(HIST("MftGlobal/SE/deltaEta_deltaPhi"), deltaPhi, deltaEta, trkeffw); + target->getPairHist()->Fill(step, fSampleIndex, vz, triggerTrack.pt(), assoTrack.pt(), deltaPhi, deltaEta, trkeffw); } // associated mft tracks } // trigger tracks } // fillCorrMftGlobal @@ -598,15 +682,21 @@ struct LongrangeCorrelation { histos.fill(HIST("QA/VtxZHist"), col.posZ()); } - void processFt0aGlobalSE(CollTable::iterator const& col, aod::FT0s const&, TrksTable const& tracks) + void processFt0aGlobalSE(CollTable::iterator const& col, aod::FT0s const&, TrksTable const& tracks, aod::BCsWithTimestamps const&) { if (!isEventSelected(col)) { return; } if (col.has_foundFT0()) { + auto bc = col.bc_as(); + loadEffCorrection(bc.timestamp()); + loadMultCorrection(bc.timestamp()); fillYield(tracks); const auto& ft0 = col.foundFT0(); auto multiplicity = col.multiplicity(); + float multw = getMultEffCorr(multiplicity); + if (isUseEffCorr) + multiplicity = multiplicity * multw; if (multiplicity < cfgMinMult || multiplicity >= cfgMaxMult) { return; } @@ -614,15 +704,21 @@ struct LongrangeCorrelation { } } // same event - void processFt0cGlobalSE(CollTable::iterator const& col, aod::FT0s const&, TrksTable const& tracks) + void processFt0cGlobalSE(CollTable::iterator const& col, aod::FT0s const&, TrksTable const& tracks, aod::BCsWithTimestamps const&) { if (!isEventSelected(col)) { return; } if (col.has_foundFT0()) { + auto bc = col.bc_as(); + loadEffCorrection(bc.timestamp()); + loadMultCorrection(bc.timestamp()); fillYield(tracks); const auto& ft0 = col.foundFT0(); auto multiplicity = col.multiplicity(); + float multw = getMultEffCorr(multiplicity); + if (isUseEffCorr) + multiplicity = multiplicity * multw; if (multiplicity < cfgMinMult || multiplicity >= cfgMaxMult) { return; } @@ -630,28 +726,39 @@ struct LongrangeCorrelation { } } // same event - void processMftGlobalSE(CollTable::iterator const& col, MftTrkTable const& mfttracks, TrksTable const& tracks) + void processMftGlobalSE(CollTable::iterator const& col, MftTrkTable const& mfttracks, TrksTable const& tracks, aod::BCsWithTimestamps const&) { if (!isEventSelected(col)) { return; } + auto bc = col.bc_as(); + loadEffCorrection(bc.timestamp()); + loadMultCorrection(bc.timestamp()); fillYield(tracks); auto multiplicity = col.multiplicity(); + float multw = getMultEffCorr(multiplicity); + if (isUseEffCorr) + multiplicity = multiplicity * multw; if (multiplicity < cfgMinMult || multiplicity >= cfgMaxMult) { return; } fillCorrMftGlobal(sameMftGlobal, tracks, mfttracks, false, col.posZ(), multiplicity); } // same event - void processFt0aMftSE(CollTable::iterator const& col, aod::FT0s const&, MftTrkTable const& mfttracks) + void processFt0aMftSE(CollTable::iterator const& col, aod::FT0s const&, MftTrkTable const& mfttracks, aod::BCsWithTimestamps const&) { if (!isEventSelected(col)) { return; } if (col.has_foundFT0()) { fillYield(mfttracks); + auto bc = col.bc_as(); + loadMultCorrection(bc.timestamp()); const auto& ft0 = col.foundFT0(); auto multiplicity = col.multiplicity(); + float multw = getMultEffCorr(multiplicity); + if (isUseEffCorr) + multiplicity = multiplicity * multw; if (multiplicity < cfgMinMult || multiplicity >= cfgMaxMult) { return; } @@ -659,15 +766,20 @@ struct LongrangeCorrelation { } } // same event - void processFt0aFt0cSE(CollTable::iterator const& col, aod::FT0s const&) + void processFt0aFt0cSE(CollTable::iterator const& col, aod::FT0s const&, aod::BCsWithTimestamps const&) { if (!isEventSelected(col)) { return; } if (col.has_foundFT0()) { + auto bc = col.bc_as(); + loadMultCorrection(bc.timestamp()); auto multiplicity = col.multiplicity(); histos.fill(HIST("Ft0aFt0c/SE/hMult"), multiplicity); const auto& ft0 = col.foundFT0(); + float multw = getMultEffCorr(multiplicity); + if (isUseEffCorr) + multiplicity = multiplicity * multw; if (multiplicity < cfgMinMult || multiplicity >= cfgMaxMult) { return; } @@ -675,7 +787,7 @@ struct LongrangeCorrelation { } } // same event - void processFt0aGlobalME(CollTable const& col, aod::FT0s const&, TrksTable const& tracks) + void processFt0aGlobalME(CollTable const& col, aod::FT0s const&, TrksTable const& tracks, aod::BCsWithTimestamps const&) { auto getTracksSize = [&tracks, this](CollTable::iterator const& collision) { auto associatedTracks = tracks.sliceByCached(o2::aod::track::collisionId, collision.globalIndex(), this->cache); @@ -692,10 +804,16 @@ struct LongrangeCorrelation { continue; } if (col1.has_foundFT0() && col2.has_foundFT0()) { + auto bc = col1.bc_as(); + loadEffCorrection(bc.timestamp()); + loadMultCorrection(bc.timestamp()); auto slicedTriggerTracks = tracks.sliceBy(perColGlobal, col1.globalIndex()); fillYield(slicedTriggerTracks); const auto& ft0 = col2.foundFT0(); auto multiplicity = col1.multiplicity(); + float multw = getMultEffCorr(multiplicity); + if (isUseEffCorr) + multiplicity = multiplicity * multw; if (multiplicity < cfgMinMult || multiplicity >= cfgMaxMult) { continue; } @@ -704,7 +822,7 @@ struct LongrangeCorrelation { } } // mixed event - void processFt0cGlobalME(CollTable const& col, aod::FT0s const&, TrksTable const& tracks) + void processFt0cGlobalME(CollTable const& col, aod::FT0s const&, TrksTable const& tracks, aod::BCsWithTimestamps const&) { auto getTracksSize = [&tracks, this](CollTable::iterator const& collision) { auto associatedTracks = tracks.sliceByCached(o2::aod::track::collisionId, collision.globalIndex(), this->cache); @@ -721,10 +839,16 @@ struct LongrangeCorrelation { continue; } if (col1.has_foundFT0() && col2.has_foundFT0()) { + auto bc = col1.bc_as(); + loadEffCorrection(bc.timestamp()); + loadMultCorrection(bc.timestamp()); auto slicedTriggerTracks = tracks.sliceBy(perColGlobal, col1.globalIndex()); fillYield(slicedTriggerTracks); const auto& ft0 = col2.foundFT0(); auto multiplicity = col1.multiplicity(); + float multw = getMultEffCorr(multiplicity); + if (isUseEffCorr) + multiplicity = multiplicity * multw; if (multiplicity < cfgMinMult || multiplicity >= cfgMaxMult) { continue; } @@ -733,7 +857,7 @@ struct LongrangeCorrelation { } } // mixed event - void processMftGlobalME(CollTable const& col, MftTrkTable const& mfttracks, TrksTable const& tracks) + void processMftGlobalME(CollTable const& col, MftTrkTable const& mfttracks, TrksTable const& tracks, aod::BCsWithTimestamps const&) { auto getTracksSize = [&tracks, this](CollTable::iterator const& collision) { auto associatedTracks = tracks.sliceByCached(o2::aod::track::collisionId, collision.globalIndex(), this->cache); @@ -748,7 +872,13 @@ struct LongrangeCorrelation { if (!isEventSelected(col1) || !isEventSelected(col2)) { continue; } + auto bc = col1.bc_as(); + loadEffCorrection(bc.timestamp()); + loadMultCorrection(bc.timestamp()); auto multiplicity = col1.multiplicity(); + float multw = getMultEffCorr(multiplicity); + if (isUseEffCorr) + multiplicity = multiplicity * multw; if (multiplicity < cfgMinMult || multiplicity >= cfgMaxMult) { continue; } @@ -756,7 +886,7 @@ struct LongrangeCorrelation { } } // mixed event - void processFt0aMftME(CollTable const& col, aod::FT0s const&, TrksTable const& tracks, MftTrkTable const& mfttracks) + void processFt0aMftME(CollTable const& col, aod::FT0s const&, TrksTable const& tracks, MftTrkTable const& mfttracks, aod::BCsWithTimestamps const&) { auto getTracksSize = [&tracks, this](CollTable::iterator const& collision) { auto associatedTracks = tracks.sliceByCached(o2::aod::track::collisionId, collision.globalIndex(), this->cache); @@ -773,10 +903,15 @@ struct LongrangeCorrelation { continue; } if (col1.has_foundFT0() && col2.has_foundFT0()) { + auto bc = col1.bc_as(); + loadMultCorrection(bc.timestamp()); auto slicedTriggerMftTracks = mfttracks.sliceBy(perColMft, col1.globalIndex()); fillYield(slicedTriggerMftTracks); const auto& ft0 = col2.foundFT0(); auto multiplicity = col1.multiplicity(); + float multw = getMultEffCorr(multiplicity); + if (isUseEffCorr) + multiplicity = multiplicity * multw; if (multiplicity < cfgMinMult || multiplicity >= cfgMaxMult) { continue; } @@ -785,7 +920,7 @@ struct LongrangeCorrelation { } } // mixed event - void processFt0aFt0cME(CollTable const& col, aod::FT0s const&, TrksTable const& tracks) + void processFt0aFt0cME(CollTable const& col, aod::FT0s const&, TrksTable const& tracks, aod::BCsWithTimestamps const&) { auto getTracksSize = [&tracks, this](CollTable::iterator const& collision) { auto associatedTracks = tracks.sliceByCached(o2::aod::track::collisionId, collision.globalIndex(), this->cache); @@ -802,10 +937,15 @@ struct LongrangeCorrelation { continue; } if (col1.has_foundFT0() && col2.has_foundFT0()) { + auto bc = col1.bc_as(); + loadMultCorrection(bc.timestamp()); auto multiplicity = col1.multiplicity(); histos.fill(HIST("Ft0aFt0c/ME/hMult"), multiplicity); const auto& ft0a = col1.foundFT0(); const auto& ft0c = col2.foundFT0(); + float multw = getMultEffCorr(multiplicity); + if (isUseEffCorr) + multiplicity = multiplicity * multw; if (multiplicity < cfgMinMult || multiplicity >= cfgMaxMult) { continue; } @@ -835,7 +975,7 @@ struct LongrangeCorrelation { } return true; } - + void processEff(aod::McCollisions::iterator const& mcCollision, CollTableMC const& RecCols, aod::McParticles const& GenParticles, TrksTableMC const& RecTracks) { if (std::abs(mcCollision.posZ()) >= cfgVtxCut) { @@ -870,12 +1010,12 @@ struct LongrangeCorrelation { auto recTracksPart = RecTracks.sliceBy(perColMC, RecCol.globalIndex()); for (const auto& Rectrack : recTracksPart) { if (Rectrack.has_mcParticle()) { - auto mcpart = Rectrack.mcParticle(); - histos.fill(HIST("hmcrecdndptRecoAll"), mcpart.eta(), mcpart.pt(), multiplicity, mcCollision.posZ()); + auto mcpart = Rectrack.mcParticle(); + histos.fill(HIST("hmcrecdndptRecoAll"), mcpart.eta(), mcpart.pt(), multiplicity, mcCollision.posZ()); if (mcpart.isPhysicalPrimary()) { histos.fill(HIST("hmcrecdndptRecoPrimary"), mcpart.eta(), mcpart.pt(), multiplicity, mcCollision.posZ()); - } - } else { + } + } else { histos.fill(HIST("hmcrecdndptFake"), Rectrack.eta(), Rectrack.pt(), multiplicity, mcCollision.posZ()); } } // track (mcrec) loop From 004e3c004cc599996a7901962ec4014692a6869f Mon Sep 17 00:00:00 2001 From: abmodak <67369858+abmodak@users.noreply.github.com> Date: Wed, 8 Oct 2025 19:27:26 +0200 Subject: [PATCH 4/5] Add PID --- .../Tasks/longrangeCorrelation.cxx | 115 ++++++++++++++++-- 1 file changed, 103 insertions(+), 12 deletions(-) diff --git a/PWGCF/TwoParticleCorrelations/Tasks/longrangeCorrelation.cxx b/PWGCF/TwoParticleCorrelations/Tasks/longrangeCorrelation.cxx index 430e70cf74f..64005e7d6c6 100644 --- a/PWGCF/TwoParticleCorrelations/Tasks/longrangeCorrelation.cxx +++ b/PWGCF/TwoParticleCorrelations/Tasks/longrangeCorrelation.cxx @@ -29,6 +29,7 @@ #include "Common/DataModel/FT0Corrected.h" #include "Common/DataModel/Multiplicity.h" #include "Common/DataModel/PIDResponse.h" +#include "Common/DataModel/PIDResponseITS.h" #include "Common/DataModel/TrackSelectionTables.h" #include "CCDB/BasicCCDBManager.h" @@ -45,6 +46,7 @@ #include "Framework/O2DatabasePDGPlugin.h" #include "Framework/StepTHn.h" #include "Framework/runDataProcessing.h" +#include "ReconstructionDataFormats/PID.h" #include "ReconstructionDataFormats/Track.h" #include @@ -90,6 +92,12 @@ enum KindOfCorrType { kFT0AFT0C }; +enum KindOfParticles { + PIONS, + KAONS, + PROTONS +}; + static constexpr std::string_view kCorrType[] = {"Ft0aGlobal/", "Ft0cGlobal/", "MftGlobal/", "Ft0aMft/", "Ft0aFt0c/"}; static constexpr std::string_view kEvntType[] = {"SE/", "ME/"}; auto static constexpr kMinFt0cCell = 96; @@ -138,6 +146,12 @@ struct LongrangeCorrelation { Configurable isApplyGoodZvtxFT0vsPV{"isApplyGoodZvtxFT0vsPV", false, "Enable GoodZvtxFT0vsPV cut"}; Configurable isReadoutCenter{"isReadoutCenter", false, "Enable Readout Center"}; Configurable isUseEffCorr{"isUseEffCorr", false, "Enable efficiency correction"}; + Configurable isUseItsPid{"isUseItsPid", false, "Use ITS PID for particle identification"}; + Configurable cfgTofPidPtCut{"cfgTofPidPtCut", 0.3f, "Minimum pt to use TOF N-sigma"}; + Configurable cfgTrackPid{"cfgTrackPid", 0, "1 = pion, 2 = kaon, 3 = proton, 0 for no PID"}; + Configurable> itsNsigmaPidCut{"itsNsigmaPidCut", std::vector{3, 2.5, 2, -3, -2.5, -2}, "ITS n-sigma cut for pions_posNsigma, kaons_posNsigma, protons_posNsigma, pions_negNsigma, kaons_negNsigma, protons_negNsigma"}; + Configurable> tpcNsigmaPidCut{"tpcNsigmaPidCut", std::vector{1.5, 1.5, 1.5, -1.5, -1.5, -1.5}, "TPC n-sigma cut for pions_posNsigma, kaons_posNsigma, protons_posNsigma, pions_negNsigma, kaons_negNsigma, protons_negNsigma"}; + Configurable> tofNsigmaPidCut{"tofNsigmaPidCut", std::vector{1.5, 1.5, 1.5, -1.5, -1.5, -1.5}, "TOF n-sigma cut for pions_posNsigma, kaons_posNsigma, protons_posNsigma, pions_negNsigma, kaons_negNsigma, protons_negNsigma"}; Configurable cfgEffccdbPath{"cfgEffccdbPath", "/alice/data/CCDB/Users/a/abmodak/OO/Efficiency", "Browse track eff object from CCDB"}; Configurable cfgMultccdbPath{"cfgMultccdbPath", "/alice/data/CCDB/Users/a/abmodak/OO/Multiplicity", "Browse mult efficiency object from CCDB"}; ConfigurableAxis axisDeltaPhi{"axisDeltaPhi", {72, -PIHalf, PIHalf * 3}, "delta phi axis for histograms"}; @@ -159,7 +173,7 @@ struct LongrangeCorrelation { ConfigurableAxis channelFt0aAxis{"channelFt0aAxis", {96, 0.0, 96.0}, "FT0A channel"}; using CollTable = soa::Join; - using TrksTable = soa::Filtered>; + using TrksTable = soa::Filtered>; using MftTrkTable = soa::Filtered; using CollTableMC = soa::SmallGroups>; using TrksTableMC = soa::Filtered>; @@ -184,6 +198,11 @@ struct LongrangeCorrelation { bool fLoadTrkEffCorr = false; bool fLoadMultEffCorr = false; + std::vector tofNsigmaCut; + std::vector itsNsigmaCut; + std::vector tpcNsigmaCut; + o2::aod::ITSResponse itsResponse; + template void addHistos() { @@ -225,6 +244,10 @@ struct LongrangeCorrelation { std::vector userAxis; + tofNsigmaCut = tofNsigmaPidCut; + itsNsigmaCut = itsNsigmaPidCut; + tpcNsigmaCut = tpcNsigmaPidCut; + if (doprocessEventStat) { histos.add("QA/EventHist", "events", kTH1F, {axisEvent}, false); histos.add("QA/VtxZHist", "v_{z} (cm)", kTH1F, {axisVtxZ}, false); @@ -343,14 +366,26 @@ struct LongrangeCorrelation { } template - void fillYield(TTracks tracks) + void fillYieldTpc(TTracks tracks) + { + histos.fill(HIST(kCorrType[corrType]) + HIST(kEvntType[evntType]) + HIST("hMult"), tracks.size()); + for (auto const& iTrk : tracks) { + if (cfgTrackPid && getTrackPID(iTrk) != cfgTrackPid) + continue; // if PID is selected, check if the track has the right PID + histos.fill(HIST(kCorrType[corrType]) + HIST(kEvntType[evntType]) + HIST("Trig_etavsphi"), iTrk.phi(), iTrk.eta()); + histos.fill(HIST(kCorrType[corrType]) + HIST(kEvntType[evntType]) + HIST("Trig_eta"), iTrk.eta()); + histos.fill(HIST(kCorrType[corrType]) + HIST(kEvntType[evntType]) + HIST("Trig_phi"), iTrk.phi()); + histos.fill(HIST(kCorrType[corrType]) + HIST(kEvntType[evntType]) + HIST("Trig_pt"), iTrk.pt()); + } + } + + template + void fillYieldMft(TTracks tracks) { histos.fill(HIST(kCorrType[corrType]) + HIST(kEvntType[evntType]) + HIST("hMult"), tracks.size()); for (auto const& iTrk : tracks) { auto phi = iTrk.phi(); - if constexpr (corrType == kFT0AMFT) { - o2::math_utils::bringTo02Pi(phi); - } + o2::math_utils::bringTo02Pi(phi); histos.fill(HIST(kCorrType[corrType]) + HIST(kEvntType[evntType]) + HIST("Trig_etavsphi"), phi, iTrk.eta()); histos.fill(HIST(kCorrType[corrType]) + HIST(kEvntType[evntType]) + HIST("Trig_eta"), iTrk.eta()); histos.fill(HIST(kCorrType[corrType]) + HIST(kEvntType[evntType]) + HIST("Trig_phi"), phi); @@ -420,6 +455,56 @@ struct LongrangeCorrelation { return eff; } + template + int getTrackPID(TTrack track) + { + // Computing Nsigma arrays for pion, kaon, and protons + std::array nSigmaTPC = {track.tpcNSigmaPi(), track.tpcNSigmaKa(), track.tpcNSigmaPr()}; + std::array nSigmaTOF = {track.tofNSigmaPi(), track.tofNSigmaKa(), track.tofNSigmaPr()}; + std::array nSigmaITS = {itsResponse.nSigmaITS(track), itsResponse.nSigmaITS(track), itsResponse.nSigmaITS(track)}; + int pid = -1; + + std::array nSigmaToUse = isUseItsPid ? nSigmaITS : nSigmaTPC; // Choose which nSigma to use: TPC or ITS + std::vector detectorNsigmaCut = isUseItsPid ? itsNsigmaCut : tpcNsigmaCut; // Choose which nSigma to use: TPC or ITS + + bool isPion, isKaon, isProton; + bool isDetectedPion = nSigmaToUse[0] < detectorNsigmaCut[0] && nSigmaToUse[0] > detectorNsigmaCut[0 + 3]; + bool isDetectedKaon = nSigmaToUse[1] < detectorNsigmaCut[1] && nSigmaToUse[1] > detectorNsigmaCut[1 + 3]; + bool isDetectedProton = nSigmaToUse[2] < detectorNsigmaCut[2] && nSigmaToUse[2] > detectorNsigmaCut[2 + 3]; + + bool isTofPion = nSigmaTOF[0] < tofNsigmaCut[0] && nSigmaTOF[0] > tofNsigmaCut[0 + 3]; + bool isTofKaon = nSigmaTOF[1] < tofNsigmaCut[1] && nSigmaTOF[1] > tofNsigmaCut[1 + 3]; + bool isTofProton = nSigmaTOF[2] < tofNsigmaCut[2] && nSigmaTOF[2] > tofNsigmaCut[2 + 3]; + + if (track.pt() > cfgTofPidPtCut && !track.hasTOF()) { + return 0; + } else if (track.pt() > cfgTofPidPtCut && track.hasTOF()) { + isPion = isTofPion && isDetectedPion; + isKaon = isTofKaon && isDetectedKaon; + isProton = isTofProton && isDetectedProton; + } else { + isPion = isDetectedPion; + isKaon = isDetectedKaon; + isProton = isDetectedProton; + } + + if ((isPion && isKaon) || (isPion && isProton) || (isKaon && isProton)) { + return 0; // more than one particle satisfy the criteria + } + + if (isPion) { + pid = PIONS; + } else if (isKaon) { + pid = KAONS; + } else if (isProton) { + pid = PROTONS; + } else { + return 0; // no particle satisfies the criteria + } + + return pid + 1; // shift the pid by 1, 1 = pion, 2 = kaon, 3 = proton + } + template void fillCorrFt0aGlobal(TTarget target, TTriggers const& triggers, TFT0s const& ft0, bool mixing, float vz, float multiplicity) { @@ -428,6 +513,8 @@ struct LongrangeCorrelation { histos.fill(HIST("Ft0aGlobal/SE/hMult_used"), multiplicity); for (auto const& triggerTrack : triggers) { + if (cfgTrackPid && getTrackPID(triggerTrack) != cfgTrackPid) + continue; // if PID is selected, check if the track has the right PID float trkeffw = 1.0f; if (isUseEffCorr) trkeffw = getTrkEffCorr(triggerTrack.eta(), triggerTrack.pt(), vz); @@ -474,6 +561,8 @@ struct LongrangeCorrelation { if (!mixing) histos.fill(HIST("Ft0cGlobal/SE/hMult_used"), multiplicity); for (auto const& triggerTrack : triggers) { + if (cfgTrackPid && getTrackPID(triggerTrack) != cfgTrackPid) + continue; // if PID is selected, check if the track has the right PID float trkeffw = 1.0f; if (isUseEffCorr) trkeffw = getTrkEffCorr(triggerTrack.eta(), triggerTrack.pt(), vz); @@ -520,6 +609,8 @@ struct LongrangeCorrelation { if (!mixing) histos.fill(HIST("MftGlobal/SE/hMult_used"), multiplicity); for (auto const& triggerTrack : triggers) { + if (cfgTrackPid && getTrackPID(triggerTrack) != cfgTrackPid) + continue; // if PID is selected, check if the track has the right PID float trkeffw = 1.0f; if (isUseEffCorr) trkeffw = getTrkEffCorr(triggerTrack.eta(), triggerTrack.pt(), vz); @@ -691,7 +782,7 @@ struct LongrangeCorrelation { auto bc = col.bc_as(); loadEffCorrection(bc.timestamp()); loadMultCorrection(bc.timestamp()); - fillYield(tracks); + fillYieldTpc(tracks); const auto& ft0 = col.foundFT0(); auto multiplicity = col.multiplicity(); float multw = getMultEffCorr(multiplicity); @@ -713,7 +804,7 @@ struct LongrangeCorrelation { auto bc = col.bc_as(); loadEffCorrection(bc.timestamp()); loadMultCorrection(bc.timestamp()); - fillYield(tracks); + fillYieldTpc(tracks); const auto& ft0 = col.foundFT0(); auto multiplicity = col.multiplicity(); float multw = getMultEffCorr(multiplicity); @@ -734,7 +825,7 @@ struct LongrangeCorrelation { auto bc = col.bc_as(); loadEffCorrection(bc.timestamp()); loadMultCorrection(bc.timestamp()); - fillYield(tracks); + fillYieldTpc(tracks); auto multiplicity = col.multiplicity(); float multw = getMultEffCorr(multiplicity); if (isUseEffCorr) @@ -751,7 +842,7 @@ struct LongrangeCorrelation { return; } if (col.has_foundFT0()) { - fillYield(mfttracks); + fillYieldMft(mfttracks); auto bc = col.bc_as(); loadMultCorrection(bc.timestamp()); const auto& ft0 = col.foundFT0(); @@ -808,7 +899,7 @@ struct LongrangeCorrelation { loadEffCorrection(bc.timestamp()); loadMultCorrection(bc.timestamp()); auto slicedTriggerTracks = tracks.sliceBy(perColGlobal, col1.globalIndex()); - fillYield(slicedTriggerTracks); + fillYieldTpc(slicedTriggerTracks); const auto& ft0 = col2.foundFT0(); auto multiplicity = col1.multiplicity(); float multw = getMultEffCorr(multiplicity); @@ -843,7 +934,7 @@ struct LongrangeCorrelation { loadEffCorrection(bc.timestamp()); loadMultCorrection(bc.timestamp()); auto slicedTriggerTracks = tracks.sliceBy(perColGlobal, col1.globalIndex()); - fillYield(slicedTriggerTracks); + fillYieldTpc(slicedTriggerTracks); const auto& ft0 = col2.foundFT0(); auto multiplicity = col1.multiplicity(); float multw = getMultEffCorr(multiplicity); @@ -906,7 +997,7 @@ struct LongrangeCorrelation { auto bc = col1.bc_as(); loadMultCorrection(bc.timestamp()); auto slicedTriggerMftTracks = mfttracks.sliceBy(perColMft, col1.globalIndex()); - fillYield(slicedTriggerMftTracks); + fillYieldMft(slicedTriggerMftTracks); const auto& ft0 = col2.foundFT0(); auto multiplicity = col1.multiplicity(); float multw = getMultEffCorr(multiplicity); From fff91edddf40aa1cd8edef1ff4235da4cb2e57ba Mon Sep 17 00:00:00 2001 From: abmodak <67369858+abmodak@users.noreply.github.com> Date: Wed, 8 Oct 2025 19:43:26 +0200 Subject: [PATCH 5/5] Fix o2 linter error --- PWGCF/TwoParticleCorrelations/Tasks/longrangeCorrelation.cxx | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/PWGCF/TwoParticleCorrelations/Tasks/longrangeCorrelation.cxx b/PWGCF/TwoParticleCorrelations/Tasks/longrangeCorrelation.cxx index 64005e7d6c6..7f01edb08b7 100644 --- a/PWGCF/TwoParticleCorrelations/Tasks/longrangeCorrelation.cxx +++ b/PWGCF/TwoParticleCorrelations/Tasks/longrangeCorrelation.cxx @@ -146,6 +146,7 @@ struct LongrangeCorrelation { Configurable isApplyGoodZvtxFT0vsPV{"isApplyGoodZvtxFT0vsPV", false, "Enable GoodZvtxFT0vsPV cut"}; Configurable isReadoutCenter{"isReadoutCenter", false, "Enable Readout Center"}; Configurable isUseEffCorr{"isUseEffCorr", false, "Enable efficiency correction"}; + Configurable cfgLowEffCut{"cfgLowEffCut", 0.001f, "Low efficiency cut"}; Configurable isUseItsPid{"isUseItsPid", false, "Use ITS PID for particle identification"}; Configurable cfgTofPidPtCut{"cfgTofPidPtCut", 0.3f, "Minimum pt to use TOF N-sigma"}; Configurable cfgTrackPid{"cfgTrackPid", 0, "1 = pion, 2 = kaon, 3 = proton, 0 for no PID"}; @@ -434,7 +435,7 @@ struct LongrangeCorrelation { } else { eff = 1.0; } - if (eff < 0.001) + if (eff < cfgLowEffCut) eff = 1.0; return eff; @@ -449,7 +450,7 @@ struct LongrangeCorrelation { } else { eff = 1.0; } - if (eff < 0.001) + if (eff < cfgLowEffCut) eff = 1.0; return eff;