From 80fae3ec0cd77d8686711aca8b7a1f92799c2b4b Mon Sep 17 00:00:00 2001 From: Marcello Di Costanzo <96481191+Marcellocosti@users.noreply.github.com> Date: Sat, 11 Jan 2025 15:07:48 +0100 Subject: [PATCH] [PWGHF] Extend information of sparses for efficiency evaluation. Refactor analysis utilities. (#9153) Co-authored-by: ALICE Action Bot --- PWGHF/Core/CentralityEstimation.h | 76 ++++ PWGHF/D2H/Tasks/taskDplus.cxx | 131 +------ PWGHF/D2H/Tasks/taskDs.cxx | 351 ++++++++++++------ PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx | 37 +- .../TableProducer/candidateCreator2Prong.cxx | 7 +- .../TableProducer/candidateCreator3Prong.cxx | 7 +- PWGHF/Tasks/taskMcValidation.cxx | 24 +- PWGHF/Utils/utilsAnalysis.h | 27 ++ PWGHF/Utils/utilsEvSelHf.h | 111 +++--- 9 files changed, 421 insertions(+), 350 deletions(-) diff --git a/PWGHF/Core/CentralityEstimation.h b/PWGHF/Core/CentralityEstimation.h index ea9d1f257b6..203cc4680c8 100644 --- a/PWGHF/Core/CentralityEstimation.h +++ b/PWGHF/Core/CentralityEstimation.h @@ -113,6 +113,82 @@ float getCentralityColl(const Coll&) return 105.0f; } +/// Get the centrality +/// \param collision is the collision with the centrality information +/// \param centEstimator integer to select the centrality estimator +/// \return collision centrality +template +float getCentralityColl(const Coll& collision, int centEstimator) +{ + switch (centEstimator) { + case CentralityEstimator::FT0A: + if constexpr (hasFT0ACent) { + return collision.centFT0A(); + } + LOG(fatal) << "Collision does not have centFT0A()."; + break; + case CentralityEstimator::FT0C: + if constexpr (hasFT0CCent) { + return collision.centFT0C(); + } + LOG(fatal) << "Collision does not have centFT0C()."; + break; + case CentralityEstimator::FT0M: + if constexpr (hasFT0MCent) { + return collision.centFT0M(); + } + LOG(fatal) << "Collision does not have centFT0M()."; + break; + case CentralityEstimator::FV0A: + if constexpr (hasFV0ACent) { + return collision.centFV0A(); + } + LOG(fatal) << "Collision does not have centFV0A()."; + break; + default: + LOG(fatal) << "Centrality estimator not valid. Possible values are V0A, T0M, T0A, T0C."; + break; + } + return -999.f; +} + +/// \brief Function to get MC collision centrality +/// \param collSlice collection of reconstructed collisions associated to a generated one +/// \return generated MC collision centrality +template +float getCentralityGenColl(CCs const& collSlice) +{ + float centrality{-1}; + float multiplicity{0.f}; + for (const auto& collision : collSlice) { + float collMult = collision.numContrib(); + if (collMult > multiplicity) { + centrality = getCentralityColl(collision); + multiplicity = collMult; + } + } + return centrality; +} + +/// \brief Function to get MC collision centrality +/// \param collSlice collection of reconstructed collisions associated to a generated one +/// \param centEstimator integer to select the centrality estimator +/// \return generated MC collision centrality +template +float getCentralityGenColl(CCs const& collSlice, int centEstimator) +{ + float centrality{-1}; + float multiplicity{0.f}; + for (const auto& collision : collSlice) { + float collMult = collision.numContrib(); + if (collMult > multiplicity) { + centrality = getCentralityColl(collision, centEstimator); + multiplicity = collMult; + } + } + return centrality; +} + } // namespace o2::hf_centrality #endif // PWGHF_CORE_CENTRALITYESTIMATION_H_ diff --git a/PWGHF/D2H/Tasks/taskDplus.cxx b/PWGHF/D2H/Tasks/taskDplus.cxx index 625e12581cc..9026e4c97d6 100644 --- a/PWGHF/D2H/Tasks/taskDplus.cxx +++ b/PWGHF/D2H/Tasks/taskDplus.cxx @@ -28,22 +28,15 @@ #include "PWGHF/Core/HfHelper.h" #include "PWGHF/DataModel/CandidateReconstructionTables.h" #include "PWGHF/DataModel/CandidateSelectionTables.h" +#include "PWGHF/Utils/utilsEvSelHf.h" +#include "PWGHF/Utils/utilsAnalysis.h" using namespace o2; using namespace o2::analysis; using namespace o2::framework; using namespace o2::framework::expressions; using namespace o2::hf_centrality; - -enum OccupancyEstimator { None = 0, - ITS, - FT0C }; - -enum BHadMothers { NotMatched = 0, - BPlus, - BZero, - Bs, - LambdaBZero }; +using namespace o2::hf_occupancy; /// D± analysis task struct HfTaskDplus { @@ -455,10 +448,10 @@ struct HfTaskDplus { if (storeCentrality || storeOccupancy) { auto collision = candidate.template collision_as(); if (storeCentrality && centEstimator != CentralityEstimator::None) { - cent = getCentrality(collision); + cent = getCentralityColl(collision, centEstimator); } if (storeOccupancy && occEstimator != OccupancyEstimator::None) { - occ = getOccupancy(collision); + occ = getOccupancyColl(collision, occEstimator); } } @@ -506,10 +499,10 @@ struct HfTaskDplus { if (storeCentrality || storeOccupancy) { auto collision = candidate.template collision_as(); if (storeCentrality && centEstimator != CentralityEstimator::None) { - cent = getCentrality(collision); + cent = getCentralityColl(collision, centEstimator); } if (storeOccupancy && occEstimator != OccupancyEstimator::None) { - occ = getOccupancy(collision); + occ = getOccupancyColl(collision, occEstimator); } } @@ -526,10 +519,10 @@ struct HfTaskDplus { } auto collision = candidate.template collision_as(); if (storeCentrality && centEstimator != CentralityEstimator::None) { - cent = getCentrality(collision); + cent = getCentralityColl(collision, centEstimator); } if (storeOccupancy && occEstimator != OccupancyEstimator::None) { - occ = getOccupancy(collision); + occ = getOccupancyColl(collision, occEstimator); } fillHistoMCRec(candidate); fillSparseML(candidate, ptBhad, flagBHad, cent, occ); @@ -556,10 +549,10 @@ struct HfTaskDplus { const auto recoCollsPerGenMcColl = mcRecoCollisions.sliceBy(recoColPerMcCollision, mcGenCollision.globalIndex()); const auto mcParticlesPerGenMcColl = mcGenParticles.sliceBy(mcParticlesPerMcCollision, mcGenCollision.globalIndex()); if (storeCentrality && centEstimator != CentralityEstimator::None) { - cent = getMcGenCollCentrality(recoCollsPerGenMcColl); + cent = getCentralityGenColl(recoCollsPerGenMcColl, centEstimator); } if (storeOccupancy && occEstimator != OccupancyEstimator::None) { - occ = getMcGenCollOccupancy(recoCollsPerGenMcColl); + occ = getOccupancyGenColl(recoCollsPerGenMcColl, occEstimator); } for (const auto& particle : mcParticlesPerGenMcColl) { @@ -582,108 +575,6 @@ struct HfTaskDplus { } } - /// Get the occupancy - /// \param collision is the collision with the occupancy information - /// \return collision occupancy - template - float getOccupancy(Coll const& collision) - { - float occupancy = -999.; - switch (occEstimator) { - case OccupancyEstimator::ITS: - occupancy = collision.trackOccupancyInTimeRange(); - break; - case OccupancyEstimator::FT0C: - occupancy = collision.ft0cOccupancyInTimeRange(); - break; - default: - LOG(warning) << "Occupancy estimator not valid. Possible values are ITS or FT0C. Fallback to ITS"; - occupancy = collision.trackOccupancyInTimeRange(); - break; - } - return occupancy; - } - - /// \brief Function to get MC collision occupancy - /// \param collSlice collection of reconstructed collisions associated to a generated one - /// \return generated MC collision occupancy - template - int getMcGenCollOccupancy(CCs const& collSlice) - { - float multiplicity{0.f}; - int occupancy = 0; - for (const auto& collision : collSlice) { - float collMult{0.f}; - collMult = collision.numContrib(); - if (collMult > multiplicity) { - occupancy = getOccupancy(collision); - multiplicity = collMult; - } - } // end loop over collisions - - return occupancy; - } - - /// Get the centrality - /// \param collision is the collision with the centrality information - /// \return collision centrality - template - float getCentrality(Coll const& collision) - { - float cent = -999.; - switch (centEstimator) { - case CentralityEstimator::FT0C: - cent = collision.centFT0C(); - break; - case CentralityEstimator::FT0M: - cent = collision.centFT0M(); - break; - default: - LOG(warning) << "Centrality estimator not valid. Possible values are FT0C, FT0M. Fallback to FT0C"; - cent = collision.centFT0C(); - break; - } - return cent; - } - - /// \brief Function to get MC collision centrality - /// \param collSlice collection of reconstructed collisions associated to a generated one - /// \return generated MC collision centrality - template - float getMcGenCollCentrality(CCs const& collSlice) - { - float centrality{-1}; - float multiplicity{0.f}; - for (const auto& collision : collSlice) { - float collMult = collision.numContrib(); - if (collMult > multiplicity) { - centrality = getCentrality(collision); - multiplicity = collMult; - } - } - return centrality; - } - - /// Convert the B hadron mother PDG for non prompt candidates to a flag - /// \param pdg of the b hadron mother - /// \return integer map to specific mothers' PDG codes - int getBHadMotherFlag(const int& flagBHad) - { - if (std::abs(flagBHad) == o2::constants::physics::kBPlus) { - return BHadMothers::BPlus; - } - if (std::abs(flagBHad) == o2::constants::physics::kB0) { - return BHadMothers::BZero; - } - if (std::abs(flagBHad) == o2::constants::physics::kBS) { - return BHadMothers::Bs; - } - if (std::abs(flagBHad) == o2::constants::physics::kLambdaB0) { - return BHadMothers::LambdaBZero; - } - return BHadMothers::NotMatched; - } - // process functions void processData(CandDplusData const& candidates, CollisionsCent const& collisions) { diff --git a/PWGHF/D2H/Tasks/taskDs.cxx b/PWGHF/D2H/Tasks/taskDs.cxx index e9c49d9fc86..4c9a966e6b9 100644 --- a/PWGHF/D2H/Tasks/taskDs.cxx +++ b/PWGHF/D2H/Tasks/taskDs.cxx @@ -34,6 +34,8 @@ #include "PWGHF/Core/CentralityEstimation.h" #include "PWGHF/DataModel/CandidateReconstructionTables.h" #include "PWGHF/DataModel/CandidateSelectionTables.h" +#include "PWGHF/Utils/utilsEvSelHf.h" +#include "PWGHF/Utils/utilsAnalysis.h" using namespace o2; using namespace o2::analysis; @@ -76,27 +78,23 @@ struct HfTaskDs { Configurable massDplusSignalMin{"massDplusSignalMin", 1.866, "min mass for Dplus signal"}; Configurable massDplusSignalMax{"massDplusSignalMax", 1.906, "max mass for Dplus signal"}; Configurable fillPercentiles{"fillPercentiles", true, "Wheter to fill multiplicity axis with percentiles or raw information"}; - - ConfigurableAxis axisPt{"axisPt", {VARIABLE_WIDTH, 1.f, 2.f, 3.f, 4.f, 5.f, 6.f, 8.f, 12.f, 24.f}, "axis for pT"}; - ConfigurableAxis axisNPvContributors{"axisNPvContributors", {200, -0.5f, 199.5f}, "axis for NPvContributors"}; - ConfigurableAxis axisMlScore0{"axisMlScore0", {100, 0., 1.}, "axis for ML output score 0"}; - ConfigurableAxis axisMlScore1{"axisMlScore1", {100, 0., 1.}, "axis for ML output score 1"}; - ConfigurableAxis axisMlScore2{"axisMlScore2", {100, 0., 1.}, "axis for ML output score 2"}; - ConfigurableAxis axisCentrality{"axisCentrality", {100, 0., 1.}, "axis for centrality/multiplicity"}; + Configurable storeOccupancy{"storeOccupancy", false, "Flag to store occupancy information"}; + Configurable occEstimator{"occEstimator", 0, "Occupancy estimation (None: 0, ITS: 1, FT0C: 2)"}; struct : ConfigurableGroup { Configurable ccdburl{"ccdburl", "http://alice-ccdb.cern.ch", "The CCDB endpoint url address"}; - Configurable ccdbPath{"ccdbpath", "Centrality/Calibration", "The CCDB path for centrality/multiplicity information"}; + Configurable ccdbPath{"ccdbPath", "Centrality/Calibration", "The CCDB path for centrality/multiplicity information"}; Configurable reconstructionPass{"reconstructionPass", "", {"Apass to use when fetching the calibration tables. Empty (default) does not check for any pass. Use `metadata` to fetch it from the AO2D metadata. Otherwise it will override the metadata."}}; } ccdbConfig; HfHelper hfHelper; + SliceCache cache; Service ccdb; - using TH1_ptr = std::shared_ptr; - using TH2_ptr = std::shared_ptr; - using THnSparse_ptr = std::shared_ptr; - using histTypes = std::variant; + using TH1Ptr = std::shared_ptr; + using TH2Ptr = std::shared_ptr; + using THnSparsePtr = std::shared_ptr; + using HistTypes = std::variant; template using MemberFunctionPointer = bool (HfTaskDs::*)(const CandDs&); @@ -115,9 +113,20 @@ struct HfTaskDs { using CandDsMcRecoWithMl = soa::Filtered>; using CandDsMcGen = soa::Join; + Filter filterDsFlag = (o2::aod::hf_track_index::hfflag & static_cast(BIT(aod::hf_cand_3prong::DecayType::DsToKKPi))) != static_cast(0); + Preslice candDsPerCollision = aod::hf_cand::collisionId; PresliceUnsorted colPerMcCollision = aod::mccollisionlabel::mcCollisionId; - SliceCache cache; + + ConfigurableAxis axisPt{"axisPt", {VARIABLE_WIDTH, 1.f, 2.f, 3.f, 4.f, 5.f, 6.f, 8.f, 12.f, 24.f}, "axis for pT"}; + ConfigurableAxis axisPtBHad{"axisPtBHad", {50, 0., 100}, "axis for pt of B hadron decayed into D candidate"}; + ConfigurableAxis axisFlagBHad{"axisFlagBHad", {5, 0, 5}, "axis for B hadron mother flag"}; + ConfigurableAxis axisNPvContributors{"axisNPvContributors", {200, -0.5f, 199.5f}, "axis for NPvContributors"}; + ConfigurableAxis axisMlScore0{"axisMlScore0", {100, 0., 1.}, "axis for ML output score 0"}; + ConfigurableAxis axisMlScore1{"axisMlScore1", {100, 0., 1.}, "axis for ML output score 1"}; + ConfigurableAxis axisMlScore2{"axisMlScore2", {100, 0., 1.}, "axis for ML output score 2"}; + ConfigurableAxis axisCentrality{"axisCentrality", {100, 0., 1.}, "axis for centrality/multiplicity"}; + ConfigurableAxis axisOccupancy{"axisOccupancy", {14, 0., 14000.}, "axis for occupancy"}; int offsetDplusDecayChannel = aod::hf_cand_3prong::DecayChannelDToKKPi::DplusToPhiPi - aod::hf_cand_3prong::DecayChannelDToKKPi::DsToPhiPi; // Offset between Dplus and Ds to use the same decay channel. See aod::hf_cand_3prong::DecayChannelDToKKPi int mRunNumber{0}; @@ -127,22 +136,20 @@ struct HfTaskDs { TProfile* hVtxZFT0C; TProfile* hVtxZNTracks; - Filter filterDsFlag = (o2::aod::hf_track_index::hfflag & static_cast(BIT(aod::hf_cand_3prong::DecayType::DsToKKPi))) != static_cast(0); - HistogramRegistry registry{"registry", {}}; std::array folders = {"Data/", "MC/Ds/Prompt/", "MC/Ds/NonPrompt/", "MC/Dplus/Prompt/", "MC/Dplus/NonPrompt/", "MC/Dplus/Bkg/", "MC/Lc/", "MC/Bkg/"}; - std::unordered_map dataHistograms = {}; - std::unordered_map mcDsPromptHistograms = {}; - std::unordered_map mcDsNonPromptHistograms = {}; - std::unordered_map mcDplusPromptHistograms = {}; - std::unordered_map mcDplusNonPromptHistograms = {}; - std::unordered_map mcDplusBkgHistograms = {}; - std::unordered_map mcLcBkgHistograms = {}; - std::unordered_map mcBkgHistograms = {}; + std::unordered_map dataHistograms = {}; + std::unordered_map mcDsPromptHistograms = {}; + std::unordered_map mcDsNonPromptHistograms = {}; + std::unordered_map mcDplusPromptHistograms = {}; + std::unordered_map mcDplusNonPromptHistograms = {}; + std::unordered_map mcDplusBkgHistograms = {}; + std::unordered_map mcLcBkgHistograms = {}; + std::unordered_map mcBkgHistograms = {}; - std::array, DataType::kDataTypes> histosPtr = {dataHistograms, mcDsPromptHistograms, mcDsNonPromptHistograms, mcDplusPromptHistograms, mcDplusNonPromptHistograms, mcDplusBkgHistograms, mcLcBkgHistograms, mcBkgHistograms}; + std::array, DataType::kDataTypes> histosPtr = {dataHistograms, mcDsPromptHistograms, mcDsNonPromptHistograms, mcDplusPromptHistograms, mcDplusNonPromptHistograms, mcDplusBkgHistograms, mcLcBkgHistograms, mcBkgHistograms}; void init(InitContext&) { @@ -156,24 +163,57 @@ struct HfTaskDs { } AxisSpec ptbins{axisPt, "#it{p}_{T} (GeV/#it{c})"}; + AxisSpec ptBHad{axisPtBHad, "#it{p}_{T}(B) (GeV/#it{c})"}; + AxisSpec flagBHad{axisFlagBHad, "B Hadron flag"}; AxisSpec ybins = {100, -5., 5, "#it{y}"}; AxisSpec massbins = {600, 1.67, 2.27, "inv. mass (KK#pi) (GeV/#it{c}^{2})"}; AxisSpec centralitybins = {100, 0., 100., "Centrality"}; + AxisSpec npvcontributorsbins = {axisNPvContributors, "NPvContributors"}; + AxisSpec mlscore0bins = {axisMlScore0, "Score 0"}; + AxisSpec mlscore1bins = {axisMlScore1, "Score 1"}; + AxisSpec mlscore2bins = {axisMlScore2, "Score 2"}; + AxisSpec occupancybins = {axisOccupancy, "Occupancy"}; histosPtr[DataType::Data]["hNPvContribAll"] = registry.add((folders[DataType::Data] + "hNPvContribAll").c_str(), "3-prong candidates;NPvContributors;Centrality;Entries", HistType::kTH2F, {axisNPvContributors, {100, 0., 100}}); + std::vector axes = {massbins, ptbins, centralitybins}; + std::vector axesMl = {massbins, ptbins, centralitybins, mlscore0bins, mlscore1bins, mlscore2bins}; + std::vector axesFd = {massbins, ptbins, centralitybins, ptBHad, flagBHad}; + std::vector axesFdMl = {massbins, ptbins, centralitybins, mlscore0bins, mlscore1bins, mlscore2bins, ptBHad, flagBHad}; + std::vector axesWithNpv = {massbins, ptbins, centralitybins, npvcontributorsbins}; + std::vector axesWithNpvMl = {massbins, ptbins, centralitybins, npvcontributorsbins, mlscore0bins, mlscore1bins, mlscore2bins}; + std::vector axesGenPrompt = {ptbins, ybins, npvcontributorsbins, centralitybins}; + std::vector axesGenFd = {ptbins, ybins, npvcontributorsbins, ptBHad, flagBHad, centralitybins}; + std::vector axesGenBkg = {ptbins, ybins, npvcontributorsbins, centralitybins}; + + if (storeOccupancy) { + axes.insert(axes.end(), {occupancybins}); + axesMl.insert(axesMl.end(), {occupancybins}); + axesFd.insert(axesFd.end(), {occupancybins}); + axesFdMl.insert(axesFdMl.end(), {occupancybins}); + axesWithNpv.insert(axesWithNpv.end(), {occupancybins}); + axesWithNpvMl.insert(axesWithNpvMl.end(), {occupancybins}); + axesGenPrompt.insert(axesGenPrompt.end(), {occupancybins}); + axesGenFd.insert(axesGenFd.end(), {occupancybins}); + axesGenBkg.insert(axesGenBkg.end(), {occupancybins}); + } + for (auto i = 0; i < DataType::kDataTypes; ++i) { if (doprocessDataWithCentFT0C || doprocessDataWithCentFT0M || doprocessDataWithCentNTracksPV || doprocessData || doprocessMcWithCentFT0C || doprocessMcWithCentFT0M || doprocessMcWithCentNTracksPV || doprocessMc) { if (i == DataType::Data) { // If data do not fill PV contributors in sparse - histosPtr[i]["hSparseMass"] = registry.add((folders[i] + "hSparseMass").c_str(), "THn for Ds", HistType::kTHnSparseF, {massbins, ptbins, centralitybins}); + histosPtr[i]["hSparseMass"] = registry.add((folders[i] + "hSparseMass").c_str(), "THn for Ds", HistType::kTHnSparseF, axes); + } else if (i == DataType::McDsNonPrompt) { // If data do not fill PV contributors in sparse + histosPtr[i]["hSparseMass"] = registry.add((folders[i] + "hSparseMass").c_str(), "THn for Ds", HistType::kTHnSparseF, axesFd); } else { - histosPtr[i]["hSparseMass"] = registry.add((folders[i] + "hSparseMass").c_str(), "THn for Ds", HistType::kTHnSparseF, {massbins, ptbins, centralitybins, axisNPvContributors}); + histosPtr[i]["hSparseMass"] = registry.add((folders[i] + "hSparseMass").c_str(), "THn for Ds", HistType::kTHnSparseF, axesWithNpv); } } else if (doprocessDataWithMlAndCentFT0C || doprocessDataWithMlAndCentFT0M || doprocessDataWithMlAndCentNTracksPV || doprocessDataWithMl || doprocessMcWithMlAndCentFT0C || doprocessMcWithMlAndCentFT0M || doprocessMcWithMlAndCentNTracksPV || doprocessMcWithMl) { if (i == DataType::Data) { // If data do not fill PV contributors in sparse - histosPtr[i]["hSparseMass"] = registry.add((folders[i] + "hSparseMass").c_str(), "THn for Ds", HistType::kTHnSparseF, {massbins, ptbins, centralitybins, axisMlScore0, axisMlScore1, axisMlScore2}); + histosPtr[i]["hSparseMass"] = registry.add((folders[i] + "hSparseMass").c_str(), "THn for Ds", HistType::kTHnSparseF, axesMl); + } else if (i == DataType::McDsNonPrompt) { // If data do not fill PV contributors in sparse + histosPtr[i]["hSparseMass"] = registry.add((folders[i] + "hSparseMass").c_str(), "THn for Ds", HistType::kTHnSparseF, axesFdMl); } else { - histosPtr[i]["hSparseMass"] = registry.add((folders[i] + "hSparseMass").c_str(), "THn for Ds", HistType::kTHnSparseF, {massbins, ptbins, centralitybins, axisMlScore0, axisMlScore1, axisMlScore2, axisNPvContributors}); + histosPtr[i]["hSparseMass"] = registry.add((folders[i] + "hSparseMass").c_str(), "THn for Ds", HistType::kTHnSparseF, axesWithNpvMl); } } histosPtr[i]["hNPvContribCands"] = registry.add((folders[i] + "hNPvContribCands").c_str(), "3-prong candidates;NPvContributors;Centrality;Entries", HistType::kTH2F, {axisNPvContributors, centralitybins}); @@ -211,13 +251,20 @@ struct HfTaskDs { for (auto i = 0; i < DataType::kDataTypes; ++i) { if (i == DataType::McDsPrompt || i == DataType::McDsNonPrompt || i == DataType::McDplusPrompt || i == DataType::McDplusNonPrompt || i == DataType::McDplusBkg || i == DataType::McLcBkg) { - histosPtr[i]["hEtaGen"] = registry.add((folders[i] + "hEtaGen").c_str(), "3-prong candidates (matched);#eta;entries", {HistType::kTH1F, {{100, -2., 2.}}}); histosPtr[i]["hPtGen"] = registry.add((folders[i] + "hPtGen").c_str(), "MC particles (unmatched);#it{p}_{T}^{gen.} (GeV/#it{c});entries", {HistType::kTH1F, {ptbins}}); histosPtr[i]["hPtVsYRecoPID"] = registry.add((folders[i] + "hPtVsYRecoPID").c_str(), "3-prong candidates (RecoPID - matched);#it{p}_{T}^{rec.}; #it{y}", {HistType::kTH2F, {ptbins, {ybins}}}); histosPtr[i]["hPtVsYRecoTopol"] = registry.add((folders[i] + "hPtVsYRecoTopol").c_str(), "3-prong candidates (RecoTopol - matched);#it{p}_{T}^{rec.}; #it{y}", {HistType::kTH2F, {ptbins, {ybins}}}); histosPtr[i]["hPtVsYRecoSkim"] = registry.add((folders[i] + "hPtVsYRecoSkim").c_str(), "3-prong candidates (RecoSkim - matched);#it{p}_{T}^{rec.}; #it{y}", {HistType::kTH2F, {ptbins, {ybins}}}); - histosPtr[i]["hPtYNPvContribGen"] = registry.add((folders[i] + "hPtYNPvContribGen").c_str(), "Thn for generated candidates", {HistType::kTHnSparseF, {ptbins, {ybins}, axisNPvContributors}}); + } + if (i == DataType::McDsPrompt || i == DataType::McDplusPrompt) { + histosPtr[i]["hSparseGen"] = registry.add((folders[i] + "hSparseGen").c_str(), "Thn for generated prompt candidates", HistType::kTHnSparseF, axesGenPrompt); + } + if (i == DataType::McDsNonPrompt || i == DataType::McDplusNonPrompt) { + histosPtr[i]["hSparseGen"] = registry.add((folders[i] + "hSparseGen").c_str(), "Thn for generated nonprompt candidates", HistType::kTHnSparseF, axesGenFd); + } + if (i == DataType::McBkg) { + histosPtr[i]["hSparseGen"] = registry.add((folders[i] + "hSparseGen").c_str(), "Thn for non-matched generated candidates", HistType::kTHnSparseF, axesGenBkg); } } } @@ -352,28 +399,28 @@ struct HfTaskDs { void fillHisto(const T1& candidate, DataType dataType) { auto pt = candidate.pt(); - std::get(histosPtr[dataType]["hPt"])->Fill(pt); - std::get(histosPtr[dataType]["hPtProng0"])->Fill(candidate.ptProng0()); - std::get(histosPtr[dataType]["hPtProng1"])->Fill(candidate.ptProng1()); - std::get(histosPtr[dataType]["hPtProng2"])->Fill(candidate.ptProng2()); - std::get(histosPtr[dataType]["hEta"])->Fill(candidate.eta(), pt); - std::get(histosPtr[dataType]["hCt"])->Fill(hfHelper.ctDs(candidate), pt); - std::get(histosPtr[dataType]["hDecayLength"])->Fill(candidate.decayLength(), pt); - std::get(histosPtr[dataType]["hDecayLengthXY"])->Fill(candidate.decayLengthXY(), pt); - std::get(histosPtr[dataType]["hNormalisedDecayLengthXY"])->Fill(candidate.decayLengthXYNormalised(), pt); - std::get(histosPtr[dataType]["hCPA"])->Fill(candidate.cpa(), pt); - std::get(histosPtr[dataType]["hCPAxy"])->Fill(candidate.cpaXY(), pt); - std::get(histosPtr[dataType]["hImpactParameterXY"])->Fill(candidate.impactParameterXY(), pt); - std::get(histosPtr[dataType]["hMaxNormalisedDeltaIP"])->Fill(candidate.maxNormalisedDeltaIP(), pt); - std::get(histosPtr[dataType]["hImpactParameterProngSqSum"])->Fill(candidate.impactParameterProngSqSum(), pt); - std::get(histosPtr[dataType]["hDecayLengthError"])->Fill(candidate.errorDecayLength(), pt); - std::get(histosPtr[dataType]["hDecayLengthXYError"])->Fill(candidate.errorDecayLengthXY(), pt); - std::get(histosPtr[dataType]["hImpactParameterError"])->Fill(candidate.errorImpactParameter0(), pt); - std::get(histosPtr[dataType]["hImpactParameterError"])->Fill(candidate.errorImpactParameter1(), pt); - std::get(histosPtr[dataType]["hImpactParameterError"])->Fill(candidate.errorImpactParameter2(), pt); - std::get(histosPtr[dataType]["hd0Prong0"])->Fill(candidate.impactParameter0(), pt); - std::get(histosPtr[dataType]["hd0Prong1"])->Fill(candidate.impactParameter1(), pt); - std::get(histosPtr[dataType]["hd0Prong2"])->Fill(candidate.impactParameter2(), pt); + std::get(histosPtr[dataType]["hPt"])->Fill(pt); + std::get(histosPtr[dataType]["hPtProng0"])->Fill(candidate.ptProng0()); + std::get(histosPtr[dataType]["hPtProng1"])->Fill(candidate.ptProng1()); + std::get(histosPtr[dataType]["hPtProng2"])->Fill(candidate.ptProng2()); + std::get(histosPtr[dataType]["hEta"])->Fill(candidate.eta(), pt); + std::get(histosPtr[dataType]["hCt"])->Fill(hfHelper.ctDs(candidate), pt); + std::get(histosPtr[dataType]["hDecayLength"])->Fill(candidate.decayLength(), pt); + std::get(histosPtr[dataType]["hDecayLengthXY"])->Fill(candidate.decayLengthXY(), pt); + std::get(histosPtr[dataType]["hNormalisedDecayLengthXY"])->Fill(candidate.decayLengthXYNormalised(), pt); + std::get(histosPtr[dataType]["hCPA"])->Fill(candidate.cpa(), pt); + std::get(histosPtr[dataType]["hCPAxy"])->Fill(candidate.cpaXY(), pt); + std::get(histosPtr[dataType]["hImpactParameterXY"])->Fill(candidate.impactParameterXY(), pt); + std::get(histosPtr[dataType]["hMaxNormalisedDeltaIP"])->Fill(candidate.maxNormalisedDeltaIP(), pt); + std::get(histosPtr[dataType]["hImpactParameterProngSqSum"])->Fill(candidate.impactParameterProngSqSum(), pt); + std::get(histosPtr[dataType]["hDecayLengthError"])->Fill(candidate.errorDecayLength(), pt); + std::get(histosPtr[dataType]["hDecayLengthXYError"])->Fill(candidate.errorDecayLengthXY(), pt); + std::get(histosPtr[dataType]["hImpactParameterError"])->Fill(candidate.errorImpactParameter0(), pt); + std::get(histosPtr[dataType]["hImpactParameterError"])->Fill(candidate.errorImpactParameter1(), pt); + std::get(histosPtr[dataType]["hImpactParameterError"])->Fill(candidate.errorImpactParameter2(), pt); + std::get(histosPtr[dataType]["hd0Prong0"])->Fill(candidate.impactParameter0(), pt); + std::get(histosPtr[dataType]["hd0Prong1"])->Fill(candidate.impactParameter1(), pt); + std::get(histosPtr[dataType]["hd0Prong2"])->Fill(candidate.impactParameter2(), pt); return; } @@ -382,7 +429,7 @@ struct HfTaskDs { /// \param candidate is candidate /// \param dataType is data class, as defined in DataType enum /// \param finalState is either KKPi or PiKK, as defined in FinalState enum - template + template void fillSparse(const Cand& candidate, DataType dataType, FinalState finalState) { auto mass = finalState == FinalState::KKPi ? hfHelper.invMassDsToKKPi(candidate) : hfHelper.invMassDsToPiKK(candidate); @@ -398,46 +445,88 @@ struct HfTaskDs { } if (dataType == DataType::Data) { // If data do not fill PV contributors in sparse - std::get(histosPtr[dataType]["hSparseMass"])->Fill(mass, pt, evaluateCentralityCand(candidate), outputMl[0], outputMl[1], outputMl[2]); - return; + if (storeOccupancy) { + std::get(histosPtr[dataType]["hSparseMass"])->Fill(mass, pt, evaluateCentralityCand(candidate), outputMl[0], outputMl[1], outputMl[2], o2::hf_occupancy::getOccupancyColl(candidate.template collision_as(), occEstimator)); + return; + } else { + std::get(histosPtr[dataType]["hSparseMass"])->Fill(mass, pt, evaluateCentralityCand(candidate), outputMl[0], outputMl[1], outputMl[2]); + return; + } + } + if constexpr (isMc) { + if (dataType == DataType::McDsNonPrompt) { // If data do not fill PV contributors in sparse + if (storeOccupancy) { + std::get(histosPtr[dataType]["hSparseMass"])->Fill(mass, pt, evaluateCentralityCand(candidate), outputMl[0], outputMl[1], outputMl[2], candidate.ptBhadMotherPart(), getBHadMotherFlag(candidate.pdgBhadMotherPart()), o2::hf_occupancy::getOccupancyColl(candidate.template collision_as(), occEstimator)); + return; + } else { + std::get(histosPtr[dataType]["hSparseMass"])->Fill(mass, pt, evaluateCentralityCand(candidate), outputMl[0], outputMl[1], outputMl[2], candidate.ptBhadMotherPart(), getBHadMotherFlag(candidate.pdgBhadMotherPart())); + return; + } + } else { + if (storeOccupancy) { + std::get(histosPtr[dataType]["hSparseMass"])->Fill(mass, pt, evaluateCentralityCand(candidate), candidate.template collision_as().numContrib(), outputMl[0], outputMl[1], outputMl[2], o2::hf_occupancy::getOccupancyColl(candidate.template collision_as(), occEstimator)); + return; + } else { + std::get(histosPtr[dataType]["hSparseMass"])->Fill(mass, pt, evaluateCentralityCand(candidate), candidate.template collision_as().numContrib(), outputMl[0], outputMl[1], outputMl[2]); + return; + } + } } - std::get(histosPtr[dataType]["hSparseMass"])->Fill(mass, pt, evaluateCentralityCand(candidate), outputMl[0], outputMl[1], outputMl[2], candidate.template collision_as().numContrib()); - - return; } /// Fill mass sparse if ML information is not present /// \param candidate is candidate /// \param dataType is data class, as defined in DataType enum /// \param finalState is either KKPi or PiKK, as defined in FinalState enum - template + template void fillSparse(const Cand& candidate, DataType dataType, FinalState finalState) { auto mass = finalState == FinalState::KKPi ? hfHelper.invMassDsToKKPi(candidate) : hfHelper.invMassDsToPiKK(candidate); auto pt = candidate.pt(); if (dataType == DataType::Data) { // If data do not fill PV contributors in sparse - std::get(histosPtr[dataType]["hSparseMass"])->Fill(mass, pt, evaluateCentralityCand(candidate)); - return; + if (storeOccupancy) { + std::get(histosPtr[dataType]["hSparseMass"])->Fill(mass, pt, evaluateCentralityCand(candidate), o2::hf_occupancy::getOccupancyColl(candidate.template collision_as(), occEstimator)); + return; + } else { + std::get(histosPtr[dataType]["hSparseMass"])->Fill(mass, pt, evaluateCentralityCand(candidate)); + return; + } + } + if constexpr (isMc) { + if (dataType == DataType::McDsNonPrompt) { // If data do not fill PV contributors in sparse + if (storeOccupancy) { + std::get(histosPtr[dataType]["hSparseMass"])->Fill(mass, pt, evaluateCentralityCand(candidate), candidate.ptBhadMotherPart(), getBHadMotherFlag(candidate.pdgBhadMotherPart()), o2::hf_occupancy::getOccupancyColl(candidate.template collision_as(), occEstimator)); + return; + } else { + std::get(histosPtr[dataType]["hSparseMass"])->Fill(mass, pt, evaluateCentralityCand(candidate), candidate.ptBhadMotherPart(), getBHadMotherFlag(candidate.pdgBhadMotherPart())); + return; + } + } else { + if (storeOccupancy) { + std::get(histosPtr[dataType]["hSparseMass"])->Fill(mass, pt, evaluateCentralityCand(candidate), candidate.template collision_as().numContrib(), o2::hf_occupancy::getOccupancyColl(candidate.template collision_as(), occEstimator)); + return; + } else { + std::get(histosPtr[dataType]["hSparseMass"])->Fill(mass, pt, evaluateCentralityCand(candidate), candidate.template collision_as().numContrib()); + return; + } + } } - std::get(histosPtr[dataType]["hSparseMass"])->Fill(mass, pt, evaluateCentralityCand(candidate), candidate.template collision_as().numContrib()); - - return; } /// Fill histograms of quantities for the KKPi daugther-mass hypothesis /// \param candidate is candidate /// \param dataType is data class, as defined in DataType enum - template + template void fillHistoKKPi(const T1& candidate, DataType dataType) { auto pt = candidate.pt(); - fillSparse(candidate, dataType, FinalState::KKPi); + fillSparse(candidate, dataType, FinalState::KKPi); - std::get(histosPtr[dataType]["hCos3PiK"])->Fill(hfHelper.cos3PiKDsToKKPi(candidate), pt); - std::get(histosPtr[dataType]["hAbsCos3PiK"])->Fill(hfHelper.absCos3PiKDsToKKPi(candidate), pt); - std::get(histosPtr[dataType]["hDeltaMassPhi"])->Fill(hfHelper.deltaMassPhiDsToKKPi(candidate), pt); - std::get(histosPtr[dataType]["hMassKK"])->Fill(hfHelper.massKKPairDsToKKPi(candidate), pt); + std::get(histosPtr[dataType]["hCos3PiK"])->Fill(hfHelper.cos3PiKDsToKKPi(candidate), pt); + std::get(histosPtr[dataType]["hAbsCos3PiK"])->Fill(hfHelper.absCos3PiKDsToKKPi(candidate), pt); + std::get(histosPtr[dataType]["hDeltaMassPhi"])->Fill(hfHelper.deltaMassPhiDsToKKPi(candidate), pt); + std::get(histosPtr[dataType]["hMassKK"])->Fill(hfHelper.massKKPairDsToKKPi(candidate), pt); return; } @@ -445,16 +534,16 @@ struct HfTaskDs { /// Fill histograms of quantities for the PiKK daugther-mass hypothesis /// \param candidate is candidate /// \param dataType is data class, as defined in DataType enum - template + template void fillHistoPiKK(const T1& candidate, DataType dataType) { auto pt = candidate.pt(); - fillSparse(candidate, dataType, FinalState::PiKK); + fillSparse(candidate, dataType, FinalState::PiKK); - std::get(histosPtr[dataType]["hCos3PiK"])->Fill(hfHelper.cos3PiKDsToPiKK(candidate), pt); - std::get(histosPtr[dataType]["hAbsCos3PiK"])->Fill(hfHelper.absCos3PiKDsToPiKK(candidate), pt); - std::get(histosPtr[dataType]["hDeltaMassPhi"])->Fill(hfHelper.deltaMassPhiDsToPiKK(candidate), pt); - std::get(histosPtr[dataType]["hMassKK"])->Fill(hfHelper.massKKPairDsToPiKK(candidate), pt); + std::get(histosPtr[dataType]["hCos3PiK"])->Fill(hfHelper.cos3PiKDsToPiKK(candidate), pt); + std::get(histosPtr[dataType]["hAbsCos3PiK"])->Fill(hfHelper.absCos3PiKDsToPiKK(candidate), pt); + std::get(histosPtr[dataType]["hDeltaMassPhi"])->Fill(hfHelper.deltaMassPhiDsToPiKK(candidate), pt); + std::get(histosPtr[dataType]["hMassKK"])->Fill(hfHelper.massKKPairDsToPiKK(candidate), pt); return; } @@ -490,30 +579,30 @@ struct HfTaskDs { if (candidate.isSelDsToKKPi() >= selectionFlagDs) { // KKPi fillHisto(candidate, dataType); - fillHistoKKPi(candidate, dataType); + fillHistoKKPi(candidate, dataType); if (TESTBIT(candidate.isSelDsToKKPi(), aod::SelectionStep::RecoSkims)) { - std::get(histosPtr[dataType]["hPtVsYRecoSkim"])->Fill(pt, yCand); + std::get(histosPtr[dataType]["hPtVsYRecoSkim"])->Fill(pt, yCand); } if (TESTBIT(candidate.isSelDsToKKPi(), aod::SelectionStep::RecoTopol)) { - std::get(histosPtr[dataType]["hPtVsYRecoTopol"])->Fill(pt, yCand); + std::get(histosPtr[dataType]["hPtVsYRecoTopol"])->Fill(pt, yCand); } if (TESTBIT(candidate.isSelDsToKKPi(), aod::SelectionStep::RecoPID)) { - std::get(histosPtr[dataType]["hPtVsYRecoPID"])->Fill(pt, yCand); + std::get(histosPtr[dataType]["hPtVsYRecoPID"])->Fill(pt, yCand); } } if (candidate.isSelDsToPiKK() >= selectionFlagDs) { // PiKK fillHisto(candidate, dataType); - fillHistoPiKK(candidate, dataType); + fillHistoPiKK(candidate, dataType); if (TESTBIT(candidate.isSelDsToPiKK(), aod::SelectionStep::RecoSkims)) { - std::get(histosPtr[dataType]["hPtVsYRecoSkim"])->Fill(pt, yCand); + std::get(histosPtr[dataType]["hPtVsYRecoSkim"])->Fill(pt, yCand); } if (TESTBIT(candidate.isSelDsToPiKK(), aod::SelectionStep::RecoTopol)) { - std::get(histosPtr[dataType]["hPtVsYRecoTopol"])->Fill(pt, yCand); + std::get(histosPtr[dataType]["hPtVsYRecoTopol"])->Fill(pt, yCand); } if (TESTBIT(candidate.isSelDsToPiKK(), aod::SelectionStep::RecoPID)) { - std::get(histosPtr[dataType]["hPtVsYRecoPID"])->Fill(pt, yCand); + std::get(histosPtr[dataType]["hPtVsYRecoPID"])->Fill(pt, yCand); } } } @@ -529,11 +618,11 @@ struct HfTaskDs { if (candidate.isSelDsToKKPi() >= selectionFlagDs) { // KKPi fillHisto(candidate, DataType::Data); - fillHistoKKPi(candidate, DataType::Data); + fillHistoKKPi(candidate, DataType::Data); } if (candidate.isSelDsToPiKK() >= selectionFlagDs) { // PiKK fillHisto(candidate, DataType::Data); - fillHistoPiKK(candidate, DataType::Data); + fillHistoPiKK(candidate, DataType::Data); } } @@ -566,11 +655,11 @@ struct HfTaskDs { if (candidate.isSelDsToKKPi() >= selectionFlagDs || candidate.isSelDsToPiKK() >= selectionFlagDs) { if (candidate.isSelDsToKKPi() >= selectionFlagDs) { // KKPi fillHisto(candidate, DataType::McBkg); - fillHistoKKPi(candidate, DataType::McBkg); + fillHistoKKPi(candidate, DataType::McBkg); } if (candidate.isSelDsToPiKK() >= selectionFlagDs) { // PiKK fillHisto(candidate, DataType::McBkg); - fillHistoPiKK(candidate, DataType::McBkg); + fillHistoPiKK(candidate, DataType::McBkg); } } } @@ -579,21 +668,28 @@ struct HfTaskDs { } template - void fillMcGenHistos(CandDsMcGen const& mcParticles, - Coll const& recoCollisions) + void fillMcGenHistosSparse(CandDsMcGen const& mcParticles, + Coll const& recoCollisions) { + // MC gen. for (const auto& particle : mcParticles) { + const auto& recoCollsPerMcColl = recoCollisions.sliceBy(colPerMcCollision, particle.mcCollision().globalIndex()); + if (std::abs(particle.flagMcMatchGen()) == 1 << aod::hf_cand_3prong::DecayType::DsToKKPi) { if (particle.flagMcDecayChanGen() == decayChannel || (fillDplusMc && particle.flagMcDecayChanGen() == (decayChannel + offsetDplusDecayChannel))) { auto pt = particle.pt(); double y{0.f}; unsigned maxNumContrib = 0; - const auto& recoCollsPerMcColl = recoCollisions.sliceBy(colPerMcCollision, particle.mcCollision().globalIndex()); for (const auto& recCol : recoCollsPerMcColl) { maxNumContrib = recCol.numContrib() > maxNumContrib ? recCol.numContrib() : maxNumContrib; } + float cent = o2::hf_centrality::getCentralityGenColl(recoCollsPerMcColl); + float occ{-1.}; + if (storeOccupancy && occEstimator != o2::hf_occupancy::OccupancyEstimator::None) { + occ = o2::hf_occupancy::getOccupancyGenColl(recoCollsPerMcColl, occEstimator); + } if (particle.flagMcDecayChanGen() == decayChannel) { y = RecoDecay::y(particle.pVector(), o2::constants::physics::MassDS); @@ -602,14 +698,25 @@ struct HfTaskDs { } if (particle.originMcGen() == RecoDecay::OriginType::Prompt) { - std::get(histosPtr[DataType::McDsPrompt]["hPtGen"])->Fill(pt); // gen. level pT - std::get(histosPtr[DataType::McDsPrompt]["hEtaGen"])->Fill(particle.eta()); - std::get(histosPtr[DataType::McDsPrompt]["hPtYNPvContribGen"])->Fill(pt, y, maxNumContrib); // gen. level pT + std::get(histosPtr[DataType::McDsPrompt]["hPtGen"])->Fill(pt); // gen. level pT + std::get(histosPtr[DataType::McDsPrompt]["hEtaGen"])->Fill(particle.eta()); + if (storeOccupancy && occEstimator != o2::hf_occupancy::OccupancyEstimator::None) { + std::get(histosPtr[DataType::McDsPrompt]["hSparseGen"])->Fill(pt, y, maxNumContrib, cent, occ); + } else { + std::get(histosPtr[DataType::McDsPrompt]["hSparseGen"])->Fill(pt, y, maxNumContrib, cent); + } } if (particle.originMcGen() == RecoDecay::OriginType::NonPrompt) { - std::get(histosPtr[DataType::McDsNonPrompt]["hPtGen"])->Fill(pt); // gen. level pT - std::get(histosPtr[DataType::McDsNonPrompt]["hEtaGen"])->Fill(particle.eta()); // gen. level pT - std::get(histosPtr[DataType::McDsNonPrompt]["hPtYNPvContribGen"])->Fill(pt, y, maxNumContrib); // gen. level pT + std::get(histosPtr[DataType::McDsNonPrompt]["hPtGen"])->Fill(pt); // gen. level pT + std::get(histosPtr[DataType::McDsNonPrompt]["hEtaGen"])->Fill(particle.eta()); + auto bHadMother = mcParticles.rawIteratorAt(particle.idxBhadMotherPart() - mcParticles.offset()); + int flagGenB = getBHadMotherFlag(bHadMother.pdgCode()); + float ptGenB = bHadMother.pt(); + if (storeOccupancy && occEstimator != o2::hf_occupancy::OccupancyEstimator::None) { + std::get(histosPtr[DataType::McDsNonPrompt]["hSparseGen"])->Fill(pt, y, maxNumContrib, ptGenB, flagGenB, cent, occ); + } else { + std::get(histosPtr[DataType::McDsNonPrompt]["hSparseGen"])->Fill(pt, y, maxNumContrib, ptGenB, flagGenB, cent); + } } } else if (fillDplusMc) { y = RecoDecay::y(particle.pVector(), o2::constants::physics::MassDPlus); @@ -617,17 +724,45 @@ struct HfTaskDs { continue; } if (particle.originMcGen() == RecoDecay::OriginType::Prompt) { - std::get(histosPtr[DataType::McDplusPrompt]["hPtGen"])->Fill(pt); // gen. level pT - std::get(histosPtr[DataType::McDplusPrompt]["hEtaGen"])->Fill(particle.eta()); - std::get(histosPtr[DataType::McDplusPrompt]["hPtYNPvContribGen"])->Fill(pt, y, maxNumContrib); // gen. level pT + std::get(histosPtr[DataType::McDplusPrompt]["hPtGen"])->Fill(pt); // gen. level pT + std::get(histosPtr[DataType::McDplusPrompt]["hEtaGen"])->Fill(particle.eta()); + if (storeOccupancy && occEstimator != o2::hf_occupancy::OccupancyEstimator::None) { + std::get(histosPtr[DataType::McDplusPrompt]["hSparseGen"])->Fill(pt, y, maxNumContrib, cent, occ); + } else { + std::get(histosPtr[DataType::McDplusPrompt]["hSparseGen"])->Fill(pt, y, maxNumContrib, cent); + } } if (particle.originMcGen() == RecoDecay::OriginType::NonPrompt) { - std::get(histosPtr[DataType::McDplusNonPrompt]["hPtGen"])->Fill(pt); // gen. level pT - std::get(histosPtr[DataType::McDplusNonPrompt]["hEtaGen"])->Fill(particle.eta()); - std::get(histosPtr[DataType::McDplusNonPrompt]["hPtYNPvContribGen"])->Fill(pt, y, maxNumContrib); // gen. level pT + std::get(histosPtr[DataType::McDplusNonPrompt]["hPtGen"])->Fill(pt); // gen. level pT + std::get(histosPtr[DataType::McDplusNonPrompt]["hEtaGen"])->Fill(particle.eta()); + auto bHadMother = mcParticles.rawIteratorAt(particle.idxBhadMotherPart() - mcParticles.offset()); + int flagGenB = getBHadMotherFlag(bHadMother.pdgCode()); + float ptGenB = bHadMother.pt(); + if (storeOccupancy && occEstimator != o2::hf_occupancy::OccupancyEstimator::None) { + std::get(histosPtr[DataType::McDplusNonPrompt]["hSparseGen"])->Fill(pt, y, maxNumContrib, ptGenB, flagGenB, cent, occ); + } else { + std::get(histosPtr[DataType::McDplusNonPrompt]["hSparseGen"])->Fill(pt, y, maxNumContrib, ptGenB, flagGenB, cent); + } } } } + } else { // not matched candidates + auto pt = particle.pt(); + double y = RecoDecay::y(particle.pVector(), o2::constants::physics::MassDS); + unsigned maxNumContrib = 0; + for (const auto& recCol : recoCollsPerMcColl) { + maxNumContrib = recCol.numContrib() > maxNumContrib ? recCol.numContrib() : maxNumContrib; + } + float cent = o2::hf_centrality::getCentralityGenColl(recoCollsPerMcColl); + float occ{-1.}; + if (storeOccupancy && occEstimator != o2::hf_occupancy::OccupancyEstimator::None) { + occ = o2::hf_occupancy::getOccupancyGenColl(recoCollsPerMcColl, occEstimator); + } + if (storeOccupancy && occEstimator != o2::hf_occupancy::OccupancyEstimator::None) { + std::get(histosPtr[DataType::McBkg]["hSparseGen"])->Fill(pt, y, maxNumContrib, cent, occ); + } else { + std::get(histosPtr[DataType::McBkg]["hSparseGen"])->Fill(pt, y, maxNumContrib, cent); + } } } } @@ -640,16 +775,16 @@ struct HfTaskDs { { int numPvContributors = collision.numContrib(); float centrality = evaluateCentralityColl(collision); - std::get(histosPtr[DataType::Data]["hNPvContribAll"])->Fill(numPvContributors, centrality); + std::get(histosPtr[DataType::Data]["hNPvContribAll"])->Fill(numPvContributors, centrality); for (int i = 0; i < DataType::kDataTypes; i++) { if (nCandsPerType[i]) { - std::get(histosPtr[i]["hNPvContribCands"])->Fill(numPvContributors, centrality); + std::get(histosPtr[i]["hNPvContribCands"])->Fill(numPvContributors, centrality); } if (nCandsInSignalRegionDsPerType[i]) { - std::get(histosPtr[i]["hNPvContribCandsInSignalRegionDs"])->Fill(numPvContributors, centrality); + std::get(histosPtr[i]["hNPvContribCandsInSignalRegionDs"])->Fill(numPvContributors, centrality); } if (nCandsInSignalRegionDplusPerType[i]) { - std::get(histosPtr[i]["hNPvContribCandsInSignalRegionDplus"])->Fill(numPvContributors, centrality); + std::get(histosPtr[i]["hNPvContribCandsInSignalRegionDplus"])->Fill(numPvContributors, centrality); } } } @@ -809,7 +944,7 @@ struct HfTaskDs { } fillNPvContribHisto(collision, nCandsPerType, nCandsInSignalRegionDsPerType, nCandsInSignalRegionDplusPerType); } - fillMcGenHistos(mcParticles, collisions); + fillMcGenHistosSparse(mcParticles, collisions); } void processDataWithCentFT0C(CollisionsWithFT0C const& collisions, diff --git a/PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx b/PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx index 2e6544d684a..a154b416b2b 100644 --- a/PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx +++ b/PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx @@ -37,6 +37,7 @@ using namespace o2::aod; using namespace o2::framework; using namespace o2::framework::expressions; using namespace o2::hf_centrality; +using namespace o2::hf_occupancy; using namespace o2::hf_evsel; enum DecayChannel { DplusToPiKPi = 0, @@ -330,32 +331,6 @@ struct HfTaskFlowCharmHadrons { } } - /// Get the centrality - /// \param collision is the collision with the centrality information - float getCentrality(CollsWithQvecs::iterator const& collision) - { - float 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; - } - /// Check if the collision is selected /// \param collision is the collision with the Q vector information /// \param bc is the bunch crossing with timestamp information @@ -366,9 +341,9 @@ struct HfTaskFlowCharmHadrons { aod::BCsWithTimestamps const&, float& centrality) { - float occupancy = hfEvSel.getOccupancy(collision, occEstimator); + float occupancy = getOccupancyColl(collision, occEstimator); const auto rejectionMask = hfEvSel.getHfCollisionRejectionMask(collision, centrality, ccdb, registry); - centrality = getCentrality(collision); + centrality = o2::hf_centrality::getCentralityColl(collision, centEstimator); /// monitor the satisfied event selections hfEvSel.fillHistograms(collision, rejectionMask, centrality, occupancy); @@ -431,14 +406,14 @@ struct HfTaskFlowCharmHadrons { void runFlowAnalysis(CollsWithQvecs::iterator const& collision, T1 const& candidates) { - float cent = getCentrality(collision); + float cent = o2::hf_centrality::getCentralityColl(collision, centEstimator); if (cent < centralityMin || cent > centralityMax) { return; } float occupancy = 0.; uint16_t hfevflag{}; if (occEstimator != 0) { - occupancy = hfEvSel.getOccupancy(collision, occEstimator); + occupancy = getOccupancyColl(collision, occEstimator); registry.fill(HIST("trackOccVsFT0COcc"), collision.trackOccupancyInTimeRange(), collision.ft0cOccupancyInTimeRange()); hfevflag = hfEvSel.getHfCollisionRejectionMask(collision, cent, ccdb, registry); } @@ -632,7 +607,7 @@ struct HfTaskFlowCharmHadrons { aod::BCsWithTimestamps const& bcs) { float centrality{-1.f}; - if (!isCollSelected(collision, bcs, centrality)) { + if (!isCollSelected(collision, bcs, centrality)) { // no selection on the centrality is applied on purpose to allow for the resolution study in post-processing return; } diff --git a/PWGHF/TableProducer/candidateCreator2Prong.cxx b/PWGHF/TableProducer/candidateCreator2Prong.cxx index a151ac4e975..d12acd58271 100644 --- a/PWGHF/TableProducer/candidateCreator2Prong.cxx +++ b/PWGHF/TableProducer/candidateCreator2Prong.cxx @@ -56,6 +56,7 @@ using namespace o2::hf_evsel; using namespace o2::hf_trkcandsel; using namespace o2::aod::hf_cand_2prong; using namespace o2::hf_centrality; +using namespace o2::hf_occupancy; using namespace o2::constants::physics; using namespace o2::framework; using namespace o2::aod::pid_tpc_tof_utils; @@ -627,7 +628,7 @@ struct HfCandidateCreator2Prong { /// bitmask with event. selection info float centrality{-1.f}; - float occupancy = hfEvSel.getOccupancy(collision); + float occupancy = getOccupancyColl(collision, OccupancyEstimator::Its); const auto rejectionMask = hfEvSel.getHfCollisionRejectionMask(collision, centrality, ccdb, registry); /// monitor the satisfied event selections @@ -645,7 +646,7 @@ struct HfCandidateCreator2Prong { /// bitmask with event. selection info float centrality{-1.f}; - float occupancy = hfEvSel.getOccupancy(collision); + float occupancy = getOccupancyColl(collision, OccupancyEstimator::Its); const auto rejectionMask = hfEvSel.getHfCollisionRejectionMask(collision, centrality, ccdb, registry); /// monitor the satisfied event selections @@ -663,7 +664,7 @@ struct HfCandidateCreator2Prong { /// bitmask with event. selection info float centrality{-1.f}; - float occupancy = hfEvSel.getOccupancy(collision); + float occupancy = getOccupancyColl(collision, OccupancyEstimator::Its); const auto rejectionMask = hfEvSel.getHfCollisionRejectionMask(collision, centrality, ccdb, registry); /// monitor the satisfied event selections diff --git a/PWGHF/TableProducer/candidateCreator3Prong.cxx b/PWGHF/TableProducer/candidateCreator3Prong.cxx index dd1b7b3e07b..d59cf2989b2 100644 --- a/PWGHF/TableProducer/candidateCreator3Prong.cxx +++ b/PWGHF/TableProducer/candidateCreator3Prong.cxx @@ -44,6 +44,7 @@ using namespace o2::hf_evsel; using namespace o2::hf_trkcandsel; using namespace o2::aod::hf_cand_3prong; using namespace o2::hf_centrality; +using namespace o2::hf_occupancy; using namespace o2::constants::physics; using namespace o2::framework; using namespace o2::framework::expressions; @@ -405,7 +406,7 @@ struct HfCandidateCreator3Prong { /// bitmask with event. selection info float centrality{-1.f}; - float occupancy = hfEvSel.getOccupancy(collision); + float occupancy = getOccupancyColl(collision, OccupancyEstimator::Its); const auto rejectionMask = hfEvSel.getHfCollisionRejectionMask(collision, centrality, ccdb, registry); /// monitor the satisfied event selections @@ -423,7 +424,7 @@ struct HfCandidateCreator3Prong { /// bitmask with event. selection info float centrality{-1.f}; - float occupancy = hfEvSel.getOccupancy(collision); + float occupancy = getOccupancyColl(collision, OccupancyEstimator::Its); const auto rejectionMask = hfEvSel.getHfCollisionRejectionMask(collision, centrality, ccdb, registry); /// monitor the satisfied event selections @@ -441,7 +442,7 @@ struct HfCandidateCreator3Prong { /// bitmask with event. selection info float centrality{-1.f}; - float occupancy = hfEvSel.getOccupancy(collision); + float occupancy = getOccupancyColl(collision, OccupancyEstimator::Its); const auto rejectionMask = hfEvSel.getHfCollisionRejectionMask(collision, centrality, ccdb, registry); /// monitor the satisfied event selections diff --git a/PWGHF/Tasks/taskMcValidation.cxx b/PWGHF/Tasks/taskMcValidation.cxx index 558fe20f90f..38666bcba6a 100644 --- a/PWGHF/Tasks/taskMcValidation.cxx +++ b/PWGHF/Tasks/taskMcValidation.cxx @@ -44,6 +44,7 @@ using namespace o2::framework; using namespace o2::framework::expressions; using namespace o2::hf_evsel; using namespace o2::hf_centrality; +using namespace o2::hf_occupancy; namespace { @@ -242,27 +243,6 @@ struct HfTaskMcValidationGen { hfEvSelMc.addHistograms(registry); // particles monitoring } - /// \brief Function to get MC collision occupancy - /// \param collSlice collection of reconstructed collisions - /// \return collision occupancy - template - int getOccupancyColl(CCs const& collSlice) - { - float multiplicity{0.f}; - int occupancy = 0; - for (const auto& collision : collSlice) { - float collMult{0.f}; - collMult = collision.numContrib(); - - if (collMult > multiplicity) { - occupancy = collision.trackOccupancyInTimeRange(); - multiplicity = collMult; - } - } // end loop over collisions - - return occupancy; - } - template void runCheckGenParticles(GenColl const& mcCollision, Particles const& mcParticles, RecoColls const& recoCollisions, BCsInfo const&, std::array& counterPrompt, std::array& counterNonPrompt) { @@ -275,7 +255,7 @@ struct HfTaskMcValidationGen { float centrality{105.f}; int occupancy = 0; if (storeOccupancy) { - occupancy = getOccupancyColl(recoCollisions); + occupancy = getOccupancyGenColl(recoCollisions, OccupancyEstimator::Its); } uint16_t rejectionMask{0}; if constexpr (centEstimator == CentralityEstimator::FT0C) { diff --git a/PWGHF/Utils/utilsAnalysis.h b/PWGHF/Utils/utilsAnalysis.h index b10a851072d..9148076e0ec 100644 --- a/PWGHF/Utils/utilsAnalysis.h +++ b/PWGHF/Utils/utilsAnalysis.h @@ -21,6 +21,33 @@ namespace o2::analysis { + +enum BHadMothers { NotMatched = 0, + BPlus, + BZero, + Bs, + LambdaBZero }; + +/// Convert the B hadron mother PDG for non prompt candidates to a flag +/// \param pdg of the b hadron mother +/// \return integer map to specific mothers' PDG codes +BHadMothers getBHadMotherFlag(const int flagBHad) +{ + if (std::abs(flagBHad) == o2::constants::physics::kBPlus) { + return BHadMothers::BPlus; + } + if (std::abs(flagBHad) == o2::constants::physics::kB0) { + return BHadMothers::BZero; + } + if (std::abs(flagBHad) == o2::constants::physics::kBS) { + return BHadMothers::Bs; + } + if (std::abs(flagBHad) == o2::constants::physics::kLambdaB0) { + return BHadMothers::LambdaBZero; + } + return BHadMothers::NotMatched; +} + /// Finds pT bin in an array. /// \param bins array of pT bins /// \param value pT diff --git a/PWGHF/Utils/utilsEvSelHf.h b/PWGHF/Utils/utilsEvSelHf.h index e1c7ced8fa8..cba0f802c3d 100644 --- a/PWGHF/Utils/utilsEvSelHf.h +++ b/PWGHF/Utils/utilsEvSelHf.h @@ -32,6 +32,51 @@ #include "PWGLF/DataModel/mcCentrality.h" #include "PWGHF/Core/CentralityEstimation.h" +namespace o2::hf_occupancy +{ +// centrality selection estimators +enum OccupancyEstimator { None = 0, + Its, + Ft0c }; + +/// Get the occupancy +/// \param collision is the collision with the occupancy information +/// \return collision occupancy +template +float getOccupancyColl(Coll const& collision, int occEstimator) +{ + switch (occEstimator) { + case OccupancyEstimator::Its: + return collision.trackOccupancyInTimeRange(); + case OccupancyEstimator::Ft0c: + return collision.ft0cOccupancyInTimeRange(); + default: + LOG(fatal) << "Occupancy estimator not valid. Possible values are ITS or FT0C."; + break; + } + return -999.f; +}; + +/// \brief Function to get MC collision occupancy +/// \param collSlice collection of reconstructed collisions associated to a generated one +/// \return generated MC collision occupancy +template +int getOccupancyGenColl(CCs const& collSlice, int occEstimator) +{ + float multiplicity{0.f}; + int occupancy = 0; + for (const auto& collision : collSlice) { + float collMult{0.f}; + collMult = collision.numContrib(); + if (collMult > multiplicity) { + occupancy = getOccupancyColl(collision, occEstimator); + multiplicity = collMult; + } + } // end loop over collisions + return occupancy; +}; +} // namespace o2::hf_occupancy + namespace o2::hf_evsel { // event rejection types @@ -166,7 +211,7 @@ struct HfEventSelection : o2::framework::ConfigurableGroup { uint16_t rejectionMask{0}; // 16 bits, in case new ev. selections will be added if constexpr (centEstimator != o2::hf_centrality::CentralityEstimator::None) { - centrality = getCentrality(collision); + centrality = o2::hf_centrality::getCentralityColl(collision, centEstimator); if (centrality < centralityMin || centrality > centralityMax) { SETBIT(rejectionMask, EventRejection::Centrality); } @@ -211,7 +256,7 @@ struct HfEventSelection : o2::framework::ConfigurableGroup { SETBIT(rejectionMask, EventRejection::NoCollInRofStandard); } if (useOccupancyCut) { - float occupancy = getOccupancy(collision, occEstimator); + float occupancy = o2::hf_occupancy::getOccupancyColl(collision, occEstimator); if (occupancy < occupancyMin || occupancy > occupancyMax) { SETBIT(rejectionMask, EventRejection::Occupancy); } @@ -260,46 +305,6 @@ struct HfEventSelection : o2::framework::ConfigurableGroup { return rejectionMask; } - /// Get the occupancy - /// \param collision is the collision with the occupancy information - /// \param occEstimator is the occupancy estimator (1: ITS, 2: FT0C) - template - float getOccupancy(Coll const& collision, int occEstimator = 1) - { - switch (occEstimator) { - case 1: // ITS - return collision.trackOccupancyInTimeRange(); - break; - case 2: // FT0c - return collision.ft0cOccupancyInTimeRange(); - break; - default: - LOG(warning) << "Occupancy estimator not valid. Possible values are ITS or FT0C. Fallback to ITS"; - return collision.trackOccupancyInTimeRange(); - break; - } - } - - /// Get the centrality - /// \param collision is the collision with the centrality information - /// \param centEstimator is the centrality estimator from hf_centrality::CentralityEstimator - template - float getCentrality(Coll const& collision) - { - if constexpr (centEstimator == o2::hf_centrality::CentralityEstimator::FT0A) { - return collision.centFT0A(); - } else if constexpr (centEstimator == o2::hf_centrality::CentralityEstimator::FT0C) { - return collision.centFT0C(); - } else if constexpr (centEstimator == o2::hf_centrality::CentralityEstimator::FT0M) { - return collision.centFT0M(); - } else if constexpr (centEstimator == o2::hf_centrality::CentralityEstimator::FV0A) { - return collision.centFV0A(); - } else { - LOG(warning) << "Centrality estimator not valid. Possible values are V0A, T0M, T0A, T0C. Fallback to FT0c"; - return collision.centFT0C(); - } - } - /// \brief Fills histograms for monitoring event selections satisfied by the collision. /// \param collision analysed collision /// \param rejectionMask bitmask storing the info about which ev. selections are not satisfied by the collision @@ -389,27 +394,7 @@ struct HfEventSelectionMc { auto bc = mcCollision.template bc_as(); if constexpr (centEstimator != o2::hf_centrality::CentralityEstimator::None) { - float multiplicity{0.f}; - for (const auto& collision : collSlice) { - float collCent{0.f}; - float collMult{0.f}; - if constexpr (centEstimator == o2::hf_centrality::CentralityEstimator::FT0A) { - collCent = collision.centFT0A(); - } else if constexpr (centEstimator == o2::hf_centrality::CentralityEstimator::FT0C) { - collCent = collision.centFT0C(); - } else if constexpr (centEstimator == o2::hf_centrality::CentralityEstimator::FT0M) { - collCent = collision.centFT0M(); - } else if constexpr (centEstimator == o2::hf_centrality::CentralityEstimator::FV0A) { - collCent = collision.centFV0A(); - } else { - LOGP(fatal, "Unsupported centrality estimator!"); - } - collMult = collision.numContrib(); - if (collMult > multiplicity) { - centrality = collCent; - multiplicity = collMult; - } - } + centrality = o2::hf_centrality::getCentralityGenColl(collSlice, centEstimator); /// centrality selection if (centrality < centralityMin || centrality > centralityMax) { SETBIT(rejectionMask, EventRejection::Centrality);