From 5b3dc40e0d335d8d78fcced55b012251653f8876 Mon Sep 17 00:00:00 2001 From: Mingyu Zhang <83645570+Mingyu3360715@users.noreply.github.com> Date: Tue, 7 May 2024 03:08:57 +0800 Subject: [PATCH 1/5] PWGHF: add ptB in THnSparse, add MC gen THnSparse --- PWGHF/D2H/Tasks/taskD0.cxx | 51 ++++++++++++++++++++++++++++++-------- 1 file changed, 40 insertions(+), 11 deletions(-) diff --git a/PWGHF/D2H/Tasks/taskD0.cxx b/PWGHF/D2H/Tasks/taskD0.cxx index 84edf213c89..64f6e9b61de 100644 --- a/PWGHF/D2H/Tasks/taskD0.cxx +++ b/PWGHF/D2H/Tasks/taskD0.cxx @@ -56,10 +56,13 @@ struct HfTaskD0 { ConfigurableAxis thnConfigAxisBkgScore{"thnConfigAxisBkgScore", {50, 0, 1}, "Bkg score bins"}; ConfigurableAxis thnConfigAxisNonPromptScore{"thnConfigAxisNonPromptScore", {50, 0, 1}, "Non-prompt score bins"}; ConfigurableAxis thnConfigAxisMass{"thnConfigAxisMass", {120, 1.5848, 2.1848}, "Cand. inv-mass bins"}; - ConfigurableAxis thnConfigAxisPt{"thnConfigAxisPt", {120, 0, 24}, "Cand. pT bins"}; + ConfigurableAxis thnConfigAxisPtB{"thnConfigAxisPtB", {1000, 0, 100}, "Cand. beauty mother pTB bins"}; + ConfigurableAxis thnConfigAxisPt{"thnConfigAxisPt", {500, 0, 50}, "Cand. pT bins"}; ConfigurableAxis thnConfigAxisY{"thnConfigAxisY", {20, -1, 1}, "Cand. rapidity bins"}; ConfigurableAxis thnConfigAxisOrigin{"thnConfigAxisOrigin", {3, -0.5, 2.5}, "Cand. origin type"}; ConfigurableAxis thnConfigAxisCandType{"thnConfigAxisCandType", {4, -0.5, 3.5}, "D0 type"}; + ConfigurableAxis thnConfigAxisGenPtD{"thnConfigAxisGenPtD", {500, 0, 50}, "Gen Pt D"}; + ConfigurableAxis thnConfigAxisGenPtB{"thnConfigAxisGenPtB", {1000, 0, 100}, "Gen Pt B"}; HfHelper hfHelper; @@ -213,11 +216,17 @@ struct HfTaskD0 { const AxisSpec thnAxisNonPromptScore{thnConfigAxisNonPromptScore, "BDT score non-prompt."}; const AxisSpec thnAxisMass{thnConfigAxisMass, "inv. mass (#pi K) (GeV/#it{c}^{2})"}; const AxisSpec thnAxisPt{thnConfigAxisPt, "#it{p}_{T} (GeV/#it{c})"}; + const AxisSpec thnAxisPtB{thnConfigAxisPtB, "#it{p}_{T}^{B} (GeV/#it{c})"}; const AxisSpec thnAxisY{thnConfigAxisY, "y"}; const AxisSpec thnAxisOrigin{thnConfigAxisOrigin, "Origin"}; const AxisSpec thnAxisCandType{thnConfigAxisCandType, "D0 type"}; + const AxisSpec thnAxisGenPtD{thnConfigAxisGenPtD, "#it{p}_{T} (GeV/#it{c})"}; + const AxisSpec thnAxisGenPtB{thnConfigAxisGenPtB, "#it{p}_{T}^{B} (GeV/#it{c})"}; - registry.add("hBdtScoreVsMassVsPtVsYVsOriginVsD0Type", "Thn for D0 candidates with BDT scores", HistType::kTHnSparseD, {thnAxisBkgScore, thnAxisNonPromptScore, thnAxisMass, thnAxisPt, thnAxisY, thnAxisOrigin, thnAxisCandType}); + registry.add("hSparseAcc", "Thn for generated D0 from charm and beauty", HistType::kTHnSparseD, {thnAxisGenPtD, thnAxisGenPtB, thnAxisY, thnAxisOrigin}); + registry.get(HIST("hSparseAcc"))->Sumw2(); + + registry.add("hBdtScoreVsMassVsPtVsYVsOriginVsD0Type", "Thn for D0 candidates with BDT scores", HistType::kTHnSparseD, {thnAxisBkgScore, thnAxisNonPromptScore, thnAxisMass, thnAxisPt, thnAxisPtB, thnAxisY, thnAxisOrigin, thnAxisCandType}); registry.get(HIST("hBdtScoreVsMassVsPtVsYVsOriginVsD0Type"))->Sumw2(); } } @@ -284,10 +293,10 @@ struct HfTaskD0 { if constexpr (applyMl) { if (candidate.isSelD0() >= selectionFlagD0) { - registry.fill(HIST("hBdtScoreVsMassVsPtVsYVsOriginVsD0Type"), candidate.mlProbD0()[0], candidate.mlProbD0()[1], massD0, ptCandidate, hfHelper.yD0(candidate), 0, SigD0); + registry.fill(HIST("hBdtScoreVsMassVsPtVsYVsOriginVsD0Type"), candidate.mlProbD0()[0], candidate.mlProbD0()[1], massD0, ptCandidate, -1, hfHelper.yD0(candidate), 0, SigD0); } if (candidate.isSelD0bar() >= selectionFlagD0bar) { - registry.fill(HIST("hBdtScoreVsMassVsPtVsYVsOriginVsD0Type"), candidate.mlProbD0bar()[0], candidate.mlProbD0bar()[1], massD0bar, ptCandidate, hfHelper.yD0(candidate), 0, SigD0bar); + registry.fill(HIST("hBdtScoreVsMassVsPtVsYVsOriginVsD0Type"), candidate.mlProbD0bar()[0], candidate.mlProbD0bar()[1], massD0bar, ptCandidate, -1, hfHelper.yD0(candidate), 0, SigD0bar); } } } @@ -341,9 +350,9 @@ struct HfTaskD0 { // Get the corresponding MC particle. auto indexMother = RecoDecay::getMother(mcParticles, candidate.template prong0_as().template mcParticle_as>(), o2::constants::physics::Pdg::kD0, true); auto particleMother = mcParticles.rawIteratorAt(indexMother); - auto ptGen = particleMother.pt(); // gen. level pT - auto yGen = RecoDecay::y(particleMother.pVector(), o2::constants::physics::MassD0); // gen. level y - registry.fill(HIST("hPtGenSig"), ptGen); // gen. level pT + auto ptGen = particleMother.pt(); // gen. level pT + auto yGen = RecoDecay::y(particleMother.pVector(), o2::constants::physics::MassD0); // gen. level y + registry.fill(HIST("hPtGenSig"), ptGen); // gen. level pT auto ptRec = candidate.pt(); auto yRec = hfHelper.yD0(candidate); if (candidate.isRecoHfFlag() >= selectionFlagHf) { @@ -427,6 +436,9 @@ struct HfTaskD0 { auto ctCandidate = hfHelper.ctD0(candidate); auto cpaCandidate = candidate.cpa(); auto cpaxyCandidate = candidate.cpaXY(); + float ptB = -1; + auto indexMother = RecoDecay::getMother(mcParticles, candidate.template prong0_as().template mcParticle_as>(), o2::constants::physics::Pdg::kD0, true); + auto particleMother = mcParticles.rawIteratorAt(indexMother); if (candidate.isSelD0() >= selectionFlagD0) { registry.fill(HIST("hMassSigBkgD0"), massD0, ptCandidate, rapidityCandidate); if (candidate.flagMcMatchRec() == (1 << aod::hf_cand_2prong::DecayType::D0ToPiK)) { @@ -454,7 +466,10 @@ struct HfTaskD0 { registry.fill(HIST("hDecLengthxyVsPtSig"), declengthxyCandidate, ptCandidate); registry.fill(HIST("hMassSigD0"), massD0, ptCandidate, rapidityCandidate); if constexpr (applyMl) { - registry.fill(HIST("hBdtScoreVsMassVsPtVsYVsOriginVsD0Type"), candidate.mlProbD0()[0], candidate.mlProbD0()[1], massD0, ptCandidate, rapidityCandidate, candidate.originMcRec(), SigD0); + if (candidate.originMcRec() == RecoDecay::OriginType::NonPrompt) { + ptB = hfHelper.getBeautyMotherPt(mcParticles, particleMother); + } + registry.fill(HIST("hBdtScoreVsMassVsPtVsYVsOriginVsD0Type"), candidate.mlProbD0()[0], candidate.mlProbD0()[1], massD0, ptCandidate, ptB, rapidityCandidate, candidate.originMcRec(), SigD0); } } else { registry.fill(HIST("hPtProng0Bkg"), ptProng0, rapidityCandidate); @@ -474,7 +489,10 @@ struct HfTaskD0 { if (candidate.flagMcMatchRec() == -(1 << aod::hf_cand_2prong::DecayType::D0ToPiK)) { registry.fill(HIST("hMassReflBkgD0"), massD0, ptCandidate, rapidityCandidate); if constexpr (applyMl) { - registry.fill(HIST("hBdtScoreVsMassVsPtVsYVsOriginVsD0Type"), candidate.mlProbD0()[0], candidate.mlProbD0()[1], massD0, ptCandidate, rapidityCandidate, candidate.originMcRec(), ReflectedD0); + if (candidate.originMcRec() == RecoDecay::OriginType::NonPrompt) { + ptB = hfHelper.getBeautyMotherPt(mcParticles, particleMother); + } + registry.fill(HIST("hBdtScoreVsMassVsPtVsYVsOriginVsD0Type"), candidate.mlProbD0()[0], candidate.mlProbD0()[1], massD0, ptCandidate, ptB, rapidityCandidate, candidate.originMcRec(), ReflectedD0); } } } @@ -484,14 +502,20 @@ struct HfTaskD0 { if (candidate.flagMcMatchRec() == -(1 << aod::hf_cand_2prong::DecayType::D0ToPiK)) { registry.fill(HIST("hMassSigD0bar"), massD0bar, ptCandidate, rapidityCandidate); if constexpr (applyMl) { - registry.fill(HIST("hBdtScoreVsMassVsPtVsYVsOriginVsD0Type"), candidate.mlProbD0bar()[0], candidate.mlProbD0bar()[1], massD0bar, ptCandidate, rapidityCandidate, candidate.originMcRec(), SigD0bar); + if (candidate.originMcRec() == RecoDecay::OriginType::NonPrompt) { + ptB = hfHelper.getBeautyMotherPt(mcParticles, particleMother); + } + registry.fill(HIST("hBdtScoreVsMassVsPtVsYVsOriginVsD0Type"), candidate.mlProbD0bar()[0], candidate.mlProbD0bar()[1], massD0bar, ptCandidate, ptB, rapidityCandidate, candidate.originMcRec(), SigD0bar); } } else { registry.fill(HIST("hMassBkgD0bar"), massD0bar, ptCandidate, rapidityCandidate); if (candidate.flagMcMatchRec() == (1 << aod::hf_cand_2prong::DecayType::D0ToPiK)) { registry.fill(HIST("hMassReflBkgD0bar"), massD0bar, ptCandidate, rapidityCandidate); if constexpr (applyMl) { - registry.fill(HIST("hBdtScoreVsMassVsPtVsYVsOriginVsD0Type"), candidate.mlProbD0bar()[0], candidate.mlProbD0bar()[1], massD0bar, ptCandidate, rapidityCandidate, candidate.originMcRec(), ReflectedD0bar); + if (candidate.originMcRec() == RecoDecay::OriginType::NonPrompt) { + ptB = hfHelper.getBeautyMotherPt(mcParticles, particleMother); + } + registry.fill(HIST("hBdtScoreVsMassVsPtVsYVsOriginVsD0Type"), candidate.mlProbD0bar()[0], candidate.mlProbD0bar()[1], massD0bar, ptCandidate, ptB, rapidityCandidate, candidate.originMcRec(), ReflectedD0bar); } } } @@ -503,18 +527,23 @@ struct HfTaskD0 { if (yCandGenMax >= 0. && std::abs(RecoDecay::y(particle.pVector(), o2::constants::physics::MassD0)) > yCandGenMax) { continue; } + float ptGenB = -1; auto ptGen = particle.pt(); auto yGen = RecoDecay::y(particle.pVector(), o2::constants::physics::MassD0); registry.fill(HIST("hPtGen"), ptGen); registry.fill(HIST("hPtVsYGen"), ptGen, yGen); + if (particle.originMcGen() == RecoDecay::OriginType::Prompt) { registry.fill(HIST("hPtGenPrompt"), ptGen); registry.fill(HIST("hYGenPrompt"), yGen); registry.fill(HIST("hPtVsYGenPrompt"), ptGen, yGen); + registry.fill(HIST("hSparseAcc"), ptGen, ptGenB, yGen, 1); } else { + ptGenB = hfHelper.getBeautyMotherPt(mcParticles, particle); registry.fill(HIST("hPtGenNonPrompt"), ptGen); registry.fill(HIST("hYGenNonPrompt"), yGen); registry.fill(HIST("hPtVsYGenNonPrompt"), ptGen, yGen); + registry.fill(HIST("hSparseAcc"), ptGen, ptGenB, yGen, 2); } registry.fill(HIST("hEtaGen"), particle.eta()); } From 3221467b999ba133bd0356cfc25261c0eb508acd Mon Sep 17 00:00:00 2001 From: Mingyu Zhang <83645570+Mingyu3360715@users.noreply.github.com> Date: Tue, 7 May 2024 03:11:44 +0800 Subject: [PATCH 2/5] PWGHF: add function for getting pt of D meson beauty mother --- PWGHF/Core/HfHelper.h | 51 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 51 insertions(+) diff --git a/PWGHF/Core/HfHelper.h b/PWGHF/Core/HfHelper.h index 2ea6afa54ad..0abd41a70e3 100644 --- a/PWGHF/Core/HfHelper.h +++ b/PWGHF/Core/HfHelper.h @@ -931,6 +931,57 @@ class HfHelper return true; } + /// Get the pt of the beauty hadron(non-prompt case), return negative value for prompt case + /// \param particlesMC table with MC particles + /// \param particle MC particle + /// \return a double corresponding to beauty mother pt + template + double getBeautyMotherPt(const T& particlesMC, + const typename T::iterator& particle) + { + int stage = 0; // mother tree level + + // vector of vectors with mother indices; each line corresponds to a "stage" + std::vector> arrayIds{}; + std::vector initVec{particle.globalIndex()}; + arrayIds.push_back(initVec); // the first vector contains the index of the original particle + + while (arrayIds[-stage].size() > 0) { + // vector of mother indices for the current stage + std::vector arrayIdsStage{}; + for (auto& iPart : arrayIds[-stage]) { // check all the particles that were the mothers at the previous stage + auto particleMother = particlesMC.rawIteratorAt(iPart - particlesMC.offset()); + if (particleMother.has_mothers()) { + for (auto iMother = particleMother.mothersIds().front(); iMother <= particleMother.mothersIds().back(); ++iMother) { // loop over the mother particles of the analysed particle + if (std::find(arrayIdsStage.begin(), arrayIdsStage.end(), iMother) != arrayIdsStage.end()) { // if a mother is still present in the vector, do not check it again + continue; + } + auto mother = particlesMC.rawIteratorAt(iMother - particlesMC.offset()); + // Check mother's PDG code. + auto PDGParticleIMother = std::abs(mother.pdgCode()); // PDG code of the mother + + if ((PDGParticleIMother > 500 && PDGParticleIMother < 600) || (PDGParticleIMother > 5000 && PDGParticleIMother < 6000)) { + return mother.pt(); + } + if (PDGParticleIMother == 4) { + return -999.; + } + if (PDGParticleIMother == 5) { + return -1.; + } + + // add mother index in the vector for the current stage + arrayIdsStage.push_back(iMother); + } + } + } + // add vector of mother indices for the current stage + arrayIds.push_back(arrayIdsStage); + stage--; + } + return -999.; + } + private: }; From 42fe643c64b665abd655486b917621ac1cce4ced Mon Sep 17 00:00:00 2001 From: Mingyu Zhang <83645570+Mingyu3360715@users.noreply.github.com> Date: Tue, 7 May 2024 03:15:42 +0800 Subject: [PATCH 3/5] fix thn name --- PWGHF/D2H/Tasks/taskD0.cxx | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/PWGHF/D2H/Tasks/taskD0.cxx b/PWGHF/D2H/Tasks/taskD0.cxx index 64f6e9b61de..a33b620465b 100644 --- a/PWGHF/D2H/Tasks/taskD0.cxx +++ b/PWGHF/D2H/Tasks/taskD0.cxx @@ -226,8 +226,8 @@ struct HfTaskD0 { registry.add("hSparseAcc", "Thn for generated D0 from charm and beauty", HistType::kTHnSparseD, {thnAxisGenPtD, thnAxisGenPtB, thnAxisY, thnAxisOrigin}); registry.get(HIST("hSparseAcc"))->Sumw2(); - registry.add("hBdtScoreVsMassVsPtVsYVsOriginVsD0Type", "Thn for D0 candidates with BDT scores", HistType::kTHnSparseD, {thnAxisBkgScore, thnAxisNonPromptScore, thnAxisMass, thnAxisPt, thnAxisPtB, thnAxisY, thnAxisOrigin, thnAxisCandType}); - registry.get(HIST("hBdtScoreVsMassVsPtVsYVsOriginVsD0Type"))->Sumw2(); + registry.add("hBdtScoreVsMassVsPtVsPtBVsYVsOriginVsD0Type", "Thn for D0 candidates with BDT scores", HistType::kTHnSparseD, {thnAxisBkgScore, thnAxisNonPromptScore, thnAxisMass, thnAxisPt, thnAxisPtB, thnAxisY, thnAxisOrigin, thnAxisCandType}); + registry.get(HIST("hBdtScoreVsMassVsPtVsPtBVsYVsOriginVsD0Type"))->Sumw2(); } } @@ -293,10 +293,10 @@ struct HfTaskD0 { if constexpr (applyMl) { if (candidate.isSelD0() >= selectionFlagD0) { - registry.fill(HIST("hBdtScoreVsMassVsPtVsYVsOriginVsD0Type"), candidate.mlProbD0()[0], candidate.mlProbD0()[1], massD0, ptCandidate, -1, hfHelper.yD0(candidate), 0, SigD0); + registry.fill(HIST("hBdtScoreVsMassVsPtVsPtBVsYVsOriginVsD0Type"), candidate.mlProbD0()[0], candidate.mlProbD0()[1], massD0, ptCandidate, -1, hfHelper.yD0(candidate), 0, SigD0); } if (candidate.isSelD0bar() >= selectionFlagD0bar) { - registry.fill(HIST("hBdtScoreVsMassVsPtVsYVsOriginVsD0Type"), candidate.mlProbD0bar()[0], candidate.mlProbD0bar()[1], massD0bar, ptCandidate, -1, hfHelper.yD0(candidate), 0, SigD0bar); + registry.fill(HIST("hBdtScoreVsMassVsPtVsPtBVsYVsOriginVsD0Type"), candidate.mlProbD0bar()[0], candidate.mlProbD0bar()[1], massD0bar, ptCandidate, -1, hfHelper.yD0(candidate), 0, SigD0bar); } } } @@ -469,7 +469,7 @@ struct HfTaskD0 { if (candidate.originMcRec() == RecoDecay::OriginType::NonPrompt) { ptB = hfHelper.getBeautyMotherPt(mcParticles, particleMother); } - registry.fill(HIST("hBdtScoreVsMassVsPtVsYVsOriginVsD0Type"), candidate.mlProbD0()[0], candidate.mlProbD0()[1], massD0, ptCandidate, ptB, rapidityCandidate, candidate.originMcRec(), SigD0); + registry.fill(HIST("hBdtScoreVsMassVsPtVsPtBVsYVsOriginVsD0Type"), candidate.mlProbD0()[0], candidate.mlProbD0()[1], massD0, ptCandidate, ptB, rapidityCandidate, candidate.originMcRec(), SigD0); } } else { registry.fill(HIST("hPtProng0Bkg"), ptProng0, rapidityCandidate); @@ -492,7 +492,7 @@ struct HfTaskD0 { if (candidate.originMcRec() == RecoDecay::OriginType::NonPrompt) { ptB = hfHelper.getBeautyMotherPt(mcParticles, particleMother); } - registry.fill(HIST("hBdtScoreVsMassVsPtVsYVsOriginVsD0Type"), candidate.mlProbD0()[0], candidate.mlProbD0()[1], massD0, ptCandidate, ptB, rapidityCandidate, candidate.originMcRec(), ReflectedD0); + registry.fill(HIST("hBdtScoreVsMassVsPtVsPtBVsYVsOriginVsD0Type"), candidate.mlProbD0()[0], candidate.mlProbD0()[1], massD0, ptCandidate, ptB, rapidityCandidate, candidate.originMcRec(), ReflectedD0); } } } @@ -505,7 +505,7 @@ struct HfTaskD0 { if (candidate.originMcRec() == RecoDecay::OriginType::NonPrompt) { ptB = hfHelper.getBeautyMotherPt(mcParticles, particleMother); } - registry.fill(HIST("hBdtScoreVsMassVsPtVsYVsOriginVsD0Type"), candidate.mlProbD0bar()[0], candidate.mlProbD0bar()[1], massD0bar, ptCandidate, ptB, rapidityCandidate, candidate.originMcRec(), SigD0bar); + registry.fill(HIST("hBdtScoreVsMassVsPtVsPtBVsYVsOriginVsD0Type"), candidate.mlProbD0bar()[0], candidate.mlProbD0bar()[1], massD0bar, ptCandidate, ptB, rapidityCandidate, candidate.originMcRec(), SigD0bar); } } else { registry.fill(HIST("hMassBkgD0bar"), massD0bar, ptCandidate, rapidityCandidate); @@ -515,7 +515,7 @@ struct HfTaskD0 { if (candidate.originMcRec() == RecoDecay::OriginType::NonPrompt) { ptB = hfHelper.getBeautyMotherPt(mcParticles, particleMother); } - registry.fill(HIST("hBdtScoreVsMassVsPtVsYVsOriginVsD0Type"), candidate.mlProbD0bar()[0], candidate.mlProbD0bar()[1], massD0bar, ptCandidate, ptB, rapidityCandidate, candidate.originMcRec(), ReflectedD0bar); + registry.fill(HIST("hBdtScoreVsMassVsPtVsPtBVsYVsOriginVsD0Type"), candidate.mlProbD0bar()[0], candidate.mlProbD0bar()[1], massD0bar, ptCandidate, ptB, rapidityCandidate, candidate.originMcRec(), ReflectedD0bar); } } } From 2d650e654f13611569b13fa13803883fbe0a7295 Mon Sep 17 00:00:00 2001 From: Mingyu Zhang <83645570+Mingyu3360715@users.noreply.github.com> Date: Tue, 7 May 2024 03:22:46 +0800 Subject: [PATCH 4/5] add missing #include --- PWGHF/Core/HfHelper.h | 1 + 1 file changed, 1 insertion(+) diff --git a/PWGHF/Core/HfHelper.h b/PWGHF/Core/HfHelper.h index 0abd41a70e3..01ba9a8b165 100644 --- a/PWGHF/Core/HfHelper.h +++ b/PWGHF/Core/HfHelper.h @@ -20,6 +20,7 @@ #include #include #include +#include // std::vector #include "CommonConstants/PhysicsConstants.h" From 7df2960102023d800c95cefb15065b01a649cfce Mon Sep 17 00:00:00 2001 From: Mingyu Zhang <83645570+Mingyu3360715@users.noreply.github.com> Date: Tue, 7 May 2024 03:24:08 +0800 Subject: [PATCH 5/5] fix format issue --- PWGHF/Core/HfHelper.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PWGHF/Core/HfHelper.h b/PWGHF/Core/HfHelper.h index 01ba9a8b165..4808d9d062d 100644 --- a/PWGHF/Core/HfHelper.h +++ b/PWGHF/Core/HfHelper.h @@ -20,7 +20,7 @@ #include #include #include -#include // std::vector +#include #include "CommonConstants/PhysicsConstants.h"