diff --git a/PWGHF/D2H/Tasks/taskCharmPolarisation.cxx b/PWGHF/D2H/Tasks/taskCharmPolarisation.cxx index 86cdbdef317..827e6396462 100644 --- a/PWGHF/D2H/Tasks/taskCharmPolarisation.cxx +++ b/PWGHF/D2H/Tasks/taskCharmPolarisation.cxx @@ -215,9 +215,8 @@ struct HfTaskCharmPolarisation { EventPlaneHelper epHelper; using CollisionsWithMcLabels = soa::SmallGroups>; + using CollisionsWithMcLabelsAndCent = soa::SmallGroups>; using CollsWithQVecs = soa::Join; - using CollsWithQVecsWithMcLabels = soa::SmallGroups>; - using GenCollisWithQvecs = soa::Join; using TracksWithMcLabels = soa::Join; using TracksWithExtra = soa::Join; @@ -250,6 +249,8 @@ struct HfTaskCharmPolarisation { Preslice lcToPKPiWithMcPerCollision = aod::hf_cand::collisionId; Preslice lcToPKPiWithMcAndMlPerCollision = aod::hf_cand::collisionId; + PresliceUnsorted colPerMcCollision = aod::mcparticle::mcCollisionId; + ConfigurableAxis configTHnAxisEulerPhi{"configTHnAxisEulerPhi", {24, -o2::constants::math::PI, o2::constants::math::PI}, "Euler polar angle #phi"}; ConfigurableAxis configThnAxisInvMass{"configThnAxisInvMass", {200, 0.139f, 0.179f}, "#it{M} (GeV/#it{c}^{2})"}; // o2-linter: disable=pdg/explicit-mass (false positive) ConfigurableAxis configThnAxisPt{"configThnAxisPt", {100, 0.f, 100.f}, "#it{p}_{T} (GeV/#it{c})"}; @@ -271,6 +272,7 @@ struct HfTaskCharmPolarisation { ConfigurableAxis configThnAxisNumItsClsMin{"configThnAxisNumItsClsMin", {4, 3.5f, 7.5f}, "min #it{N}_{cls ITS}"}; ConfigurableAxis configThnAxisNumTpcClsMin{"configThnAxisNumTpcClsMin", {3, 79.5f, 140.5f}, "min #it{N}_{cls TPC}"}; ConfigurableAxis configThnAxisCharge{"configThnAxisCharge", {2, -2.f, 2.f}, "electric charge"}; + ConfigurableAxis configThnAxisCentrality{"configThnAxisCentrality", {100, 0.f, 100.f}, "centrality (%)"}; HistogramRegistry registry{"registry", {}}; @@ -307,12 +309,12 @@ struct HfTaskCharmPolarisation { } } - if (activatePartRecoDstar && !(doprocessDstarMc || doprocessDstarMcWithMl)) { + if (activatePartRecoDstar && !(doprocessDstarMc || doprocessDstarMcWithMl || doprocessDstarMcInPbPb || doprocessDstarMcWithMlInPbPb)) { LOGP(fatal, "Check on partly reconstructed D* mesons only possible for processDstarMc and processDstarMcWithMl"); } // check bkg rotation for MC (not supported currently) - if (nBkgRotations > 0 && (doprocessDstarMc || doprocessDstarMcWithMl || doprocessLcToPKPiMc || doprocessLcToPKPiMcWithMl || doprocessLcToPKPiBackgroundMcWithMl)) { + if (nBkgRotations > 0 && (doprocessDstarMc || doprocessDstarMcWithMl || doprocessLcToPKPiMc || doprocessLcToPKPiMcWithMl || doprocessLcToPKPiBackgroundMcWithMl || doprocessDstarMcInPbPb || doprocessDstarMcWithMlInPbPb)) { LOGP(fatal, "No background rotation supported for MC."); } @@ -349,6 +351,7 @@ struct HfTaskCharmPolarisation { const AxisSpec thnAxisInvMass2PKLcMonitoring{lcPKPiChannels.configThnAxisInvMass2PKLcMonitoring, "#it{M}^{2}(pK) from #Lambda_{c}^{+} (GeV/#it{c}^{2})"}; const AxisSpec thnAxisInvMass2PPiLcMonitoring{lcPKPiChannels.configThnAxisInvMass2PPiLcMonitoring, "#it{M}^{2}(p#pi) from #Lambda_{c}^{+} (GeV/#it{c}^{2})"}; const AxisSpec thnAxisTHnAxisEulerPhi{configTHnAxisEulerPhi, "Euler polar angle #phi"}; + const AxisSpec thnAxisCentrality{configThnAxisCentrality, "centrality (%)"}; auto invMassBins = thnAxisInvMass.binEdges; minInvMass = invMassBins.front(); @@ -357,6 +360,9 @@ struct HfTaskCharmPolarisation { registry.add("hNumPvContributorsAll", "Number of PV contributors for all events ;num. PV contributors; counts", HistType::kTH1F, {thnAxisNumPvContributors}); registry.add("hNumPvContributorsCand", "Number of PV contributors for events with candidates;num. PV contributors; counts", HistType::kTH1F, {thnAxisNumPvContributors}); registry.add("hNumPvContributorsCandInMass", "Number of PV contributors for events with candidates in the signal region;num. PV contributors; counts", HistType::kTH1F, {thnAxisNumPvContributors}); + if (doprocessDstarInPbPb || doprocessDstarMcInPbPb || doprocessDstarWithMlInPbPb || doprocessDstarMcWithMlInPbPb) { + registry.add("hCentrality", "Centrality distribution for D*+ candidates;centrality (%); counts", HistType::kTH1D, {thnAxisCentrality}); + } if (activateTHnSparseCosThStarHelicity) { std::vector hHelicityaxes = {thnAxisInvMass, thnAxisPt, thnAxisNumPvContributors, thnAxisY}; @@ -426,8 +432,8 @@ struct HfTaskCharmPolarisation { registry.add("hGenPromptHelicity", "THn for polarisation studies with cosThStar w.r.t. helicity axis and BDT scores for generated prompt D*+ candidates", HistType::kTHnSparseF, hgenPromptAxes); registry.add("hGenNonPromptHelicity", "THn for polarisation studies with cosThStar w.r.t. helicity axis and BDT scores for generated non-prompt D*+ candidates", HistType::kTHnSparseF, hgenNonPromptAxes); if (activatePartRecoDstar) { - registry.add("hPartRecoGenPromptHelicity", "THn for polarisation studies with cosThStar w.r.t. helicity axis and BDT scores for partially reconstructed generated prompt D*+ candidates", HistType::kTHnSparseF, hgenPromptAxes); - registry.add("hPartRecoGenNonPromptHelicity", "THn for polarisation studies with cosThStar w.r.t. helicity axis and BDT scores for partially reconstructed generated non-prompt D*+ candidates", HistType::kTHnSparseF, hgenNonPromptAxes); + registry.add("hGenPartRecoPromptHelicity", "THn for polarisation studies with cosThStar w.r.t. helicity axis and BDT scores for partially reconstructed generated prompt D*+ candidates", HistType::kTHnSparseF, hgenPromptAxes); + registry.add("hGenPartRecoNonPromptHelicity", "THn for polarisation studies with cosThStar w.r.t. helicity axis and BDT scores for partially reconstructed generated non-prompt D*+ candidates", HistType::kTHnSparseF, hgenNonPromptAxes); } } } @@ -500,8 +506,8 @@ struct HfTaskCharmPolarisation { registry.add("hGenPromptProduction", "THn for polarisation studies with cosThStar w.r.t. production axis and BDT scores for generated prompt D*+ candidates", HistType::kTHnSparseF, hgenPromptAxes); registry.add("hGenNonPromptProduction", "THn for polarisation studies with cosThStar w.r.t. production axis and BDT scores for generated non-prompt D*+ candidates", HistType::kTHnSparseF, hgenNonPromptAxes); if (activatePartRecoDstar) { - registry.add("hPartRecoGenPromptProduction", "THn for polarisation studies with cosThStar w.r.t. production axis and BDT scores for partially reconstructed generated prompt D*+ candidates", HistType::kTHnSparseF, hgenPromptAxes); - registry.add("hPartRecoGenNonPromptProduction", "THn for polarisation studies with cosThStar w.r.t. production axis and BDT scores for partially reconstructed generated non-prompt D*+ candidates", HistType::kTHnSparseF, hgenNonPromptAxes); + registry.add("hGenPartRecoPromptProduction", "THn for polarisation studies with cosThStar w.r.t. production axis and BDT scores for partially reconstructed generated prompt D*+ candidates", HistType::kTHnSparseF, hgenPromptAxes); + registry.add("hGenPartRecoNonPromptProduction", "THn for polarisation studies with cosThStar w.r.t. production axis and BDT scores for partially reconstructed generated non-prompt D*+ candidates", HistType::kTHnSparseF, hgenNonPromptAxes); } } } @@ -573,8 +579,8 @@ struct HfTaskCharmPolarisation { registry.add("hGenPromptBeam", "THn for polarisation studies with cosThStar w.r.t. beam axis and BDT scores for generated prompt D*+ candidates", HistType::kTHnSparseF, hgenPromptAxes); registry.add("hGenNonPromptBeam", "THn for polarisation studies with cosThStar w.r.t. beam axis and BDT scores for generated non-prompt D*+ candidates", HistType::kTHnSparseF, hgenNonPromptAxes); if (activatePartRecoDstar) { - registry.add("hPartRecoGenPromptBeam", "THn for polarisation studies with cosThStar w.r.t. beam axis and BDT scores for partially reconstructed generated prompt D*+ candidates", HistType::kTHnSparseF, hgenPromptAxes); - registry.add("hPartRecoGenNonPromptBeam", "THn for polarisation studies with cosThStar w.r.t. beam axis and BDT scores for partially reconstructed generated non-prompt D*+ candidates", HistType::kTHnSparseF, hgenNonPromptAxes); + registry.add("hGenPartRecoPromptBeam", "THn for polarisation studies with cosThStar w.r.t. beam axis and BDT scores for partially reconstructed generated prompt D*+ candidates", HistType::kTHnSparseF, hgenPromptAxes); + registry.add("hGenPartRecoNonPromptBeam", "THn for polarisation studies with cosThStar w.r.t. beam axis and BDT scores for partially reconstructed generated non-prompt D*+ candidates", HistType::kTHnSparseF, hgenNonPromptAxes); } } } @@ -631,52 +637,28 @@ struct HfTaskCharmPolarisation { registry.add("hGenPromptRandom", "THn for polarisation studies with cosThStar w.r.t. random axis and BDT scores for generated prompt D*+ candidates", HistType::kTHnSparseF, hgenPromptAxes); registry.add("hGenNonPromptRandom", "THn for polarisation studies with cosThStar w.r.t. random axis and BDT scores for generated non-prompt D*+ candidates", HistType::kTHnSparseF, hgenNonPromptAxes); if (activatePartRecoDstar) { - registry.add("hPartRecoGenPromptRandom", "THn for polarisation studies with cosThStar w.r.t. random axis and BDT scores for partially reconstructed generated prompt D*+ candidates", HistType::kTHnSparseF, hgenPromptAxes); - registry.add("hPartRecoGenNonPromptRandom", "THn for polarisation studies with cosThStar w.r.t. random axis and BDT scores for partially reconstructed generated non-prompt D*+ candidates", HistType::kTHnSparseF, hgenNonPromptAxes); + registry.add("hGenPartRecoPromptRandom", "THn for polarisation studies with cosThStar w.r.t. random axis and BDT scores for partially reconstructed generated prompt D*+ candidates", HistType::kTHnSparseF, hgenPromptAxes); + registry.add("hGenPartRecoNonPromptRandom", "THn for polarisation studies with cosThStar w.r.t. random axis and BDT scores for partially reconstructed generated non-prompt D*+ candidates", HistType::kTHnSparseF, hgenNonPromptAxes); } } } - if (activateTHnSparseCosThStarEP && (!doprocessDstarInPbPb || !doprocessDstarMcInPbPb || !doprocessDstarWithMlInPbPb || !doprocessDstarMcWithMlInPbPb)) { + if (activateTHnSparseCosThStarEP && !(doprocessDstarInPbPb || doprocessDstarWithMlInPbPb)) { LOGP(fatal, "THnSparse with cosThStar w.r.t. event plane axis is not supported for pp analysis, please check the configuration!"); } else if (activateTHnSparseCosThStarEP) { std::vector hEPaxes = {thnAxisInvMass, thnAxisPt, thnAxisNumPvContributors, thnAxisY}; - if (doprocessDstar || doprocessDstarMc || doprocessDstarWithMl || doprocessDstarMcWithMl || doprocessDstarInPbPb || doprocessDstarMcInPbPb || doprocessDstarWithMlInPbPb || doprocessDstarMcWithMlInPbPb) { - + if (doprocessDstarInPbPb || doprocessDstarWithMlInPbPb) { hEPaxes.insert(hEPaxes.end(), {thnAxisInvMassD0, thnAxisCosThetaStarEP}); - if (doprocessDstarWithMl || doprocessDstarMcWithMl || doprocessDstarWithMlInPbPb || doprocessDstarMcWithMlInPbPb) { + if (doprocessDstarWithMlInPbPb) { hEPaxes.insert(hEPaxes.end(), {thnAxisMlBkg, thnAxisMlNonPrompt}); } if (activateTrackingSys) { hEPaxes.insert(hEPaxes.end(), {thnAxisAbsEtaTrackMin, thnAxisNumItsClsMin, thnAxisNumTpcClsMin}); } - if (doprocessDstarMc || doprocessDstarMcWithMl || doprocessDstarMcInPbPb || doprocessDstarMcWithMlInPbPb) { - std::vector hRecoPromptEPAxes(hEPaxes); - hRecoPromptEPAxes.insert(hRecoPromptEPAxes.end(), {thnAxisDauToMuons}); - std::vector hRecoNonPromptEPAxes(hEPaxes); - hRecoNonPromptEPAxes.insert(hRecoNonPromptEPAxes.end(), {thnAxisDauToMuons, thnAxisPtB}); - registry.add("hRecoPromptEP", "THn for polarisation studies with cosThStar w.r.t. event plane axis and BDT scores for reconstructed prompt D*+ candidates", HistType::kTHnSparseF, hRecoPromptEPAxes); - registry.add("hRecoNonPromptEP", "THn for polarisation studies with cosThStar w.r.t. event plane axis and BDT scores for reconstructed non-prompt D*+ candidates", HistType::kTHnSparseF, hRecoNonPromptEPAxes); - if (activatePartRecoDstar) { - registry.add("hPartRecoPromptEP", "THn for polarisation studies with cosThStar w.r.t. event plane axis and BDT scores for partially reconstructed prompt D*+ candidates", HistType::kTHnSparseF, hRecoPromptEPAxes); - registry.add("hPartRecoNonPromptEP", "THn for polarisation studies with cosThStar w.r.t. event plane axis and BDT scores for partially reconstructed non-prompt D*+ candidates", HistType::kTHnSparseF, hRecoNonPromptEPAxes); - } - } else { if (nBkgRotations > 0) { hEPaxes.push_back(thnAxisIsRotatedCandidate); } registry.add("hEP", "THn for polarisation studies with cosThStar w.r.t. event plane axis and BDT scores", HistType::kTHnSparseF, hEPaxes); - } - } - if (doprocessDstarMc || doprocessDstarMcWithMl || doprocessDstarMcInPbPb || doprocessDstarMcWithMlInPbPb || doprocessLcToPKPiMc || doprocessLcToPKPiMcWithMl || doprocessLcToPKPiBackgroundMcWithMl) { - std::vector hgenPromptAxes = {thnAxisPt, thnAxisNumPvContributors, thnAxisY, thnAxisCosThetaStarEP, thnAxisDausAcc, thnAxisResoChannelLc, thnAxisCharge}; - std::vector hgenNonPromptAxes = {thnAxisPt, thnAxisNumPvContributors, thnAxisY, thnAxisCosThetaStarEP, thnAxisPtB, thnAxisDausAcc, thnAxisResoChannelLc, thnAxisCharge}; - registry.add("hGenPromptEP", "THn for polarisation studies with cosThStar w.r.t. event plane axis and BDT scores for generated prompt D*+ candidates", HistType::kTHnSparseF, hgenPromptAxes); - registry.add("hGenNonPromptEP", "THn for polarisation studies with cosThStar w.r.t. event plane axis and BDT scores for generated non-prompt D*+ candidates", HistType::kTHnSparseF, hgenNonPromptAxes); - if (activatePartRecoDstar) { - registry.add("hPartRecoGenPromptEP", "THn for polarisation studies with cosThStar w.r.t. event plane axis and BDT scores for partially reconstructed generated prompt D*+ candidates", HistType::kTHnSparseF, hgenPromptAxes); - registry.add("hPartRecoGenNonPromptEP", "THn for polarisation studies with cosThStar w.r.t. event plane axis and BDT scores for partially reconstructed generated non-prompt D*+ candidates", HistType::kTHnSparseF, hgenNonPromptAxes); - } } } @@ -1492,8 +1474,7 @@ struct HfTaskCharmPolarisation { /// Get the Q vector /// \param collision is the collision with the Q vector information - template - std::vector getQVec(CollsWithQVecs const& collision) + std::vector getQVec(CollsWithQVecs::iterator const& collision) { float xQVec = -999.; float yQVec = -999.; @@ -2068,9 +2049,18 @@ struct HfTaskCharmPolarisation { /// \param mcParticle is the Mc particle /// \param mcParticles is the table of Mc particles /// \param numPvContributors is the number of PV contributors in the associated reco collision - template - void runMcGenPolarisationAnalysis(Part const& mcParticle, Particles const& mcParticles, int numPvContributors) + template + void runMcGenPolarisationAnalysis(Part const& mcParticle, Particles const& mcParticles, int numPvContributors, Cent const* centrality = nullptr) { + if constexpr (withCent) { + assert(qVecs && "Centrality analysis requested but Cent == nullptr"); + } + if constexpr (withCent) { + if (*centrality < centralityMin || *centrality > centralityMax) { + return; // skip this collision if outside of the centrality range + } + } + int8_t origin{RecoDecay::OriginType::None}; std::vector listDaughters{}; float massDau{0.f}, massCharmHad{0.f}; @@ -2314,7 +2304,7 @@ struct HfTaskCharmPolarisation { PROCESS_SWITCH(HfTaskCharmPolarisation, processDstarMcWithMl, "Process Dstar candidates in MC with ML", false); void processDstarInPbPb(CollsWithQVecs::iterator const& collision, - FilteredCandDstarWSelFlagAndMl const& dstarCandidates, + FilteredCandDstarWSelFlag const& dstarCandidates, TracksWithExtra const& tracks) { float centrality = {-1.f}; @@ -2322,6 +2312,8 @@ struct HfTaskCharmPolarisation { if (centrality < centralityMin || centrality > centralityMax) { return; // skip this collision if outside of the centrality range } + registry.fill(HIST("hCentrality"), centrality); + auto thisCollId = collision.globalIndex(); int numPvContributors = collision.numContrib(); auto groupedDstarCandidates = dstarCandidates.sliceBy(dstarPerCollision, thisCollId); @@ -2348,6 +2340,8 @@ struct HfTaskCharmPolarisation { if (centrality < centralityMin || centrality > centralityMax) { return; // skip this collision if outside of the centrality range } + registry.fill(HIST("hCentrality"), centrality); + auto thisCollId = collision.globalIndex(); int numPvContributors = collision.numContrib(); auto groupedDstarCandidates = dstarCandidates.sliceBy(dstarWithMlPerCollision, thisCollId); @@ -2365,22 +2359,22 @@ struct HfTaskCharmPolarisation { } PROCESS_SWITCH(HfTaskCharmPolarisation, processDstarWithMlInPbPb, "Process Dstar candidates with ML in PbPb collisions", false); - void processDstarMcInPbPb(GenCollisWithQvecs::iterator const& collision, + void processDstarMcInPbPb(aod::McCollisions::iterator const&, McParticlesDstarMatched const& mcParticles, - CollsWithQVecsWithMcLabels const& collisions, // this is grouped with SmallGroupsCollisionsWithMcLabels const& collisions, + CollisionsWithMcLabelsAndCent const& collisions, // this is grouped with SmallGroupsCollisionsWithMcLabels const& collisions, FilteredCandDstarWSelFlagAndMc const& dstarCandidates, TracksWithExtra const& tracks) { float centrality = {-1.f}; - centrality = o2::hf_centrality::getCentralityColl(collision, centEstimator); - if (centrality < centralityMin || centrality > centralityMax) { - return; // skip this collision if outside of the centrality range - } int numPvContributorsGen{0}; - std::vector qVecs = getQVec(collision); - for (const auto& collision : collisions) { // loop over reco collisions associated to this gen collision + centrality = o2::hf_centrality::getCentralityColl(collision, centEstimator); + if (centrality < centralityMin || centrality > centralityMax) { + return; // skip this collision if outside of the centrality range + } + registry.fill(HIST("hCentrality"), centrality); + auto thisCollId = collision.globalIndex(); int numPvContributors = collision.numContrib(); auto groupedDstarCandidates = dstarCandidates.sliceBy(dstarWithMcPerCollision, thisCollId); @@ -2392,34 +2386,36 @@ struct HfTaskCharmPolarisation { for (const auto& dstarCandidate : groupedDstarCandidates) { nCands++; - if (runPolarisationAnalysis(dstarCandidate, 0, numPvContributors, -1 /*MC particles*/, tracks, &qVecs)) { + if (runPolarisationAnalysis(dstarCandidate, 0, numPvContributors, -1 /*MC particles*/, tracks)) { nCandsInSignalRegion++; } } fillMultHistos(numPvContributors, nCands, nCandsInSignalRegion); } for (const auto& mcParticle : mcParticles) { - runMcGenPolarisationAnalysis(mcParticle, mcParticles, numPvContributorsGen); + const auto& recoCollsPerMcColl = collisions.sliceBy(colPerMcCollision, mcParticle.mcCollision().globalIndex()); + float cent = o2::hf_centrality::getCentralityGenColl(recoCollsPerMcColl, centEstimator); + runMcGenPolarisationAnalysis(mcParticle, mcParticles, numPvContributorsGen, ¢); } } PROCESS_SWITCH(HfTaskCharmPolarisation, processDstarMcInPbPb, "Process Dstar candidates in PbPb MC without ML", false); - void processDstarMcWithMlInPbPb(GenCollisWithQvecs::iterator const& collision, + void processDstarMcWithMlInPbPb(aod::McCollisions::iterator const&, McParticlesDstarMatched const& mcParticles, - CollsWithQVecsWithMcLabels const& collisions, // this is grouped with SmallGroupsCollisionsWithMcLabels const& collisions, + CollisionsWithMcLabelsAndCent const& collisions, // this is grouped with SmallGroupsCollisionsWithMcLabels const& collisions, FilteredCandDstarWSelFlagAndMcAndMl const& dstarCandidates, TracksWithExtra const& tracks) { float centrality = {-1.f}; - centrality = o2::hf_centrality::getCentralityColl(collision, centEstimator); - if (centrality < centralityMin || centrality > centralityMax) { - return; // skip this collision if outside of the centrality range - } int numPvContributorsGen{0}; - std::vector qVecs = getQVec(collision); - for (const auto& collision : collisions) { // loop over reco collisions associated to this gen collision + centrality = o2::hf_centrality::getCentralityColl(collision, centEstimator); + if (centrality < centralityMin || centrality > centralityMax) { + return; // skip this collision if outside of the centrality range + } + registry.fill(HIST("hCentrality"), centrality); + auto thisCollId = collision.globalIndex(); int numPvContributors = collision.numContrib(); auto groupedDstarCandidates = dstarCandidates.sliceBy(dstarWithMcAndMlPerCollision, thisCollId); @@ -2431,14 +2427,16 @@ struct HfTaskCharmPolarisation { for (const auto& dstarCandidate : groupedDstarCandidates) { nCands++; - if (runPolarisationAnalysis(dstarCandidate, 0, numPvContributors, -1 /*MC particles*/, tracks, &qVecs)) { + if (runPolarisationAnalysis(dstarCandidate, 0, numPvContributors, -1 /*MC particles*/, tracks)) { nCandsInSignalRegion++; } } fillMultHistos(numPvContributors, nCands, nCandsInSignalRegion); } for (const auto& mcParticle : mcParticles) { - runMcGenPolarisationAnalysis(mcParticle, mcParticles, numPvContributorsGen); + const auto& recoCollsPerMcColl = collisions.sliceBy(colPerMcCollision, mcParticle.mcCollision().globalIndex()); + float cent = o2::hf_centrality::getCentralityGenColl(recoCollsPerMcColl, centEstimator); + runMcGenPolarisationAnalysis(mcParticle, mcParticles, numPvContributorsGen, ¢); } } PROCESS_SWITCH(HfTaskCharmPolarisation, processDstarMcWithMlInPbPb, "Process Dstar candidates in PbPb MC with ML", false);