From 822b50a479ece2c7af9d004798b829b94eadceee Mon Sep 17 00:00:00 2001 From: wulbyu Date: Thu, 23 Apr 2026 15:37:31 +0200 Subject: [PATCH 01/16] annalysis task for the v0 pt of the charm-bulk --- PWGHF/D2H/Tasks/CMakeLists.txt | 5 + PWGHF/D2H/Tasks/taskV0PtCharmBulk.cxx | 335 ++++++++++++++++++++++++++ 2 files changed, 340 insertions(+) create mode 100644 PWGHF/D2H/Tasks/taskV0PtCharmBulk.cxx diff --git a/PWGHF/D2H/Tasks/CMakeLists.txt b/PWGHF/D2H/Tasks/CMakeLists.txt index 608145a8e51..13f12fcdabb 100644 --- a/PWGHF/D2H/Tasks/CMakeLists.txt +++ b/PWGHF/D2H/Tasks/CMakeLists.txt @@ -158,3 +158,8 @@ o2physics_add_dpl_workflow(task-hidden-charm SOURCES taskHiddenCharm.cxx PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(task-v0-pt-charm-bulk + SOURCES taskV0PtCharmBulk.cxx + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2Physics::EventFilteringUtils + COMPONENT_NAME Analysis) diff --git a/PWGHF/D2H/Tasks/taskV0PtCharmBulk.cxx b/PWGHF/D2H/Tasks/taskV0PtCharmBulk.cxx new file mode 100644 index 00000000000..9d5f07ef0ce --- /dev/null +++ b/PWGHF/D2H/Tasks/taskV0PtCharmBulk.cxx @@ -0,0 +1,335 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file taskV0PtCharmBulk.cxx +/// \brief v0 pt for the charm-bulk correlation analysis task +/// \author Wu Chuntai, UNIPD, CCNU, and INFN Padova +/// \author Andrea Rossi, INFN Padova + +#include "PWGHF/Core/CentralityEstimation.h" +#include "PWGHF/Core/HfHelper.h" +#include "PWGHF/DataModel/AliasTables.h" +#include "PWGHF/DataModel/CandidateReconstructionTables.h" +#include "PWGHF/DataModel/CandidateSelectionTables.h" +#include "PWGHF/Utils/utilsEvSelHf.h" + +#include "Common/Core/RecoDecay.h" +#include "Common/DataModel/Centrality.h" +#include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/Multiplicity.h" +#include "Common/DataModel/TrackSelectionTables.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include +#include +#include +#include +#include + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; +using namespace o2::hf_centrality; +using namespace o2::hf_evsel; + +enum DecayChannel { + D0ToPiK = 0, + D0ToKPi +}; + +struct HfTaskV0PtCharmBulk +{ + // General configuration + Configurable etaAMin{"etaAMin", -0.8, "eta min for A subevent"}; + Configurable etaAMax{"etaAMax", -0.2, "eta max for A subevent"}; + Configurable etaBMin{"etaBMin", 0.2, "eta min for B subevent"}; + Configurable etaBMax{"etaBMax", 0.8, "eta max for B subevent"}; + Configurable centEstimator{"centEstimator", 2, "Centrality estimation (FT0A: 1, FT0C: 2, FT0M: 3, FV0A: 4)"}; + Configurable ccdbUrl{"ccdbUrl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; + Configurable> classMl{"classMl", {0, 2}, "Indexes of BDT scores to be stored. Two indexes max."}; + + // Track configuration + Configurable tpcNClsCrossedRowsMin{"tpcNClsCrossedRowsMin", 70, "min. TPC crossed rows for associated tracks"}; + Configurable etaTrkMax{"etaTrkMax", 1., "max. track eta"}; + Configurable ptTrkMin{"ptTrkMin", 0.2, "min. track pT"}; + Configurable ptTrkMax{"ptTrkMax", 5., "max. track pT"}; + Configurable dcaXYTrkMax{"dcaXYTrkMax", 1., "max. track DCA XY"}; + Configurable dcaZTrkMax{"dcaZTrkMax", 1., "max. track DCA Z"}; + Configurable usePtDiffDcaXYCut{"usePtDiffDcaXYCut", true, "Use pt-differential DCAxy cut for associated tracks"}; + Configurable dcaXYTrkNSigmaMax{"dcaXYTrkNSigmaMax", 7, "Cut on number of sigma deviations from expected DCA in the transverse direction"}; + Configurable dcaXYPtPrimTrkFunc{"dcaXYPtPrimTrkFunc", "(0.0026+0.005/(x^1.01))", "Functional form of pt-dependent DCAxy cut"}; + + // Candidate configuration + Configurable selectionFlag{"selectionFlag", 1, "Selection Flag for hadron (ML score tables are required to run the task)"}; + + // Event configuration + Configurable forceCharmInCollision{"forceCharmInCollision", true, "Flag to force charm in collision"}; + HfEventSelection hfEvSel; // event selection and monitoring + + o2::framework::Service ccdb{}; + SliceCache cache; + TF1* funcDcaXYPtCutPrimTrk = nullptr; + + // Subcribe and join the tables + using TracksTable = soa::Filtered>; + using D0CandTable = soa::Filtered>; + using CollsWithCentMult = soa::Join; + + // Select tracks and candidates + Filter filterSelectTracks = (nabs(aod::track::eta) < etaTrkMax) && (aod::track::pt > ptTrkMin) && (aod::track::pt < ptTrkMax) && (nabs(aod::track::dcaXY) < dcaXYTrkMax) && (nabs(aod::track::dcaZ) < dcaZTrkMax); + Filter filterSelectD0Candidates = aod::hf_sel_candidate_d0::isSelD0 >= selectionFlag || aod::hf_sel_candidate_d0::isSelD0bar >= selectionFlag; + + // pre-slice by collision for tracks + Preslice tracksTablePerColl = aod::track::collisionId; + Preslice candD0TablePerColl = aod::hf_cand::collisionId; + + // Partitions for selected candidates + Partition selectedD0ToPiK = aod::hf_sel_candidate_d0::isSelD0 >= selectionFlag; + Partition selectedD0ToKPi = aod::hf_sel_candidate_d0::isSelD0bar >= selectionFlag; + + // Partition selectedTracks = (nabs(aod::track::eta) < etaTrkMax) && (aod::track::pt > ptTrkMin) && (aod::track::pt < ptTrkMax) && (nabs(aod::track::dcaXY) < dcaXYTrkMax) && (nabs(aod::track::dcaZ) < dcaZTrkMax); + + // THnSparse configuration + ConfigurableAxis thnConfigAxisCent{"thnConfigAxisCent", {100, 0, 100}, ""}; + ConfigurableAxis thnConfigAxisCandMass{"thnConfigAxisCandMass", {200, 1.68, 2.08}, ""}; + ConfigurableAxis thnConfigAxisCandPt{"thnConfigAxisCandPt", {10, 0., 10.}, ""}; + ConfigurableAxis thnConfigAxisCandEta{"thnConfigAxisCandEta", {200, -1., 1.}, ""}; + ConfigurableAxis thnConfigAxisMPtTrkA{"thnConfigAxisMPtTrkA", {100, 0., 5.}, ""}; + ConfigurableAxis thnConfigAxisMPtTrkB{"thnConfigAxisMPtTrkB", {100, 0., 5.}, ""}; + ConfigurableAxis thnConfigAxisMlOne{"thnConfigAxisMlOne", {100, 0., 1.}, ""}; + ConfigurableAxis thnConfigAxisMlTwo{"thnConfigAxisMlTwo", {100, 0., 1.}, ""}; + + HistogramRegistry registry{"registry", {}}; + + void init(InitContext&) { + std::array doprocess{doprocessD0WCentMult}; + if ((std::accumulate(doprocess.begin(), doprocess.end(), 0)) == 0) { + LOGP(fatal, "At least one process function should be enabled at a time."); + } + + hfEvSel.addHistograms(registry); // collision monitoring + ccdb->setURL(ccdbUrl); + ccdb->setCaching(true); + ccdb->setLocalObjectValidityChecking(); + + // Define the axes for the THnSparse + const AxisSpec axisCent = {thnConfigAxisCent, "Centrality"}; + const AxisSpec axisCandMass = {thnConfigAxisCandMass, "Inv. mass (GeV/#it{c}^{2})"}; + const AxisSpec axisCandPt = {thnConfigAxisCandPt, "#it{p}_{T} (GeV/#it{c})"}; + const AxisSpec axisCandEta = {thnConfigAxisCandEta, "Eta"}; + const AxisSpec axisMPtTrkA = {thnConfigAxisMPtTrkA, "Mean pT of tracks in subevent A (GeV/#it{c})"}; + const AxisSpec axisMPtTrkB = {thnConfigAxisMPtTrkB, "Mean pT of tracks in subevent B (GeV/#it{c})"}; + const AxisSpec axisMlOne = {thnConfigAxisMlOne, "ML score 1"}; + const AxisSpec axisMlTwo = {thnConfigAxisMlTwo, "ML score 2"}; + + std::vector axes = {axisCent, axisCandMass, axisCandPt, axisCandEta, axisMPtTrkA, axisMPtTrkB, axisMlOne, axisMlTwo}; + registry.add("hSparseMeanTrkPtCharm", "THn for mean track pT", HistType::kTHnSparseF, axes); + + registry.add("hMeanTrkPtAVsCent", "Mean track pT A vs centrality", HistType::kTH2F, {axisCent, axisMPtTrkA}); + registry.add("hMeanTrkPtBVsCent", "Mean track pT B vs centrality", HistType::kTH2F, {axisCent, axisMPtTrkB}); + + if (usePtDiffDcaXYCut) { + funcDcaXYPtCutPrimTrk = new TF1("funcDcaXYPtCutPrimTrk", Form("[0]*%s", dcaXYPtPrimTrkFunc.value.data()), 0.001, 100); + funcDcaXYPtCutPrimTrk->SetParameter(0, dcaXYTrkNSigmaMax); + LOGF(info, "DCAxy pt-dependence function: %s", Form("[0]*%s", dcaXYPtPrimTrkFunc.value.data())); + } + } // End init + + /// Function to calculate mean pT of tracks in a given eta range (subevent) for a specific collision + /// \param tracks are the tracks to be used for mean pT calculation + /// \param candidates are the D0 candidates to be analyzed + /// \return a pair of pairs: {{meanPtA, countA}, {meanPtB, countB}} + template + std::pair calculateMeanPt(TracksTable const& tracks, CandT const& candidates) + { + float sumPtA = 0.f; + int countA = 0; + float sumPtB = 0.f; + int countB = 0; + std::vector candProngsA; + std::vector candProngsB; + + // Gather the global indices of the candidate daughters + if constexpr (std::is_same_v) + { + for (const auto& cand : candidates) { + if (cand.eta() > etaAMin && cand.eta() < etaAMax) { + candProngsA.push_back(cand.prong0Id()); + candProngsA.push_back(cand.prong1Id()); + } else if (cand.eta() > etaBMin && cand.eta() < etaBMax) { + candProngsB.push_back(cand.prong0Id()); + candProngsB.push_back(cand.prong1Id()); + } + } + } + + // Loop over tracks and calculate sum of pT and count for subevent A and B, excluding candidate daughters + for (const auto& track : tracks) + { + float eta = track.eta(); + float pt = track.pt(); + + // Select only global tracks with DCA information or with sufficient TPC crossed rows + if (track.isGlobalTrackWoDCA() || track.tpcNClsCrossedRows() < tpcNClsCrossedRowsMin) { + continue; + } + + // Apply DCA cuts + if (usePtDiffDcaXYCut) + { + float const dcaXYTrkCut = funcDcaXYPtCutPrimTrk->Eval(pt); + if (std::fabs(track.dcaXY()) > dcaXYTrkCut) + { + continue; + } + } + + int const trackGlobalIndex = track.globalIndex(); + // A side + if (eta < etaAMax && eta > etaAMin) + { + if (std::find(candProngsB.begin(), candProngsB.end(), trackGlobalIndex) != candProngsB.end()) + { + continue; // skip tracks that are daughters of the candidate in the opposite subevent + /*TODO: Considering remove the overlapped daughter tracks by recalculating the mean pT + excluding the candidate daughters on the fly for each candidate*/ + } + sumPtA += pt; + countA++; + } + + // B side + if (eta < etaBMax && eta > etaBMin) + { + if (std::find(candProngsA.begin(), candProngsA.end(), trackGlobalIndex) != candProngsA.end()) + { + continue; // skip tracks that are daughters of the candidate in the opposite subevent + /*TODO: Considering remove the overlapped daughter tracks by recalculating the mean pT + excluding the candidate daughters on the fly for each candidate*/ + } + sumPtB += pt; + countB++; + } + } + + if (countA == 0 && countB == 0) + { + return {std::numeric_limits::quiet_NaN(), std::numeric_limits::quiet_NaN()}; + // return NaN if no tracks in either subevent + } else if (countA == 0) + { + return {std::numeric_limits::quiet_NaN(), sumPtB / countB}; // NaN for subevent A if no tracks + } else if (countB == 0) + { + return {sumPtA / countA, std::numeric_limits::quiet_NaN()}; // NaN for subevent B if no tracks + } + return {sumPtA / countA, sumPtB / countB}; + } + + /// Calculate mean pT of tracks in subevent A and B, and fill the THnSparse + /// \param collision is the collision with the centrality and multiplicity information + /// \param tracks are the tracks to be used for mean pT calculation + /// \param candidates are the D0 candidates to be analyzed + template + void runCharmBulkAnalysis(CandT const& candidates, TracksTable const& tracks, float cent) + { + auto [meanPtA, meanPtB] = calculateMeanPt(tracks, candidates); + // Loop over candidates and fill the THnSparse + for (const auto& cand : candidates) + { + float invMass = 0.f; + std::vector outputMl = {-999., -999.}; + if constexpr (std::is_same_v) { + switch (Channel) { + case DecayChannel::D0ToPiK: + invMass = HfHelper::invMassD0ToPiK(cand); + for (unsigned int iclass = 0; iclass < classMl->size(); iclass++) { + outputMl[iclass] = cand.mlProbD0()[classMl->at(iclass)]; + } + break; + case DecayChannel::D0ToKPi: + invMass = HfHelper::invMassD0barToKPi(cand); + for (unsigned int iclass = 0; iclass < classMl->size(); iclass++) { + outputMl[iclass] = cand.mlProbD0bar()[classMl->at(iclass)]; + } + break; + } + } + registry.fill(HIST("hSparseMeanTrkPtCharm"), cent, invMass, cand.pt(), cand.eta(), meanPtA, meanPtB, outputMl[0], outputMl[1]); + registry.fill(HIST("hMeanTrkPtAVsCent"), cent, meanPtA); + registry.fill(HIST("hMeanTrkPtBVsCent"), cent, meanPtB); + } + } + + /// Check event selections for collision and fill event selection histograms + /// \param collision is the collision + template + bool isSelectedHfCollision(Coll const& collision, float& cent) + { + o2::hf_evsel::HfCollisionRejectionMask collRejMask{}; + if (centEstimator == CentralityEstimator::FT0A) + { + collRejMask = hfEvSel.getHfCollisionRejectionMask(collision, cent, ccdb, registry); + } else if (centEstimator == CentralityEstimator::FT0C) + { + collRejMask = hfEvSel.getHfCollisionRejectionMask(collision, cent, ccdb, registry); + } else if (centEstimator == CentralityEstimator::FT0M) + { + collRejMask = hfEvSel.getHfCollisionRejectionMask(collision, cent, ccdb, registry); + } else if (centEstimator == CentralityEstimator::FV0A) + { + collRejMask = hfEvSel.getHfCollisionRejectionMask(collision, cent, ccdb, registry); + } else + { + LOG(fatal) << "Centrality estimator not recognized for collision selection"; + std::abort(); + } + hfEvSel.fillHistograms(collision, collRejMask, cent); + return collRejMask == 0; + } + + void processD0WCentMult(CollsWithCentMult::iterator const& collision, TracksTable const& tracks, D0CandTable const& /* D0 candidates */) { + auto tableD0ToPiK = selectedD0ToPiK->sliceByCached(aod::hf_cand::collisionId, collision.globalIndex(), cache); + auto tableD0ToKPi = selectedD0ToKPi->sliceByCached(aod::hf_cand::collisionId, collision.globalIndex(), cache); + if (forceCharmInCollision && tableD0ToPiK.size() < 1 && tableD0ToKPi.size() < 1) { + return; + } + float cent = -1.f; + if (!isSelectedHfCollision(collision, cent)) { + return; + } + runCharmBulkAnalysis(tableD0ToPiK, tracks, cent); + runCharmBulkAnalysis(tableD0ToKPi, tracks, cent); + } + PROCESS_SWITCH(HfTaskV0PtCharmBulk, processD0WCentMult, "Process D0 candidates for tracks' mean-pT analysis", true); +}; // End struct HfTaskV0PtCharmBulk + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{adaptAnalysisTask(cfgc)}; +} From 0cc8b1dcf2988b55b612b9326fdc0e8a95bc996c Mon Sep 17 00:00:00 2001 From: wulbyu Date: Tue, 28 Apr 2026 14:10:25 +0200 Subject: [PATCH 02/16] replace the top. selections of track by filter Co-authored-by: Copilot --- PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx b/PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx index d0695616c84..187cc974c0c 100644 --- a/PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx +++ b/PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx @@ -107,6 +107,7 @@ struct HfTaskPtFlucCharmHadrons { Filter filterSelectDplusCandidates = aod::hf_sel_candidate_dplus::isSelDplusToPiKPi >= selectionFlag; Filter filterSelectD0Candidates = aod::hf_sel_candidate_d0::isSelD0 >= selectionFlag || aod::hf_sel_candidate_d0::isSelD0bar >= selectionFlag; + Filter filterSelectTracks = (nabs(aod::track::eta) < 0.8f) && (aod::track::pt > ptTrkMin) && (aod::track::pt < ptTrkMax) && (nabs(aod::track::dcaXY) < cfgCutDCAxy) && (nabs(aod::track::dcaZ) < cfgCutDCAz) && (aod::track::itsNCls >= cfgITScluster) && (aod::track::itsNClsInnerBarrel >= cfgITSbarrel) && (aod::track::tpcNClsFound >= cfgTPCcluster); Partition selectedD0ToPiKWMl = aod::hf_sel_candidate_d0::isSelD0 >= selectionFlag; Partition selectedD0ToKPiWMl = aod::hf_sel_candidate_d0::isSelD0bar >= selectionFlag; @@ -191,9 +192,9 @@ struct HfTaskPtFlucCharmHadrons { } template - bool selectionTrack(const T& candidate) const + bool selectionTrack(const T& track) const { - if (!(candidate.isGlobalTrack() && candidate.isPVContributor() && candidate.itsNCls() > cfgITScluster.value && candidate.tpcNClsFound() > cfgTPCcluster.value && candidate.itsNClsInnerBarrel() >= cfgITSbarrel.value && std::abs(candidate.dcaXY()) < cfgCutDCAxy.value && std::abs(candidate.dcaZ()) < cfgCutDCAz.value)) { + if (!(track.isGlobalTrack() && track.isPVContributor())) { return false; } return true; From 34cef08e7ca680d6da693481d743e4f924a74442 Mon Sep 17 00:00:00 2001 From: wulbyu Date: Tue, 28 Apr 2026 16:50:13 +0200 Subject: [PATCH 03/16] correlation of charm and bulk Co-authored-by: Copilot --- PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx | 237 +++++++++++++++------ 1 file changed, 166 insertions(+), 71 deletions(-) diff --git a/PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx b/PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx index 187cc974c0c..31c585d657a 100644 --- a/PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx +++ b/PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx @@ -64,6 +64,7 @@ struct HfTaskPtFlucCharmHadrons { Configurable centralityMax{"centralityMax", 100., "Maximum centrality accepted in SP computation"}; Configurable ccdbUrl{"ccdbUrl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; Configurable> classMl{"classMl", {0, 2}, "Indices of BDT scores to be stored. Two indexes max."}; + Configurable saveCharmBulkCorrelations{"saveCharmBulkCorrelations", false, "Flag to save charm-bulk correlations in THnSparse"}; Configurable cfgITScluster{"cfgITScluster", 5, "Number of ITS cluster"}; Configurable cfgITSbarrel{"cfgITSbarrel", 1, "Number of ITS inner barrel"}; @@ -77,6 +78,7 @@ struct HfTaskPtFlucCharmHadrons { Configurable etaAMax{"etaAMax", 0.0f, "A: negative eta max (0)"}; Configurable etaBMin{"etaBMin", 0.4f, "B: positive eta min (eta_min)"}; // or set from etaMinGap in code Configurable etaBMax{"etaBMax", 0.8f, "B: positive eta max"}; + Configurable etaTrkMax{"etaTrkMax", 0.8f, "Track eta max"}; Configurable ptTrkMin{"ptTrkMin", 0.2f, "Track pT min for (charged hadrons)"}; Configurable ptTrkMax{"ptTrkMax", 2.0f, "Track pT max for (charged hadrons)"}; @@ -107,7 +109,7 @@ struct HfTaskPtFlucCharmHadrons { Filter filterSelectDplusCandidates = aod::hf_sel_candidate_dplus::isSelDplusToPiKPi >= selectionFlag; Filter filterSelectD0Candidates = aod::hf_sel_candidate_d0::isSelD0 >= selectionFlag || aod::hf_sel_candidate_d0::isSelD0bar >= selectionFlag; - Filter filterSelectTracks = (nabs(aod::track::eta) < 0.8f) && (aod::track::pt > ptTrkMin) && (aod::track::pt < ptTrkMax) && (nabs(aod::track::dcaXY) < cfgCutDCAxy) && (nabs(aod::track::dcaZ) < cfgCutDCAz) && (aod::track::itsNCls >= cfgITScluster) && (aod::track::itsNClsInnerBarrel >= cfgITSbarrel) && (aod::track::tpcNClsFound >= cfgTPCcluster); + Filter filterSelectTracks = (nabs(aod::track::eta) < etaTrkMax) && (aod::track::pt > ptTrkMin) && (aod::track::pt < ptTrkMax) && (nabs(aod::track::dcaXY) < cfgCutDCAxy) && (nabs(aod::track::dcaZ) < cfgCutDCAz) && (aod::track::itsNCls >= cfgITScluster) && (aod::track::itsNClsInnerBarrel >= cfgITSbarrel) && (aod::track::tpcNClsFound >= cfgTPCcluster); Partition selectedD0ToPiKWMl = aod::hf_sel_candidate_d0::isSelD0 >= selectionFlag; Partition selectedD0ToKPiWMl = aod::hf_sel_candidate_d0::isSelD0bar >= selectionFlag; @@ -119,8 +121,15 @@ struct HfTaskPtFlucCharmHadrons { ConfigurableAxis axisPt{"axisPt", {VARIABLE_WIDTH, 0.2, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0, 6.5, 8.0, 10.0}, "Candidate pT axis"}; ConfigurableAxis axisCent{"axisCent", {VARIABLE_WIDTH, 0.0, 10.0, 40.0, 80.0}, "Centrality axis"}; ConfigurableAxis axisSign{"axisSign", {VARIABLE_WIDTH, -1.0, 1.0}, "Sign axis"}; - ConfigurableAxis thnConfigAxisMlOne{"thnConfigAxisMlOne", {1000, 0., 1.}, ""}; - ConfigurableAxis thnConfigAxisMlTwo{"thnConfigAxisMlTwo", {1000, 0., 1.}, ""}; + ConfigurableAxis axisMlOne{"axisMlOne", {1000, 0., 1.}, ""}; + ConfigurableAxis axisMlTwo{"axisMlTwo", {100, 0., 1.}, ""}; + ConfigurableAxis axisCandEta{"axisCandEta", {20, -1., 1.}, ""}; + ConfigurableAxis axisMPtTrkA{"axisMPtTrkA", {50, 0., 5.}, ""}; + ConfigurableAxis axisMPtTrkB{"axisMPtTrkB", {50, 0., 5.}, ""}; + ConfigurableAxis axisPtCandProduct{"axisPtCandProduct", {50, 0., 50.}, ""}; + ConfigurableAxis axisPtTrkProduct{"axisPtTrkProduct", {50, 0., 25.}, ""}; + ConfigurableAxis axisNTrkA{"axisNTrkA", {2000, 0., 2000.}, ""}; + ConfigurableAxis axisNTrkB{"axisNTrkB", {2000, 0., 2000.}, ""}; HistogramRegistry registry{"registry", {}, OutputObjHandlingPolicy::AnalysisObject}; @@ -140,8 +149,13 @@ struct HfTaskPtFlucCharmHadrons { const AxisSpec aPt{axisPt, "#it{p}_{T} (GeV/#it{c})"}; const AxisSpec aCent{axisCent, "Centrality"}; const AxisSpec aSign{axisSign, "Sign"}; - const AxisSpec thnAxisMlOne{thnConfigAxisMlOne, "Bkg score"}; - const AxisSpec thnAxisMlTwo{thnConfigAxisMlTwo, "FD score"}; + const AxisSpec aMlOne{axisMlOne, "Bkg score"}; + const AxisSpec aMlTwo{axisMlTwo, "FD score"}; + const AxisSpec aCandEta{axisCandEta, "Eta"}; + const AxisSpec aMPtTrkA{axisMPtTrkA, "Mean pT of tracks in subevent A (GeV/#it{c})"}; + const AxisSpec aMPtTrkB{axisMPtTrkB, "Mean pT of tracks in subevent B (GeV/#it{c})"}; + const AxisSpec aPtCandProduct{axisPtCandProduct, "#it{p}_{T}^{cand} * <#it{p}_{T}>^{trks} (GeV^{2}/#it{c}^{2})"}; + const AxisSpec aPtTrkProduct{axisPtTrkProduct, "#it{p}_{T}^{trks}(A) * #it{p}_{T}^{trks}(B) (GeV^{2}/#it{c}^{2})"}; // Event-level accumulators (charged hadrons only!) registry.add("hEvents", "events vs cent", HistType::kTH1F, {aCent}, true); @@ -156,15 +170,24 @@ struct HfTaskPtFlucCharmHadrons { registry.add("pNB", " vs cent (charged hadrons)", HistType::kTProfile, {aCent}, true); // Candidate mass distributions - registry.add("hD_mass", "D candidate mass (weight=1)", HistType::kTHnSparseF, {aInvMass, aCent, aPt, aSign, thnAxisMlOne, thnAxisMlTwo}, true); - registry.add("hD_f", "hD_f", HistType::kTHnSparseF, {aInvMass, aCent, aPt, aSign, thnAxisMlOne, thnAxisMlTwo}, true); - registry.add("hD_f_wPtB", "hD_f_wPtB", HistType::kTHnSparseF, {aInvMass, aCent, aPt, aSign, thnAxisMlOne, thnAxisMlTwo}, true); + registry.add("hD_mass", "D candidate mass (weight=1)", HistType::kTHnSparseF, {aInvMass, aCent, aPt, aSign, aMlOne, aMlTwo}, true); + registry.add("hD_f", "hD_f", HistType::kTHnSparseF, {aInvMass, aCent, aPt, aSign, aMlOne, aMlTwo}, true); + registry.add("hD_f_wPtB", "hD_f_wPtB", HistType::kTHnSparseF, {aInvMass, aCent, aPt, aSign, aMlOne, aMlTwo}, true); + + // charm-bulk correlations (optional) + if (saveCharmBulkCorrelations) { + registry.add("hCharmBulkCorrelations", "Charm-bulk correlations", HistType::kTHnSparseF, {aInvMass, aCent, aPt, aSign, aMlOne, aMlTwo, aCandEta, aMPtTrkA, aMPtTrkB, aPtCandProduct}, true); + registry.add("hMeanPtTrkAllColls", "Mean pT of charged hadrons for all collisions", HistType::kTHnSparseF, {aCent, aMPtTrkA, aMPtTrkB, aPtTrkProduct, aNTrkA, aNTrkB}, true); + } ccdb->setURL(ccdbUrl); ccdb->setCaching(true); ccdb->setLocalObjectValidityChecking(); }; // end init + // --------------------------------------------------------------------------- + // Helper functions + // --------------------------------------------------------------------------- /// Get the centrality /// \param collision is the collision with the centrality information double getCentrality(Colls::iterator const& collision) @@ -191,6 +214,22 @@ struct HfTaskPtFlucCharmHadrons { return cent; } + /// Get candidate mass + template + std::pair getCandMassAndSign(const CandT& cand, DecayChannel channel, Trk const& /*tracks*/) + { + if constexpr (std::is_same_v) { + return {HfHelper::invMassDplusToPiKPi(cand), cand.prong0_as().sign()}; + } else if constexpr (std::is_same_v) { + if (channel == DecayChannel::D0ToPiK) { + return {HfHelper::invMassD0ToPiK(cand), candidate.isSelD0bar() ? 3 : 1}; // 3: reflected D0bar, 1: pure D0 excluding reflected D0bar + } else if (channel == DecayChannel::D0ToKPi) { + return {HfHelper::invMassD0barToKPi(cand), candidate.isSelD0() ? 3 : 2}; // 3: reflected D0, 2: pure D0bar excluding reflected D0 + } + } + return {0., 0.}; // default return value for unsupported types + } + template bool selectionTrack(const T& track) const { @@ -200,9 +239,33 @@ struct HfTaskPtFlucCharmHadrons { return true; } - // ------------------------- + /// remove candidate daughters from the mean pT of tracks in A and B (if they are in the respective subevent) + template + float removeDaughtersFromMeanPt(const Cand& cand, float& meanPt, int& n, const std::vector& trkIDs) + { + if constexpr (Channel == DecayChannel::DplusToPiKPi) { + for (int iProng = 0; iProng < 3; ++iProng) { // for 3-prong + int trkID = cand.template prong_as(iProng).globalIndex(); + float ptDaughter = cand.template prong_as(iProng).pt(); + if (std::find(trkIDs.begin(), trkIDs.end(), trkID) != trkIDs.end()) { + meanPt = (meanPt * n - ptDaughter) / (n - 1); + } + } + } else if constexpr (Channel == DecayChannel::D0ToPiK || Channel == DecayChannel::D0ToKPi) { + for (int iProng = 0; iProng < 2; ++iProng) { // for 2-prong + int trkID = cand.template prong_as(iProng).globalIndex(); + float ptDaughter = cand.template prong_as(iProng).pt(); + if (std::find(trkIDs.begin(), trkIDs.end(), trkID) != trkIDs.end()) { + meanPt = (meanPt * n - ptDaughter) / (n - 1); + } + } + } + return meanPt; + } + + // --------------------------------------------------------------------------- // Event-wise mean pT from charged hadron tracks - // ------------------------- + // --------------------------------------------------------------------------- template std::pair computeMeanPtCharged(Trk const& tracks, float etaMin, float etaMax) const { @@ -218,9 +281,6 @@ struct HfTaskPtFlucCharmHadrons { continue; } float pt = trk.pt(); - if (pt < ptTrkMin.value || pt >= ptTrkMax.value) { - continue; - } sumPt += pt; ++n; } @@ -292,90 +352,125 @@ struct HfTaskPtFlucCharmHadrons { } // Compute in two eta-separated subevents using CHARGED tracks only - auto [meanPtA, nA] = computeMeanPtCharged(tracks, etaAMin.value, etaAMax.value); - auto [meanPtB, nB] = computeMeanPtCharged(tracks, etaBMin.value, etaBMax.value); + auto [RawMeanPtA, nA] = computeMeanPtCharged(tracks, etaAMin.value, etaAMax.value); + auto [RawMeanPtB, nB] = computeMeanPtCharged(tracks, etaBMin.value, etaBMax.value); // QA registry.fill(HIST("pNA"), cent, static_cast(nA)); registry.fill(HIST("pNB"), cent, static_cast(nB)); - if (!std::isfinite(meanPtA) || !std::isfinite(meanPtB)) { + if (!std::isfinite(RawMeanPtA) || !std::isfinite(RawMeanPtB)) { return; } registry.fill(HIST("hEvents"), cent, 1.f); - registry.fill(HIST("pMeanPtA"), cent, meanPtA); - registry.fill(HIST("pMeanPtB"), cent, meanPtB); - registry.fill(HIST("pMeanPtAmeanPtB"), cent, meanPtA * meanPtB); + registry.fill(HIST("pMeanPtA"), cent, RawMeanPtA); + registry.fill(HIST("pMeanPtB"), cent, RawMeanPtB); + registry.fill(HIST("pMeanPtAmeanPtB"), cent, RawMeanPtA * RawMeanPtB); + if (saveCharmBulkCorrelations) { + registry.fill(HIST("hMeanPtTrkAllColls"), cent, RawMeanPtA, RawMeanPtB, RawMeanPtA * RawMeanPtB, nA, nB); + + std::vector trkIDA; + std::vector trkIDB; + // collect track IDs in A and B + for (const auto& trk : tracks) { + if (!selectionTrack(trk)) { + continue; + } + float eta = trk.eta(); + if (eta > etaAMin.value && eta < etaAMax.value) { + trkIDA.push_back(trk.globalIndex()); + } else if (eta > etaBMin.value && eta < etaBMax.value) { + trkIDB.push_back(trk.globalIndex()); + } + } - int nDcandTotA = 0; + for (const auto& cand : candidates) + { + // apply ML selection + auto [ml1, ml2] = getMlScores(cand); + if (!passMlCut(ml1, ml2)) + { + continue; + } - for (const auto& cand : candidates) { - if (!passCandInA(cand)) { - continue; - } + // remove the daughters from the mean pT of tracks and calculate pt product + float eta = cand.eta(); + float candPtProduct{0.f}; + float meanPtA{0.f}; + float meanPtB{0.f}; + if (eta > etaAMin.value && eta < etaAMax.value) { + meanPtB = removeDaughtersFromMeanPt(cand, RawMeanPtB, nB, trkIDB); + meanPtA = RawMeanPtA; // no need to remove daughters from A if candidate is in A + candPtProduct = cand.pt() * meanPtB; + } else if (eta > etaBMin.value && eta < etaBMax.value) { + meanPtA = removeDaughtersFromMeanPt(cand, RawMeanPtA, nA, trkIDA); + meanPtB = RawMeanPtB; // no need to remove daughters from B if candidate is in B + candPtProduct = cand.pt() * meanPtA; + } - // compute ML exactly like numerator - const auto [ml1, ml2] = getMlScores(cand); + // get candidate mass and sign + auto [invMass, sign] = getCandMassAndSign(cand, Channel, tracks); - // apply same ML cut - if (!passMlCut(ml1, ml2)) { - continue; + // fill charm-bulk correlation thnsparse + registry.fill(HIST("hCharmBulkCorrelations"), invMass, cent, cand.pt(), sign, ml1, ml2, eta, meanPtA, meanPtB, candPtProduct); } - float pt = cand.pt(); - if (pt < ptTrkMinD.value || pt >= ptTrkMaxD.value) { - continue; - } - - ++nDcandTotA; - } + } else + { + int nDcandTotA = 0; + for (const auto& cand : candidates) { + if (!passCandInA(cand)) { + continue; + } - if (nDcandTotA <= 0) { - return; // cannot build fraction - } + // compute ML exactly like numerator + const auto [ml1, ml2] = getMlScores(cand); - registry.fill(HIST("hEventsD"), cent, 1.f); - const float invND = 1.f / static_cast(nDcandTotA); + // apply same ML cut + if (!passMlCut(ml1, ml2)) { + continue; + } + float pt = cand.pt(); + if (pt < ptTrkMinD.value || pt >= ptTrkMaxD.value) { + continue; + } - for (const auto& candidate : candidates) { - if (!passCandInA(candidate)) { - continue; + ++nDcandTotA; } - // ML first (same definition) - const auto [ml1, ml2] = getMlScores(candidate); - if (!passMlCut(ml1, ml2)) { - continue; + if (nDcandTotA <= 0) { + return; // cannot build fraction } - // compute mass - float massCand = 0.; - float signCand = 0.; - if constexpr (std::is_same_v) { - massCand = HfHelper::invMassDplusToPiKPi(candidate); - auto trackprong0 = candidate.template prong0_as(); - signCand = trackprong0.sign(); - } else { - if constexpr (Channel == DecayChannel::D0ToPiK) { - massCand = HfHelper::invMassD0ToPiK(candidate); - signCand = candidate.isSelD0bar() ? 3 : 1; // 3: reflected D0bar, 1: pure D0 excluding reflected D0bar - } else if constexpr (Channel == DecayChannel::D0ToKPi) { - massCand = HfHelper::invMassD0barToKPi(candidate); - signCand = candidate.isSelD0() ? 3 : 2; // 3: reflected D0, 2: pure D0bar excluding reflected D0 + registry.fill(HIST("hEventsD"), cent, 1.f); + const float invND = 1.f / static_cast(nDcandTotA); + + for (const auto& candidate : candidates) { + if (!passCandInA(candidate)) { + continue; } - } - const double ptCand = candidate.pt(); - if (ptCand < ptTrkMinD.value || ptCand >= ptTrkMaxD.value) { - continue; - } + // ML first (same definition) + const auto [ml1, ml2] = getMlScores(candidate); + if (!passMlCut(ml1, ml2)) { + continue; + } - registry.fill(HIST("hD_mass"), massCand, cent, ptCand, signCand, ml1, ml2, 1.f); - registry.fill(HIST("hD_f"), massCand, cent, ptCand, signCand, ml1, ml2, invND); - registry.fill(HIST("hD_f_wPtB"), massCand, cent, ptCand, signCand, ml1, ml2, meanPtB * invND); + // compute mass + auto [massCand, signCand] = getCandMassAndSign(candidate, Channel, tracks); + + const double ptCand = candidate.pt(); + if (ptCand < ptTrkMinD.value || ptCand >= ptTrkMaxD.value) { + continue; + } + + registry.fill(HIST("hD_mass"), massCand, cent, ptCand, signCand, ml1, ml2, 1.f); + registry.fill(HIST("hD_f"), massCand, cent, ptCand, signCand, ml1, ml2, invND); + registry.fill(HIST("hD_f_wPtB"), massCand, cent, ptCand, signCand, ml1, ml2, RawMeanPtB * invND); + } } - } + } // end runPtFlucAnalysis // D0 with ML void processD0Ml(Colls::iterator const& collision, From a7d8a4b787cb9563c7d8fb45ec2dec3e525ff9cf Mon Sep 17 00:00:00 2001 From: wulbyu Date: Wed, 29 Apr 2026 13:01:36 +0200 Subject: [PATCH 04/16] correlation of charm and bulk Co-authored-by: Copilot --- PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx | 108 +++++++++++---------- 1 file changed, 58 insertions(+), 50 deletions(-) diff --git a/PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx b/PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx index 31c585d657a..28411609283 100644 --- a/PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx +++ b/PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx @@ -109,8 +109,7 @@ struct HfTaskPtFlucCharmHadrons { Filter filterSelectDplusCandidates = aod::hf_sel_candidate_dplus::isSelDplusToPiKPi >= selectionFlag; Filter filterSelectD0Candidates = aod::hf_sel_candidate_d0::isSelD0 >= selectionFlag || aod::hf_sel_candidate_d0::isSelD0bar >= selectionFlag; - Filter filterSelectTracks = (nabs(aod::track::eta) < etaTrkMax) && (aod::track::pt > ptTrkMin) && (aod::track::pt < ptTrkMax) && (nabs(aod::track::dcaXY) < cfgCutDCAxy) && (nabs(aod::track::dcaZ) < cfgCutDCAz) && (aod::track::itsNCls >= cfgITScluster) && (aod::track::itsNClsInnerBarrel >= cfgITSbarrel) && (aod::track::tpcNClsFound >= cfgTPCcluster); - + Filter filterSelectTracks = (nabs(aod::track::eta) < etaTrkMax) && (aod::track::pt > ptTrkMin) && (aod::track::pt < ptTrkMax) && (nabs(aod::track::dcaXY) < cfgCutDCAxy) && (nabs(aod::track::dcaZ) < cfgCutDCAz); Partition selectedD0ToPiKWMl = aod::hf_sel_candidate_d0::isSelD0 >= selectionFlag; Partition selectedD0ToKPiWMl = aod::hf_sel_candidate_d0::isSelD0bar >= selectionFlag; @@ -156,6 +155,8 @@ struct HfTaskPtFlucCharmHadrons { const AxisSpec aMPtTrkB{axisMPtTrkB, "Mean pT of tracks in subevent B (GeV/#it{c})"}; const AxisSpec aPtCandProduct{axisPtCandProduct, "#it{p}_{T}^{cand} * <#it{p}_{T}>^{trks} (GeV^{2}/#it{c}^{2})"}; const AxisSpec aPtTrkProduct{axisPtTrkProduct, "#it{p}_{T}^{trks}(A) * #it{p}_{T}^{trks}(B) (GeV^{2}/#it{c}^{2})"}; + const AxisSpec aNTrkA{axisNTrkA, "N_{tracks} in subevent A"}; + const AxisSpec aNTrkB{axisNTrkB, "N_{tracks} in subevent B"}; // Event-level accumulators (charged hadrons only!) registry.add("hEvents", "events vs cent", HistType::kTH1F, {aCent}, true); @@ -180,6 +181,8 @@ struct HfTaskPtFlucCharmHadrons { registry.add("hMeanPtTrkAllColls", "Mean pT of charged hadrons for all collisions", HistType::kTHnSparseF, {aCent, aMPtTrkA, aMPtTrkB, aPtTrkProduct, aNTrkA, aNTrkB}, true); } + hfEvSel.addHistograms(registry); // collision monitoring + ccdb->setURL(ccdbUrl); ccdb->setCaching(true); ccdb->setLocalObjectValidityChecking(); @@ -188,30 +191,31 @@ struct HfTaskPtFlucCharmHadrons { // --------------------------------------------------------------------------- // Helper functions // --------------------------------------------------------------------------- - /// Get the centrality - /// \param collision is the collision with the centrality information - double getCentrality(Colls::iterator const& collision) + /// Check event selections for collision and fill event selection histograms, and revalculate centrality + /// \param collision is the collision + template + bool isSelectedHfCollision(Coll const& collision, float& cent) { - double cent = -999.; - switch (centEstimator) { - case CentralityEstimator::FV0A: - cent = collision.centFV0A(); - break; - case CentralityEstimator::FT0M: - cent = collision.centFT0M(); - break; - case CentralityEstimator::FT0A: - cent = collision.centFT0A(); - break; - case CentralityEstimator::FT0C: - cent = collision.centFT0C(); - break; - default: - LOG(warning) << "Centrality estimator not valid. Possible values are V0A, T0M, T0A, T0C. Fallback to V0A"; - cent = collision.centFV0A(); - break; - } - return cent; + o2::hf_evsel::HfCollisionRejectionMask collRejMask{}; + if (centEstimator == CentralityEstimator::FT0A) + { + collRejMask = hfEvSel.getHfCollisionRejectionMask(collision, cent, ccdb, registry); + } else if (centEstimator == CentralityEstimator::FT0C) + { + collRejMask = hfEvSel.getHfCollisionRejectionMask(collision, cent, ccdb, registry); + } else if (centEstimator == CentralityEstimator::FT0M) + { + collRejMask = hfEvSel.getHfCollisionRejectionMask(collision, cent, ccdb, registry); + } else if (centEstimator == CentralityEstimator::FV0A) + { + collRejMask = hfEvSel.getHfCollisionRejectionMask(collision, cent, ccdb, registry); + } else + { + LOG(fatal) << "Centrality estimator not recognized for collision selection"; + std::abort(); + } + hfEvSel.fillHistograms(collision, collRejMask, cent); + return collRejMask == 0; } /// Get candidate mass @@ -219,12 +223,12 @@ struct HfTaskPtFlucCharmHadrons { std::pair getCandMassAndSign(const CandT& cand, DecayChannel channel, Trk const& /*tracks*/) { if constexpr (std::is_same_v) { - return {HfHelper::invMassDplusToPiKPi(cand), cand.prong0_as().sign()}; - } else if constexpr (std::is_same_v) { + return {HfHelper::invMassDplusToPiKPi(cand), cand.template prong0_as().sign()}; + } else if constexpr (std::is_same_v) { if (channel == DecayChannel::D0ToPiK) { - return {HfHelper::invMassD0ToPiK(cand), candidate.isSelD0bar() ? 3 : 1}; // 3: reflected D0bar, 1: pure D0 excluding reflected D0bar + return {HfHelper::invMassD0ToPiK(cand), cand.isSelD0bar() ? 3 : 1}; // 3: reflected D0bar, 1: pure D0 excluding reflected D0bar } else if (channel == DecayChannel::D0ToKPi) { - return {HfHelper::invMassD0barToKPi(cand), candidate.isSelD0() ? 3 : 2}; // 3: reflected D0, 2: pure D0bar excluding reflected D0 + return {HfHelper::invMassD0barToKPi(cand), cand.isSelD0() ? 3 : 2}; // 3: reflected D0, 2: pure D0bar excluding reflected D0 } } return {0., 0.}; // default return value for unsupported types @@ -233,30 +237,31 @@ struct HfTaskPtFlucCharmHadrons { template bool selectionTrack(const T& track) const { - if (!(track.isGlobalTrack() && track.isPVContributor())) { + if (!(track.isGlobalTrack() && track.isPVContributor() && track.itsNCls() > cfgITScluster.value && track.tpcNClsFound() > cfgTPCcluster.value && track.itsNClsInnerBarrel() >= cfgITSbarrel.value)) { return false; } return true; } /// remove candidate daughters from the mean pT of tracks in A and B (if they are in the respective subevent) - template - float removeDaughtersFromMeanPt(const Cand& cand, float& meanPt, int& n, const std::vector& trkIDs) + template + float removeDaughtersFromMeanPt(const CandT& cand, const float& rawMeanPt, const int& n, const std::vector& trkIDs) { + float meanPt{0.f}; if constexpr (Channel == DecayChannel::DplusToPiKPi) { + std::array daugIDs = {cand.prong0Id(), cand.prong1Id(), cand.prong2Id()}; + std::array daugPts = {cand.ptProng0(), cand.ptProng1(), cand.ptProng2()}; for (int iProng = 0; iProng < 3; ++iProng) { // for 3-prong - int trkID = cand.template prong_as(iProng).globalIndex(); - float ptDaughter = cand.template prong_as(iProng).pt(); - if (std::find(trkIDs.begin(), trkIDs.end(), trkID) != trkIDs.end()) { - meanPt = (meanPt * n - ptDaughter) / (n - 1); + if (std::find(trkIDs.begin(), trkIDs.end(), daugIDs[iProng]) != trkIDs.end()) { + meanPt = (rawMeanPt * n - daugPts[iProng]) / (n - 1); } } } else if constexpr (Channel == DecayChannel::D0ToPiK || Channel == DecayChannel::D0ToKPi) { + std::array daugIDs = {cand.prong0Id(), cand.prong1Id()}; + std::array daugPts = {cand.ptProng0(), cand.ptProng1()}; for (int iProng = 0; iProng < 2; ++iProng) { // for 2-prong - int trkID = cand.template prong_as(iProng).globalIndex(); - float ptDaughter = cand.template prong_as(iProng).pt(); - if (std::find(trkIDs.begin(), trkIDs.end(), trkID) != trkIDs.end()) { - meanPt = (meanPt * n - ptDaughter) / (n - 1); + if (std::find(trkIDs.begin(), trkIDs.end(), daugIDs[iProng]) != trkIDs.end()) { + meanPt = (rawMeanPt * n - daugPts[iProng]) / (n - 1); } } } @@ -338,16 +343,17 @@ struct HfTaskPtFlucCharmHadrons { ml2 >= mlTwoMin.value && ml2 < mlTwoMax.value); } - /// Compute the scalar product - /// \param collision is the collision with the Q vector information and event plane - /// \param candidates are the selected candidates + /// Compute the mean pT of the event using charged tracks in two eta-separated subevents, and correlate with D candidates + /// \param collision is the collision with centrality and event selection + /// \param candidates are the D candidates to correlate with the mean pT of the event + /// \param tracks are the tracks used to compute the mean pT of the event template void runPtFlucAnalysis(Colls::iterator const& collision, T1 const& candidates, Trk const& tracks) { - float cent = getCentrality(collision); - if (cent < centralityMin || cent > centralityMax) { + float cent{0.f}; + if (!isSelectedHfCollision(collision, cent)) { return; } @@ -397,24 +403,26 @@ struct HfTaskPtFlucCharmHadrons { // remove the daughters from the mean pT of tracks and calculate pt product float eta = cand.eta(); + float pt = cand.pt(); + float candPtProduct{0.f}; float meanPtA{0.f}; float meanPtB{0.f}; if (eta > etaAMin.value && eta < etaAMax.value) { - meanPtB = removeDaughtersFromMeanPt(cand, RawMeanPtB, nB, trkIDB); + meanPtB = removeDaughtersFromMeanPt(cand, RawMeanPtB, nB, trkIDB); meanPtA = RawMeanPtA; // no need to remove daughters from A if candidate is in A - candPtProduct = cand.pt() * meanPtB; + candPtProduct = pt * meanPtB; } else if (eta > etaBMin.value && eta < etaBMax.value) { - meanPtA = removeDaughtersFromMeanPt(cand, RawMeanPtA, nA, trkIDA); + meanPtA = removeDaughtersFromMeanPt(cand, RawMeanPtA, nA, trkIDA); meanPtB = RawMeanPtB; // no need to remove daughters from B if candidate is in B - candPtProduct = cand.pt() * meanPtA; + candPtProduct = pt * meanPtA; } // get candidate mass and sign auto [invMass, sign] = getCandMassAndSign(cand, Channel, tracks); // fill charm-bulk correlation thnsparse - registry.fill(HIST("hCharmBulkCorrelations"), invMass, cent, cand.pt(), sign, ml1, ml2, eta, meanPtA, meanPtB, candPtProduct); + registry.fill(HIST("hCharmBulkCorrelations"), invMass, cent, pt, sign, ml1, ml2, eta, meanPtA, meanPtB, candPtProduct); } } else { From 9d2c9aa7e3774b302b1e6e0e368045f32e4471a2 Mon Sep 17 00:00:00 2001 From: wulbyu Date: Wed, 29 Apr 2026 13:08:33 +0200 Subject: [PATCH 05/16] remove the old task --- PWGHF/D2H/Tasks/CMakeLists.txt | 5 - PWGHF/D2H/Tasks/taskV0PtCharmBulk.cxx | 335 -------------------------- 2 files changed, 340 deletions(-) delete mode 100644 PWGHF/D2H/Tasks/taskV0PtCharmBulk.cxx diff --git a/PWGHF/D2H/Tasks/CMakeLists.txt b/PWGHF/D2H/Tasks/CMakeLists.txt index 13f12fcdabb..608145a8e51 100644 --- a/PWGHF/D2H/Tasks/CMakeLists.txt +++ b/PWGHF/D2H/Tasks/CMakeLists.txt @@ -158,8 +158,3 @@ o2physics_add_dpl_workflow(task-hidden-charm SOURCES taskHiddenCharm.cxx PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore COMPONENT_NAME Analysis) - -o2physics_add_dpl_workflow(task-v0-pt-charm-bulk - SOURCES taskV0PtCharmBulk.cxx - PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2Physics::EventFilteringUtils - COMPONENT_NAME Analysis) diff --git a/PWGHF/D2H/Tasks/taskV0PtCharmBulk.cxx b/PWGHF/D2H/Tasks/taskV0PtCharmBulk.cxx deleted file mode 100644 index 9d5f07ef0ce..00000000000 --- a/PWGHF/D2H/Tasks/taskV0PtCharmBulk.cxx +++ /dev/null @@ -1,335 +0,0 @@ -// Copyright 2019-2020 CERN and copyright holders of ALICE O2. -// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. -// All rights not expressly granted are reserved. -// -// This software is distributed under the terms of the GNU General Public -// License v3 (GPL Version 3), copied verbatim in the file "COPYING". -// -// In applying this license CERN does not waive the privileges and immunities -// granted to it by virtue of its status as an Intergovernmental Organization -// or submit itself to any jurisdiction. - -/// \file taskV0PtCharmBulk.cxx -/// \brief v0 pt for the charm-bulk correlation analysis task -/// \author Wu Chuntai, UNIPD, CCNU, and INFN Padova -/// \author Andrea Rossi, INFN Padova - -#include "PWGHF/Core/CentralityEstimation.h" -#include "PWGHF/Core/HfHelper.h" -#include "PWGHF/DataModel/AliasTables.h" -#include "PWGHF/DataModel/CandidateReconstructionTables.h" -#include "PWGHF/DataModel/CandidateSelectionTables.h" -#include "PWGHF/Utils/utilsEvSelHf.h" - -#include "Common/Core/RecoDecay.h" -#include "Common/DataModel/Centrality.h" -#include "Common/DataModel/EventSelection.h" -#include "Common/DataModel/Multiplicity.h" -#include "Common/DataModel/TrackSelectionTables.h" - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include - -#include -#include -#include -#include -#include -#include - -using namespace o2; -using namespace o2::framework; -using namespace o2::framework::expressions; -using namespace o2::hf_centrality; -using namespace o2::hf_evsel; - -enum DecayChannel { - D0ToPiK = 0, - D0ToKPi -}; - -struct HfTaskV0PtCharmBulk -{ - // General configuration - Configurable etaAMin{"etaAMin", -0.8, "eta min for A subevent"}; - Configurable etaAMax{"etaAMax", -0.2, "eta max for A subevent"}; - Configurable etaBMin{"etaBMin", 0.2, "eta min for B subevent"}; - Configurable etaBMax{"etaBMax", 0.8, "eta max for B subevent"}; - Configurable centEstimator{"centEstimator", 2, "Centrality estimation (FT0A: 1, FT0C: 2, FT0M: 3, FV0A: 4)"}; - Configurable ccdbUrl{"ccdbUrl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; - Configurable> classMl{"classMl", {0, 2}, "Indexes of BDT scores to be stored. Two indexes max."}; - - // Track configuration - Configurable tpcNClsCrossedRowsMin{"tpcNClsCrossedRowsMin", 70, "min. TPC crossed rows for associated tracks"}; - Configurable etaTrkMax{"etaTrkMax", 1., "max. track eta"}; - Configurable ptTrkMin{"ptTrkMin", 0.2, "min. track pT"}; - Configurable ptTrkMax{"ptTrkMax", 5., "max. track pT"}; - Configurable dcaXYTrkMax{"dcaXYTrkMax", 1., "max. track DCA XY"}; - Configurable dcaZTrkMax{"dcaZTrkMax", 1., "max. track DCA Z"}; - Configurable usePtDiffDcaXYCut{"usePtDiffDcaXYCut", true, "Use pt-differential DCAxy cut for associated tracks"}; - Configurable dcaXYTrkNSigmaMax{"dcaXYTrkNSigmaMax", 7, "Cut on number of sigma deviations from expected DCA in the transverse direction"}; - Configurable dcaXYPtPrimTrkFunc{"dcaXYPtPrimTrkFunc", "(0.0026+0.005/(x^1.01))", "Functional form of pt-dependent DCAxy cut"}; - - // Candidate configuration - Configurable selectionFlag{"selectionFlag", 1, "Selection Flag for hadron (ML score tables are required to run the task)"}; - - // Event configuration - Configurable forceCharmInCollision{"forceCharmInCollision", true, "Flag to force charm in collision"}; - HfEventSelection hfEvSel; // event selection and monitoring - - o2::framework::Service ccdb{}; - SliceCache cache; - TF1* funcDcaXYPtCutPrimTrk = nullptr; - - // Subcribe and join the tables - using TracksTable = soa::Filtered>; - using D0CandTable = soa::Filtered>; - using CollsWithCentMult = soa::Join; - - // Select tracks and candidates - Filter filterSelectTracks = (nabs(aod::track::eta) < etaTrkMax) && (aod::track::pt > ptTrkMin) && (aod::track::pt < ptTrkMax) && (nabs(aod::track::dcaXY) < dcaXYTrkMax) && (nabs(aod::track::dcaZ) < dcaZTrkMax); - Filter filterSelectD0Candidates = aod::hf_sel_candidate_d0::isSelD0 >= selectionFlag || aod::hf_sel_candidate_d0::isSelD0bar >= selectionFlag; - - // pre-slice by collision for tracks - Preslice tracksTablePerColl = aod::track::collisionId; - Preslice candD0TablePerColl = aod::hf_cand::collisionId; - - // Partitions for selected candidates - Partition selectedD0ToPiK = aod::hf_sel_candidate_d0::isSelD0 >= selectionFlag; - Partition selectedD0ToKPi = aod::hf_sel_candidate_d0::isSelD0bar >= selectionFlag; - - // Partition selectedTracks = (nabs(aod::track::eta) < etaTrkMax) && (aod::track::pt > ptTrkMin) && (aod::track::pt < ptTrkMax) && (nabs(aod::track::dcaXY) < dcaXYTrkMax) && (nabs(aod::track::dcaZ) < dcaZTrkMax); - - // THnSparse configuration - ConfigurableAxis thnConfigAxisCent{"thnConfigAxisCent", {100, 0, 100}, ""}; - ConfigurableAxis thnConfigAxisCandMass{"thnConfigAxisCandMass", {200, 1.68, 2.08}, ""}; - ConfigurableAxis thnConfigAxisCandPt{"thnConfigAxisCandPt", {10, 0., 10.}, ""}; - ConfigurableAxis thnConfigAxisCandEta{"thnConfigAxisCandEta", {200, -1., 1.}, ""}; - ConfigurableAxis thnConfigAxisMPtTrkA{"thnConfigAxisMPtTrkA", {100, 0., 5.}, ""}; - ConfigurableAxis thnConfigAxisMPtTrkB{"thnConfigAxisMPtTrkB", {100, 0., 5.}, ""}; - ConfigurableAxis thnConfigAxisMlOne{"thnConfigAxisMlOne", {100, 0., 1.}, ""}; - ConfigurableAxis thnConfigAxisMlTwo{"thnConfigAxisMlTwo", {100, 0., 1.}, ""}; - - HistogramRegistry registry{"registry", {}}; - - void init(InitContext&) { - std::array doprocess{doprocessD0WCentMult}; - if ((std::accumulate(doprocess.begin(), doprocess.end(), 0)) == 0) { - LOGP(fatal, "At least one process function should be enabled at a time."); - } - - hfEvSel.addHistograms(registry); // collision monitoring - ccdb->setURL(ccdbUrl); - ccdb->setCaching(true); - ccdb->setLocalObjectValidityChecking(); - - // Define the axes for the THnSparse - const AxisSpec axisCent = {thnConfigAxisCent, "Centrality"}; - const AxisSpec axisCandMass = {thnConfigAxisCandMass, "Inv. mass (GeV/#it{c}^{2})"}; - const AxisSpec axisCandPt = {thnConfigAxisCandPt, "#it{p}_{T} (GeV/#it{c})"}; - const AxisSpec axisCandEta = {thnConfigAxisCandEta, "Eta"}; - const AxisSpec axisMPtTrkA = {thnConfigAxisMPtTrkA, "Mean pT of tracks in subevent A (GeV/#it{c})"}; - const AxisSpec axisMPtTrkB = {thnConfigAxisMPtTrkB, "Mean pT of tracks in subevent B (GeV/#it{c})"}; - const AxisSpec axisMlOne = {thnConfigAxisMlOne, "ML score 1"}; - const AxisSpec axisMlTwo = {thnConfigAxisMlTwo, "ML score 2"}; - - std::vector axes = {axisCent, axisCandMass, axisCandPt, axisCandEta, axisMPtTrkA, axisMPtTrkB, axisMlOne, axisMlTwo}; - registry.add("hSparseMeanTrkPtCharm", "THn for mean track pT", HistType::kTHnSparseF, axes); - - registry.add("hMeanTrkPtAVsCent", "Mean track pT A vs centrality", HistType::kTH2F, {axisCent, axisMPtTrkA}); - registry.add("hMeanTrkPtBVsCent", "Mean track pT B vs centrality", HistType::kTH2F, {axisCent, axisMPtTrkB}); - - if (usePtDiffDcaXYCut) { - funcDcaXYPtCutPrimTrk = new TF1("funcDcaXYPtCutPrimTrk", Form("[0]*%s", dcaXYPtPrimTrkFunc.value.data()), 0.001, 100); - funcDcaXYPtCutPrimTrk->SetParameter(0, dcaXYTrkNSigmaMax); - LOGF(info, "DCAxy pt-dependence function: %s", Form("[0]*%s", dcaXYPtPrimTrkFunc.value.data())); - } - } // End init - - /// Function to calculate mean pT of tracks in a given eta range (subevent) for a specific collision - /// \param tracks are the tracks to be used for mean pT calculation - /// \param candidates are the D0 candidates to be analyzed - /// \return a pair of pairs: {{meanPtA, countA}, {meanPtB, countB}} - template - std::pair calculateMeanPt(TracksTable const& tracks, CandT const& candidates) - { - float sumPtA = 0.f; - int countA = 0; - float sumPtB = 0.f; - int countB = 0; - std::vector candProngsA; - std::vector candProngsB; - - // Gather the global indices of the candidate daughters - if constexpr (std::is_same_v) - { - for (const auto& cand : candidates) { - if (cand.eta() > etaAMin && cand.eta() < etaAMax) { - candProngsA.push_back(cand.prong0Id()); - candProngsA.push_back(cand.prong1Id()); - } else if (cand.eta() > etaBMin && cand.eta() < etaBMax) { - candProngsB.push_back(cand.prong0Id()); - candProngsB.push_back(cand.prong1Id()); - } - } - } - - // Loop over tracks and calculate sum of pT and count for subevent A and B, excluding candidate daughters - for (const auto& track : tracks) - { - float eta = track.eta(); - float pt = track.pt(); - - // Select only global tracks with DCA information or with sufficient TPC crossed rows - if (track.isGlobalTrackWoDCA() || track.tpcNClsCrossedRows() < tpcNClsCrossedRowsMin) { - continue; - } - - // Apply DCA cuts - if (usePtDiffDcaXYCut) - { - float const dcaXYTrkCut = funcDcaXYPtCutPrimTrk->Eval(pt); - if (std::fabs(track.dcaXY()) > dcaXYTrkCut) - { - continue; - } - } - - int const trackGlobalIndex = track.globalIndex(); - // A side - if (eta < etaAMax && eta > etaAMin) - { - if (std::find(candProngsB.begin(), candProngsB.end(), trackGlobalIndex) != candProngsB.end()) - { - continue; // skip tracks that are daughters of the candidate in the opposite subevent - /*TODO: Considering remove the overlapped daughter tracks by recalculating the mean pT - excluding the candidate daughters on the fly for each candidate*/ - } - sumPtA += pt; - countA++; - } - - // B side - if (eta < etaBMax && eta > etaBMin) - { - if (std::find(candProngsA.begin(), candProngsA.end(), trackGlobalIndex) != candProngsA.end()) - { - continue; // skip tracks that are daughters of the candidate in the opposite subevent - /*TODO: Considering remove the overlapped daughter tracks by recalculating the mean pT - excluding the candidate daughters on the fly for each candidate*/ - } - sumPtB += pt; - countB++; - } - } - - if (countA == 0 && countB == 0) - { - return {std::numeric_limits::quiet_NaN(), std::numeric_limits::quiet_NaN()}; - // return NaN if no tracks in either subevent - } else if (countA == 0) - { - return {std::numeric_limits::quiet_NaN(), sumPtB / countB}; // NaN for subevent A if no tracks - } else if (countB == 0) - { - return {sumPtA / countA, std::numeric_limits::quiet_NaN()}; // NaN for subevent B if no tracks - } - return {sumPtA / countA, sumPtB / countB}; - } - - /// Calculate mean pT of tracks in subevent A and B, and fill the THnSparse - /// \param collision is the collision with the centrality and multiplicity information - /// \param tracks are the tracks to be used for mean pT calculation - /// \param candidates are the D0 candidates to be analyzed - template - void runCharmBulkAnalysis(CandT const& candidates, TracksTable const& tracks, float cent) - { - auto [meanPtA, meanPtB] = calculateMeanPt(tracks, candidates); - // Loop over candidates and fill the THnSparse - for (const auto& cand : candidates) - { - float invMass = 0.f; - std::vector outputMl = {-999., -999.}; - if constexpr (std::is_same_v) { - switch (Channel) { - case DecayChannel::D0ToPiK: - invMass = HfHelper::invMassD0ToPiK(cand); - for (unsigned int iclass = 0; iclass < classMl->size(); iclass++) { - outputMl[iclass] = cand.mlProbD0()[classMl->at(iclass)]; - } - break; - case DecayChannel::D0ToKPi: - invMass = HfHelper::invMassD0barToKPi(cand); - for (unsigned int iclass = 0; iclass < classMl->size(); iclass++) { - outputMl[iclass] = cand.mlProbD0bar()[classMl->at(iclass)]; - } - break; - } - } - registry.fill(HIST("hSparseMeanTrkPtCharm"), cent, invMass, cand.pt(), cand.eta(), meanPtA, meanPtB, outputMl[0], outputMl[1]); - registry.fill(HIST("hMeanTrkPtAVsCent"), cent, meanPtA); - registry.fill(HIST("hMeanTrkPtBVsCent"), cent, meanPtB); - } - } - - /// Check event selections for collision and fill event selection histograms - /// \param collision is the collision - template - bool isSelectedHfCollision(Coll const& collision, float& cent) - { - o2::hf_evsel::HfCollisionRejectionMask collRejMask{}; - if (centEstimator == CentralityEstimator::FT0A) - { - collRejMask = hfEvSel.getHfCollisionRejectionMask(collision, cent, ccdb, registry); - } else if (centEstimator == CentralityEstimator::FT0C) - { - collRejMask = hfEvSel.getHfCollisionRejectionMask(collision, cent, ccdb, registry); - } else if (centEstimator == CentralityEstimator::FT0M) - { - collRejMask = hfEvSel.getHfCollisionRejectionMask(collision, cent, ccdb, registry); - } else if (centEstimator == CentralityEstimator::FV0A) - { - collRejMask = hfEvSel.getHfCollisionRejectionMask(collision, cent, ccdb, registry); - } else - { - LOG(fatal) << "Centrality estimator not recognized for collision selection"; - std::abort(); - } - hfEvSel.fillHistograms(collision, collRejMask, cent); - return collRejMask == 0; - } - - void processD0WCentMult(CollsWithCentMult::iterator const& collision, TracksTable const& tracks, D0CandTable const& /* D0 candidates */) { - auto tableD0ToPiK = selectedD0ToPiK->sliceByCached(aod::hf_cand::collisionId, collision.globalIndex(), cache); - auto tableD0ToKPi = selectedD0ToKPi->sliceByCached(aod::hf_cand::collisionId, collision.globalIndex(), cache); - if (forceCharmInCollision && tableD0ToPiK.size() < 1 && tableD0ToKPi.size() < 1) { - return; - } - float cent = -1.f; - if (!isSelectedHfCollision(collision, cent)) { - return; - } - runCharmBulkAnalysis(tableD0ToPiK, tracks, cent); - runCharmBulkAnalysis(tableD0ToKPi, tracks, cent); - } - PROCESS_SWITCH(HfTaskV0PtCharmBulk, processD0WCentMult, "Process D0 candidates for tracks' mean-pT analysis", true); -}; // End struct HfTaskV0PtCharmBulk - -WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) -{ - return WorkflowSpec{adaptAnalysisTask(cfgc)}; -} From bd1d3c4f932e1537bc8f24ad6dd813c6337a2da1 Mon Sep 17 00:00:00 2001 From: wulbyu Date: Wed, 29 Apr 2026 13:14:52 +0200 Subject: [PATCH 06/16] magic number Co-authored-by: Copilot --- PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx b/PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx index 28411609283..02edb798a8c 100644 --- a/PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx +++ b/PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx @@ -13,6 +13,8 @@ /// \brief Analysis task for charm hadron pt fluctuation /// /// \author Prottay Das, prottay.das@cern.ch +/// \author Wu Chuntai, UNIPD, CCNU, and INFN Padova +/// \author Andrea Rossi, INFN Padova #include "PWGHF/Core/CentralityEstimation.h" #include "PWGHF/Core/HfHelper.h" @@ -251,7 +253,7 @@ struct HfTaskPtFlucCharmHadrons { if constexpr (Channel == DecayChannel::DplusToPiKPi) { std::array daugIDs = {cand.prong0Id(), cand.prong1Id(), cand.prong2Id()}; std::array daugPts = {cand.ptProng0(), cand.ptProng1(), cand.ptProng2()}; - for (int iProng = 0; iProng < 3; ++iProng) { // for 3-prong + for (int iProng = 0; iProng < 3; ++iProng) { // o2-linter: disable=magic-number (for 3-prong) if (std::find(trkIDs.begin(), trkIDs.end(), daugIDs[iProng]) != trkIDs.end()) { meanPt = (rawMeanPt * n - daugPts[iProng]) / (n - 1); } @@ -259,7 +261,7 @@ struct HfTaskPtFlucCharmHadrons { } else if constexpr (Channel == DecayChannel::D0ToPiK || Channel == DecayChannel::D0ToKPi) { std::array daugIDs = {cand.prong0Id(), cand.prong1Id()}; std::array daugPts = {cand.ptProng0(), cand.ptProng1()}; - for (int iProng = 0; iProng < 2; ++iProng) { // for 2-prong + for (int iProng = 0; iProng < 2; ++iProng) { // o2-linter: disable=magic-number (for 2-prong) if (std::find(trkIDs.begin(), trkIDs.end(), daugIDs[iProng]) != trkIDs.end()) { meanPt = (rawMeanPt * n - daugPts[iProng]) / (n - 1); } From 741fa95b933936bc96f40faaffb14caa932a615b Mon Sep 17 00:00:00 2001 From: ALICE Action Bot Date: Wed, 29 Apr 2026 11:15:44 +0000 Subject: [PATCH 07/16] Please consider the following formatting changes --- PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx | 44 +++++++++------------- 1 file changed, 18 insertions(+), 26 deletions(-) diff --git a/PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx b/PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx index 02edb798a8c..cae7bf1d295 100644 --- a/PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx +++ b/PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx @@ -198,26 +198,21 @@ struct HfTaskPtFlucCharmHadrons { template bool isSelectedHfCollision(Coll const& collision, float& cent) { - o2::hf_evsel::HfCollisionRejectionMask collRejMask{}; - if (centEstimator == CentralityEstimator::FT0A) - { - collRejMask = hfEvSel.getHfCollisionRejectionMask(collision, cent, ccdb, registry); - } else if (centEstimator == CentralityEstimator::FT0C) - { - collRejMask = hfEvSel.getHfCollisionRejectionMask(collision, cent, ccdb, registry); - } else if (centEstimator == CentralityEstimator::FT0M) - { - collRejMask = hfEvSel.getHfCollisionRejectionMask(collision, cent, ccdb, registry); - } else if (centEstimator == CentralityEstimator::FV0A) - { - collRejMask = hfEvSel.getHfCollisionRejectionMask(collision, cent, ccdb, registry); - } else - { - LOG(fatal) << "Centrality estimator not recognized for collision selection"; - std::abort(); - } - hfEvSel.fillHistograms(collision, collRejMask, cent); - return collRejMask == 0; + o2::hf_evsel::HfCollisionRejectionMask collRejMask{}; + if (centEstimator == CentralityEstimator::FT0A) { + collRejMask = hfEvSel.getHfCollisionRejectionMask(collision, cent, ccdb, registry); + } else if (centEstimator == CentralityEstimator::FT0C) { + collRejMask = hfEvSel.getHfCollisionRejectionMask(collision, cent, ccdb, registry); + } else if (centEstimator == CentralityEstimator::FT0M) { + collRejMask = hfEvSel.getHfCollisionRejectionMask(collision, cent, ccdb, registry); + } else if (centEstimator == CentralityEstimator::FV0A) { + collRejMask = hfEvSel.getHfCollisionRejectionMask(collision, cent, ccdb, registry); + } else { + LOG(fatal) << "Centrality estimator not recognized for collision selection"; + std::abort(); + } + hfEvSel.fillHistograms(collision, collRejMask, cent); + return collRejMask == 0; } /// Get candidate mass @@ -394,12 +389,10 @@ struct HfTaskPtFlucCharmHadrons { } } - for (const auto& cand : candidates) - { + for (const auto& cand : candidates) { // apply ML selection auto [ml1, ml2] = getMlScores(cand); - if (!passMlCut(ml1, ml2)) - { + if (!passMlCut(ml1, ml2)) { continue; } @@ -426,8 +419,7 @@ struct HfTaskPtFlucCharmHadrons { // fill charm-bulk correlation thnsparse registry.fill(HIST("hCharmBulkCorrelations"), invMass, cent, pt, sign, ml1, ml2, eta, meanPtA, meanPtB, candPtProduct); } - } else - { + } else { int nDcandTotA = 0; for (const auto& cand : candidates) { if (!passCandInA(cand)) { From 694069fd2ac31729fe05fcf09db2868c855979a1 Mon Sep 17 00:00:00 2001 From: wulbyu Date: Wed, 29 Apr 2026 14:59:00 +0200 Subject: [PATCH 08/16] Implemetn comments Co-authored-by: Copilot --- PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx | 22 +++++++++++++++------- 1 file changed, 15 insertions(+), 7 deletions(-) diff --git a/PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx b/PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx index cae7bf1d295..2f6bb4c6b8a 100644 --- a/PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx +++ b/PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx @@ -209,7 +209,6 @@ struct HfTaskPtFlucCharmHadrons { collRejMask = hfEvSel.getHfCollisionRejectionMask(collision, cent, ccdb, registry); } else { LOG(fatal) << "Centrality estimator not recognized for collision selection"; - std::abort(); } hfEvSel.fillHistograms(collision, collRejMask, cent); return collRejMask == 0; @@ -221,10 +220,12 @@ struct HfTaskPtFlucCharmHadrons { { if constexpr (std::is_same_v) { return {HfHelper::invMassDplusToPiKPi(cand), cand.template prong0_as().sign()}; - } else if constexpr (std::is_same_v) { + } + if constexpr (std::is_same_v) { if (channel == DecayChannel::D0ToPiK) { return {HfHelper::invMassD0ToPiK(cand), cand.isSelD0bar() ? 3 : 1}; // 3: reflected D0bar, 1: pure D0 excluding reflected D0bar - } else if (channel == DecayChannel::D0ToKPi) { + } + if (channel == DecayChannel::D0ToKPi) { return {HfHelper::invMassD0barToKPi(cand), cand.isSelD0() ? 3 : 2}; // 3: reflected D0, 2: pure D0bar excluding reflected D0 } } @@ -242,14 +243,14 @@ struct HfTaskPtFlucCharmHadrons { /// remove candidate daughters from the mean pT of tracks in A and B (if they are in the respective subevent) template - float removeDaughtersFromMeanPt(const CandT& cand, const float& rawMeanPt, const int& n, const std::vector& trkIDs) + float removeDaughtersFromMeanPt(const CandT& cand, float rawMeanPt, int n, const std::vector& trkIDs) { float meanPt{0.f}; if constexpr (Channel == DecayChannel::DplusToPiKPi) { std::array daugIDs = {cand.prong0Id(), cand.prong1Id(), cand.prong2Id()}; std::array daugPts = {cand.ptProng0(), cand.ptProng1(), cand.ptProng2()}; for (int iProng = 0; iProng < 3; ++iProng) { // o2-linter: disable=magic-number (for 3-prong) - if (std::find(trkIDs.begin(), trkIDs.end(), daugIDs[iProng]) != trkIDs.end()) { + if (std::binary_search(trkIDs.begin(), trkIDs.end(), daugIDs[iProng])) { meanPt = (rawMeanPt * n - daugPts[iProng]) / (n - 1); } } @@ -257,7 +258,7 @@ struct HfTaskPtFlucCharmHadrons { std::array daugIDs = {cand.prong0Id(), cand.prong1Id()}; std::array daugPts = {cand.ptProng0(), cand.ptProng1()}; for (int iProng = 0; iProng < 2; ++iProng) { // o2-linter: disable=magic-number (for 2-prong) - if (std::find(trkIDs.begin(), trkIDs.end(), daugIDs[iProng]) != trkIDs.end()) { + if (std::binary_search(trkIDs.begin(), trkIDs.end(), daugIDs[iProng])) { meanPt = (rawMeanPt * n - daugPts[iProng]) / (n - 1); } } @@ -376,6 +377,10 @@ struct HfTaskPtFlucCharmHadrons { std::vector trkIDA; std::vector trkIDB; + + trkIDA.reserve(tracks.size() / 2); + trkIDB.reserve(tracks.size() / 2); + // collect track IDs in A and B for (const auto& trk : tracks) { if (!selectionTrack(trk)) { @@ -389,6 +394,9 @@ struct HfTaskPtFlucCharmHadrons { } } + std::sort(trkIDA.begin(), trkIDA.end()); + std::sort(trkIDB.begin(), trkIDB.end()); + for (const auto& cand : candidates) { // apply ML selection auto [ml1, ml2] = getMlScores(cand); @@ -441,7 +449,7 @@ struct HfTaskPtFlucCharmHadrons { ++nDcandTotA; } - if (nDcandTotA <= 0) { + if (nDcandTotA == 0) { return; // cannot build fraction } From 71de245c0b0f24725bee72fd75e87a57b48d9b6b Mon Sep 17 00:00:00 2001 From: ALICE Action Bot Date: Wed, 29 Apr 2026 12:59:54 +0000 Subject: [PATCH 09/16] Please consider the following formatting changes --- PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx b/PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx index 2f6bb4c6b8a..7d5d14a1c1b 100644 --- a/PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx +++ b/PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx @@ -220,7 +220,7 @@ struct HfTaskPtFlucCharmHadrons { { if constexpr (std::is_same_v) { return {HfHelper::invMassDplusToPiKPi(cand), cand.template prong0_as().sign()}; - } + } if constexpr (std::is_same_v) { if (channel == DecayChannel::D0ToPiK) { return {HfHelper::invMassD0ToPiK(cand), cand.isSelD0bar() ? 3 : 1}; // 3: reflected D0bar, 1: pure D0 excluding reflected D0bar @@ -378,7 +378,7 @@ struct HfTaskPtFlucCharmHadrons { std::vector trkIDA; std::vector trkIDB; - trkIDA.reserve(tracks.size() / 2); + trkIDA.reserve(tracks.size() / 2); trkIDB.reserve(tracks.size() / 2); // collect track IDs in A and B From 9aa2f471926d35e1edd65c307b3c3e27eba07ffa Mon Sep 17 00:00:00 2001 From: wulbyu Date: Wed, 29 Apr 2026 15:26:39 +0200 Subject: [PATCH 10/16] head file --- PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx | 1 + 1 file changed, 1 insertion(+) diff --git a/PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx b/PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx index 7d5d14a1c1b..320c43d71aa 100644 --- a/PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx +++ b/PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx @@ -40,6 +40,7 @@ #include #include +#include #include #include #include From 4ad6a73772cf3dc39d071a924ad125161feb66c2 Mon Sep 17 00:00:00 2001 From: ALICE Action Bot Date: Wed, 29 Apr 2026 13:28:23 +0000 Subject: [PATCH 11/16] Please consider the following formatting changes --- PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx b/PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx index 320c43d71aa..38f6e167dad 100644 --- a/PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx +++ b/PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx @@ -39,8 +39,8 @@ #include #include -#include #include +#include #include #include #include From be6729846046c47c55da69d304c5092860a5189e Mon Sep 17 00:00:00 2001 From: wulbyu Date: Wed, 29 Apr 2026 17:39:38 +0200 Subject: [PATCH 12/16] D0 candidate type Co-authored-by: Copilot --- PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx b/PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx index 38f6e167dad..1712c5701ea 100644 --- a/PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx +++ b/PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx @@ -60,6 +60,12 @@ enum DecayChannel { DplusToPiKPi = 0, D0ToPiK, D0ToKPi }; +enum CandD0Type { PureD0 = 0, + PureD0bar, + ReflectedD0, + ReflectedD0bar +}; + struct HfTaskPtFlucCharmHadrons { Configurable centEstimator{"centEstimator", 2, "Centrality estimation (FT0A: 1, FT0C: 2, FT0M: 3, FV0A: 4)"}; Configurable selectionFlag{"selectionFlag", 1, "Selection Flag for hadron (e.g. 1 for skimming, 3 for topo. and kine., 7 for PID)"}; @@ -224,10 +230,10 @@ struct HfTaskPtFlucCharmHadrons { } if constexpr (std::is_same_v) { if (channel == DecayChannel::D0ToPiK) { - return {HfHelper::invMassD0ToPiK(cand), cand.isSelD0bar() ? 3 : 1}; // 3: reflected D0bar, 1: pure D0 excluding reflected D0bar + return {HfHelper::invMassD0ToPiK(cand), cand.isSelD0bar() ? CandD0Type::ReflectedD0bar : CandD0Type::PureD0}; } if (channel == DecayChannel::D0ToKPi) { - return {HfHelper::invMassD0barToKPi(cand), cand.isSelD0() ? 3 : 2}; // 3: reflected D0, 2: pure D0bar excluding reflected D0 + return {HfHelper::invMassD0barToKPi(cand), cand.isSelD0() ? CandD0Type::ReflectedD0 : CandD0Type::PureD0bar}; } } return {0., 0.}; // default return value for unsupported types From f7e83c00104795f55338970ac674cd1c542dfc39 Mon Sep 17 00:00:00 2001 From: wulbyu Date: Thu, 30 Apr 2026 14:15:15 +0200 Subject: [PATCH 13/16] apply the filter on tracks Co-authored-by: Copilot --- PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx b/PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx index 1712c5701ea..8e90aaa5a28 100644 --- a/PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx +++ b/PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx @@ -492,7 +492,7 @@ struct HfTaskPtFlucCharmHadrons { // D0 with ML void processD0Ml(Colls::iterator const& collision, CandD0DataWMl const& /*candidatesD0*/, - TracksWithExtra const& tracks) + soa::Filtered const& tracks) { auto candsD0ToPiKWMl = selectedD0ToPiKWMl->sliceByCached(aod::hf_cand::collisionId, collision.globalIndex(), cache); auto candsD0ToKPiWMl = selectedD0ToKPiWMl->sliceByCached(aod::hf_cand::collisionId, collision.globalIndex(), cache); @@ -504,7 +504,7 @@ struct HfTaskPtFlucCharmHadrons { // Dplus with ML void processDplusMl(Colls::iterator const& collision, CandDplusDataWMl const& candidatesDplus, - TracksWithExtra const& tracks) + soa::Filtered const& tracks) { runPtFlucAnalysis(collision, candidatesDplus, tracks); } From 6f23ff8114779ec67a2f746af204ef68ef93160e Mon Sep 17 00:00:00 2001 From: wulbyu Date: Tue, 5 May 2026 14:01:29 +0200 Subject: [PATCH 14/16] update the algorithm of removing daughter tracks --- PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx | 23 ++++++++++++++-------- 1 file changed, 15 insertions(+), 8 deletions(-) diff --git a/PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx b/PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx index 8e90aaa5a28..3ee5e1f9691 100644 --- a/PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx +++ b/PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx @@ -252,24 +252,31 @@ struct HfTaskPtFlucCharmHadrons { template float removeDaughtersFromMeanPt(const CandT& cand, float rawMeanPt, int n, const std::vector& trkIDs) { - float meanPt{0.f}; + float meanPt = rawMeanPt; + int removedCount = 0; + float removedSumPt = 0.f; + auto removeDaug = [&] (int daugID, float daugPt) { + if (std::binary_search(trkIDs.begin(), trkIDs.end(), daugID)) { + removedSumPt += daugPt; + ++removedCount; + } + }; if constexpr (Channel == DecayChannel::DplusToPiKPi) { std::array daugIDs = {cand.prong0Id(), cand.prong1Id(), cand.prong2Id()}; std::array daugPts = {cand.ptProng0(), cand.ptProng1(), cand.ptProng2()}; - for (int iProng = 0; iProng < 3; ++iProng) { // o2-linter: disable=magic-number (for 3-prong) - if (std::binary_search(trkIDs.begin(), trkIDs.end(), daugIDs[iProng])) { - meanPt = (rawMeanPt * n - daugPts[iProng]) / (n - 1); - } + for (int iProng = 0; iProng < 3; ++iProng) { + removeDaug(daugIDs[iProng], daugPts[iProng]); } } else if constexpr (Channel == DecayChannel::D0ToPiK || Channel == DecayChannel::D0ToKPi) { std::array daugIDs = {cand.prong0Id(), cand.prong1Id()}; std::array daugPts = {cand.ptProng0(), cand.ptProng1()}; for (int iProng = 0; iProng < 2; ++iProng) { // o2-linter: disable=magic-number (for 2-prong) - if (std::binary_search(trkIDs.begin(), trkIDs.end(), daugIDs[iProng])) { - meanPt = (rawMeanPt * n - daugPts[iProng]) / (n - 1); - } + removeDaug(daugIDs[iProng], daugPts[iProng]); } } + if (removedCount > 0) { + meanPt = (rawMeanPt * n - removedSumPt) / (n - removedCount); + } return meanPt; } From f6880e660bd9c4f15160b6c8f6549e481d3fb253 Mon Sep 17 00:00:00 2001 From: wulbyu Date: Tue, 5 May 2026 14:11:42 +0200 Subject: [PATCH 15/16] refactor the application --- PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx | 31 +++++++++------------- 1 file changed, 13 insertions(+), 18 deletions(-) diff --git a/PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx b/PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx index 3ee5e1f9691..e2109e3bc35 100644 --- a/PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx +++ b/PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx @@ -252,32 +252,27 @@ struct HfTaskPtFlucCharmHadrons { template float removeDaughtersFromMeanPt(const CandT& cand, float rawMeanPt, int n, const std::vector& trkIDs) { - float meanPt = rawMeanPt; int removedCount = 0; float removedSumPt = 0.f; - auto removeDaug = [&] (int daugID, float daugPt) { - if (std::binary_search(trkIDs.begin(), trkIDs.end(), daugID)) { - removedSumPt += daugPt; - ++removedCount; - } + auto removeDaug = [&](int daugID, float daugPt) { + if (std::binary_search(trkIDs.begin(), trkIDs.end(), daugID)) { + removedSumPt += daugPt; + ++removedCount; + } }; if constexpr (Channel == DecayChannel::DplusToPiKPi) { - std::array daugIDs = {cand.prong0Id(), cand.prong1Id(), cand.prong2Id()}; - std::array daugPts = {cand.ptProng0(), cand.ptProng1(), cand.ptProng2()}; - for (int iProng = 0; iProng < 3; ++iProng) { - removeDaug(daugIDs[iProng], daugPts[iProng]); - } + removeDaug(cand.prong0Id(), cand.ptProng0()); + removeDaug(cand.prong1Id(), cand.ptProng1()); + removeDaug(cand.prong2Id(), cand.ptProng2()); } else if constexpr (Channel == DecayChannel::D0ToPiK || Channel == DecayChannel::D0ToKPi) { - std::array daugIDs = {cand.prong0Id(), cand.prong1Id()}; - std::array daugPts = {cand.ptProng0(), cand.ptProng1()}; - for (int iProng = 0; iProng < 2; ++iProng) { // o2-linter: disable=magic-number (for 2-prong) - removeDaug(daugIDs[iProng], daugPts[iProng]); - } + removeDaug(cand.prong0Id(), cand.ptProng0()); + removeDaug(cand.prong1Id(), cand.ptProng1()); } if (removedCount > 0) { - meanPt = (rawMeanPt * n - removedSumPt) / (n - removedCount); + double totalSum = static_cast(rawMeanPt) * n; + return static_cast((totalSum - removedSumPt) / (n - removedCount)); } - return meanPt; + return rawMeanPt; } // --------------------------------------------------------------------------- From 1c430ce22d382c2106c7ad09ad7bd175e69b6f28 Mon Sep 17 00:00:00 2001 From: ALICE Action Bot Date: Tue, 5 May 2026 12:12:42 +0000 Subject: [PATCH 16/16] Please consider the following formatting changes --- PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx b/PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx index e2109e3bc35..f705e0ea31d 100644 --- a/PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx +++ b/PWGHF/D2H/Tasks/taskPtFlucCharmHadrons.cxx @@ -255,18 +255,18 @@ struct HfTaskPtFlucCharmHadrons { int removedCount = 0; float removedSumPt = 0.f; auto removeDaug = [&](int daugID, float daugPt) { - if (std::binary_search(trkIDs.begin(), trkIDs.end(), daugID)) { - removedSumPt += daugPt; - ++removedCount; - } + if (std::binary_search(trkIDs.begin(), trkIDs.end(), daugID)) { + removedSumPt += daugPt; + ++removedCount; + } }; if constexpr (Channel == DecayChannel::DplusToPiKPi) { - removeDaug(cand.prong0Id(), cand.ptProng0()); - removeDaug(cand.prong1Id(), cand.ptProng1()); - removeDaug(cand.prong2Id(), cand.ptProng2()); + removeDaug(cand.prong0Id(), cand.ptProng0()); + removeDaug(cand.prong1Id(), cand.ptProng1()); + removeDaug(cand.prong2Id(), cand.ptProng2()); } else if constexpr (Channel == DecayChannel::D0ToPiK || Channel == DecayChannel::D0ToKPi) { - removeDaug(cand.prong0Id(), cand.ptProng0()); - removeDaug(cand.prong1Id(), cand.ptProng1()); + removeDaug(cand.prong0Id(), cand.ptProng0()); + removeDaug(cand.prong1Id(), cand.ptProng1()); } if (removedCount > 0) { double totalSum = static_cast(rawMeanPt) * n;