diff --git a/Common/Core/RecoDecay.h b/Common/Core/RecoDecay.h index 8186128329d..1bb5e2fb2c1 100644 --- a/Common/Core/RecoDecay.h +++ b/Common/Core/RecoDecay.h @@ -49,6 +49,11 @@ class RecoDecay /// Default destructor ~RecoDecay() = default; + // mapping of charm-hadron origin type + enum OriginType { None = 0, + Prompt, + NonPrompt }; + // Auxiliary functions /// Sums numbers. @@ -913,6 +918,71 @@ class RecoDecay return true; } + /// Finds the origin (from charm hadronisation or beauty-hadron decay) of charm hadrons. + /// \param particlesMC table with MC particles + /// \param particle MC particle + /// \param searchUpToQuark if true tag origin based on charm/beauty quark otherwise on b-hadron + /// \return an integer corresponding to the origin (0: none, 1: prompt, 2: nonprompt) as in OriginType + template + static int getCharmHadronOrigin(const T& particlesMC, + const typename T::iterator& particle, + const bool searchUpToQuark = false) + { + int stage = 0; // mother tree level (just for debugging) + + // 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 + // printf("getMother: "); + // for (int i = stage; i < 0; i++) // Indent to make the tree look nice. + // printf(" "); + // printf("Stage %d: Mother PDG: %d, Index: %d\n", stage, PDGParticleIMother, iMother); + + if (searchUpToQuark) { + if (PDGParticleIMother == 5) { // b quark + return OriginType::NonPrompt; + } + if (PDGParticleIMother == 4) { // c quark + return OriginType::Prompt; + } + } else { + if ( + (PDGParticleIMother / 100 == 5 || // b mesons + PDGParticleIMother / 1000 == 5) // b baryons + ) { + return OriginType::NonPrompt; + } + } + // 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--; + } + if (!searchUpToQuark) { + return OriginType::Prompt; + } + return OriginType::None; + } + private: static std::vector> mListMass; ///< list of particle masses in form (PDG code, mass) }; diff --git a/EventFiltering/PWGHF/HFFilter.cxx b/EventFiltering/PWGHF/HFFilter.cxx index 23763edfd8f..51213c80f87 100644 --- a/EventFiltering/PWGHF/HFFilter.cxx +++ b/EventFiltering/PWGHF/HFFilter.cxx @@ -39,7 +39,6 @@ using namespace o2; using namespace o2::framework; using namespace o2::framework::expressions; -using namespace o2::aod::hf_cand; using namespace hf_cuts_single_track; using namespace hf_cuts_bdt_multiclass; @@ -982,10 +981,10 @@ struct HfFilter { // Main struct for HF triggers auto indexRec = RecoDecay::getMatchedMCRec(particlesMC, std::array{trackPos, trackNeg}, pdg::Code::kD0, array{+kPiPlus, -kKPlus}, true, &sign); if (indexRec > -1) { auto particle = particlesMC.rawIteratorAt(indexRec); - origin = (RecoDecay::getMother(particlesMC, particle, kBottom, true) > -1 ? OriginType::NonPrompt : OriginType::Prompt); - if (origin == OriginType::NonPrompt) { + origin = RecoDecay::getCharmHadronOrigin(particlesMC, particle); + if (origin == RecoDecay::OriginType::NonPrompt) { flag = kNonPrompt; - } else { + } else if (origin == RecoDecay::OriginType::Prompt) { flag = kPrompt; } } else { @@ -1041,8 +1040,8 @@ struct HfFilter { // Main struct for HF triggers if (indexRec > -1) { auto particle = particlesMC.rawIteratorAt(indexRec); - origin = (RecoDecay::getMother(particlesMC, particle, kBottom, true) > -1 ? OriginType::NonPrompt : OriginType::Prompt); - if (origin == OriginType::NonPrompt) { + origin = RecoDecay::getCharmHadronOrigin(particlesMC, particle); + if (origin == RecoDecay::OriginType::NonPrompt) { flag = kNonPrompt; } else { flag = kPrompt; diff --git a/PWGHF/DataModel/HFSecondaryVertex.h b/PWGHF/DataModel/HFSecondaryVertex.h index 0a4648eae61..27b10083359 100644 --- a/PWGHF/DataModel/HFSecondaryVertex.h +++ b/PWGHF/DataModel/HFSecondaryVertex.h @@ -290,10 +290,6 @@ DECLARE_SOA_DYNAMIC_COLUMN(Ct, ct, //! [](float xVtxP, float yVtxP, float zVtxP, float xVtxS, float yVtxS, float zVtxS, float px, float py, float pz, double m) -> float { return RecoDecay::ct(array{px, py, pz}, RecoDecay::distance(array{xVtxP, yVtxP, zVtxP}, array{xVtxS, yVtxS, zVtxS}), m); }); DECLARE_SOA_DYNAMIC_COLUMN(ImpactParameterXY, impactParameterXY, //! [](float xVtxP, float yVtxP, float zVtxP, float xVtxS, float yVtxS, float zVtxS, float px, float py, float pz) -> float { return RecoDecay::impParXY(array{xVtxP, yVtxP, zVtxP}, array{xVtxS, yVtxS, zVtxS}, array{px, py, pz}); }); - -// mapping of origin type -enum OriginType { Prompt = 1, - NonPrompt }; } // namespace hf_cand // specific 2-prong decay properties diff --git a/PWGHF/TableProducer/HFCandidateCreator2Prong.cxx b/PWGHF/TableProducer/HFCandidateCreator2Prong.cxx index d5b459020ac..b0fc48a20bd 100644 --- a/PWGHF/TableProducer/HFCandidateCreator2Prong.cxx +++ b/PWGHF/TableProducer/HFCandidateCreator2Prong.cxx @@ -210,7 +210,7 @@ struct HFCandidateCreator2ProngExpressions { // Check whether the particle is non-prompt (from a b quark). if (flag != 0) { auto particle = particlesMC.rawIteratorAt(indexRec); - origin = (RecoDecay::getMother(particlesMC, particle, kBottom, true) > -1 ? OriginType::NonPrompt : OriginType::Prompt); + origin = RecoDecay::getCharmHadronOrigin(particlesMC, particle); } rowMCMatchRec(flag, origin); @@ -246,7 +246,7 @@ struct HFCandidateCreator2ProngExpressions { // Check whether the particle is non-prompt (from a b quark). if (flag != 0) { - origin = (RecoDecay::getMother(particlesMC, particle, kBottom, true) > -1 ? OriginType::NonPrompt : OriginType::Prompt); + origin = RecoDecay::getCharmHadronOrigin(particlesMC, particle); } rowMCMatchGen(flag, origin); diff --git a/PWGHF/TableProducer/HFCandidateCreator3Prong.cxx b/PWGHF/TableProducer/HFCandidateCreator3Prong.cxx index c48d70c8abb..61853729578 100644 --- a/PWGHF/TableProducer/HFCandidateCreator3Prong.cxx +++ b/PWGHF/TableProducer/HFCandidateCreator3Prong.cxx @@ -251,7 +251,7 @@ struct HFCandidateCreator3ProngExpressions { // Check whether the particle is non-prompt (from a b quark). if (flag != 0) { auto particle = particlesMC.rawIteratorAt(indexRec); - origin = (RecoDecay::getMother(particlesMC, particle, kBottom, true) > -1 ? OriginType::NonPrompt : OriginType::Prompt); + origin = RecoDecay::getCharmHadronOrigin(particlesMC, particle); } rowMCMatchRec(flag, origin, swapping, channel); @@ -305,7 +305,7 @@ struct HFCandidateCreator3ProngExpressions { // Check whether the particle is non-prompt (from a b quark). if (flag != 0) { - origin = (RecoDecay::getMother(particlesMC, particle, kBottom, true) > -1 ? OriginType::NonPrompt : OriginType::Prompt); + origin = RecoDecay::getCharmHadronOrigin(particlesMC, particle); } rowMCMatchGen(flag, origin, channel); diff --git a/PWGHF/TableProducer/HFCandidateCreatorChic.cxx b/PWGHF/TableProducer/HFCandidateCreatorChic.cxx index 1173eae2cde..07bfe5f0f89 100644 --- a/PWGHF/TableProducer/HFCandidateCreatorChic.cxx +++ b/PWGHF/TableProducer/HFCandidateCreatorChic.cxx @@ -249,7 +249,7 @@ struct HFCandidateCreatorChicMC { } if (flag != 0) { auto particle = particlesMC.rawIteratorAt(indexRec); - origin = (RecoDecay::getMother(particlesMC, particle, kBottom, true) > -1 ? NonPrompt : Prompt); + origin = RecoDecay::getCharmHadronOrigin(particlesMC, particle); } rowMCMatchRec(flag, origin, channel); } diff --git a/PWGHF/TableProducer/HFCandidateCreatorX.cxx b/PWGHF/TableProducer/HFCandidateCreatorX.cxx index 159ffa34084..96f37e8b619 100644 --- a/PWGHF/TableProducer/HFCandidateCreatorX.cxx +++ b/PWGHF/TableProducer/HFCandidateCreatorX.cxx @@ -305,7 +305,7 @@ struct HFCandidateCreatorXMC { // Check whether the particle is non-prompt (from a b quark). if (flag != 0) { auto particle = particlesMC.rawIteratorAt(indexRec); - origin = (RecoDecay::getMother(particlesMC, particle, 5, true) > -1 ? NonPrompt : Prompt); + origin = RecoDecay::getCharmHadronOrigin(particlesMC, particle); } rowMCMatchRec(flag, origin, channel); @@ -338,7 +338,7 @@ struct HFCandidateCreatorXMC { // Check whether the particle is non-prompt (from a b quark). if (flag != 0) { - origin = (RecoDecay::getMother(particlesMC, particle, 5, true) > -1 ? NonPrompt : Prompt); + origin = RecoDecay::getCharmHadronOrigin(particlesMC, particle); } rowMCMatchGen(flag, origin, channel); diff --git a/PWGHF/Tasks/HFMCValidation.cxx b/PWGHF/Tasks/HFMCValidation.cxx index eddad85c5dd..99cb7d73149 100644 --- a/PWGHF/Tasks/HFMCValidation.cxx +++ b/PWGHF/Tasks/HFMCValidation.cxx @@ -26,7 +26,6 @@ using namespace o2; using namespace o2::framework; using namespace o2::framework::expressions; using namespace o2::aod; -using namespace o2::aod::hf_cand; using namespace o2::aod::hf_cand_prong2; using namespace o2::aod::hf_cand_prong3; @@ -157,10 +156,10 @@ struct ValidationGenLevel { std::size_t arrayPDGsize = arrPDGFinal[iD].size() - std::count(arrPDGFinal[iD].begin(), arrPDGFinal[iD].end(), 0); int origin = -1; if (listDaughters.size() == arrayPDGsize) { - origin = (RecoDecay::getMother(particlesMC, particle, kBottom, true) > -1 ? OriginType::NonPrompt : OriginType::Prompt); - if (origin == OriginType::Prompt) { + origin = RecoDecay::getCharmHadronOrigin(particlesMC, particle); + if (origin == RecoDecay::OriginType::Prompt) { counterPrompt[iD]++; - } else if (origin == OriginType::NonPrompt) { + } else if (origin == RecoDecay::OriginType::NonPrompt) { counterNonPrompt[iD]++; } } @@ -190,13 +189,13 @@ struct ValidationGenLevel { registry.fill(HIST("hPzDiffMotherDaughterGen"), pzDiff); registry.fill(HIST("hPDiffMotherDaughterGen"), pDiff); registry.fill(HIST("hPtDiffMotherDaughterGen"), ptDiff); - if (origin == OriginType::Prompt) { + if (origin == RecoDecay::OriginType::Prompt) { if (std::abs(particle.y()) < 0.5) { hPromptCharmHadronsPtDistr->Fill(whichHadron, particle.pt()); } hPromptCharmHadronsYDistr->Fill(whichHadron, particle.y()); hPromptCharmHadronsDecLenDistr->Fill(whichHadron, decayLength * 10000); - } else if (origin == OriginType::NonPrompt) { + } else if (origin == RecoDecay::OriginType::NonPrompt) { if (std::abs(particle.y()) < 0.5) { hNonPromptCharmHadronsPtDistr->Fill(whichHadron, particle.pt()); } @@ -288,7 +287,7 @@ struct ValidationRecLevel { whichHad = 6; } int whichOrigin = -1; - if (cand2Prong.originMCRec() == OriginType::Prompt) { + if (cand2Prong.originMCRec() == RecoDecay::OriginType::Prompt) { whichOrigin = 0; } else { whichOrigin = 1; @@ -353,7 +352,7 @@ struct ValidationRecLevel { whichHad = 5; } int whichOrigin = -1; - if (cand3Prong.originMCRec() == OriginType::Prompt) { + if (cand3Prong.originMCRec() == RecoDecay::OriginType::Prompt) { whichOrigin = 0; } else { whichOrigin = 1; diff --git a/PWGHF/Tasks/HFSelOptimisation.cxx b/PWGHF/Tasks/HFSelOptimisation.cxx index da48bfdd3e6..188418871cb 100644 --- a/PWGHF/Tasks/HFSelOptimisation.cxx +++ b/PWGHF/Tasks/HFSelOptimisation.cxx @@ -23,7 +23,6 @@ using namespace o2; using namespace o2::framework; using namespace o2::framework::expressions; -using namespace o2::aod::hf_cand; namespace { @@ -273,7 +272,7 @@ struct HfSelOptimisation { for (int iDecay{0}; iDecay < n2Prong; ++iDecay) { if (TESTBIT(cand2Prong.hfflag(), iDecay)) { if (std::abs(cand2Prong.flagMCMatchRec()) == BIT(iDecay)) { - if (cand2Prong.originMCRec() == OriginType::Prompt) { + if (cand2Prong.originMCRec() == RecoDecay::OriginType::Prompt) { isPrompt = true; switch (iDecay) { case o2::aod::hf_cand_prong2::DecayType::D0ToPiK: @@ -283,7 +282,7 @@ struct HfSelOptimisation { testSelections2Prong(cand2Prong, tracks); break; } - } else if (cand2Prong.originMCRec() == OriginType::NonPrompt) { + } else if (cand2Prong.originMCRec() == RecoDecay::OriginType::NonPrompt) { isNonPrompt = true; switch (iDecay) { case o2::aod::hf_cand_prong2::DecayType::D0ToPiK: @@ -329,7 +328,7 @@ struct HfSelOptimisation { for (int iDecay{0}; iDecay < n3Prong; ++iDecay) { if (TESTBIT(cand3Prong.hfflag(), iDecay)) { if (std::abs(cand3Prong.flagMCMatchRec()) == BIT(iDecay)) { - if (cand3Prong.originMCRec() == OriginType::Prompt) { + if (cand3Prong.originMCRec() == RecoDecay::OriginType::Prompt) { isPrompt = true; switch (iDecay) { case o2::aod::hf_cand_prong3::DecayType::DPlusToPiKPi: @@ -345,7 +344,7 @@ struct HfSelOptimisation { testSelections3Prong(cand3Prong, tracks); break; } - } else if (cand3Prong.originMCRec() == OriginType::NonPrompt) { + } else if (cand3Prong.originMCRec() == RecoDecay::OriginType::NonPrompt) { isNonPrompt = true; switch (iDecay) { case o2::aod::hf_cand_prong3::DecayType::DPlusToPiKPi: diff --git a/PWGHF/Tasks/taskD0.cxx b/PWGHF/Tasks/taskD0.cxx index 20aa1376e64..efa3a22f29b 100644 --- a/PWGHF/Tasks/taskD0.cxx +++ b/PWGHF/Tasks/taskD0.cxx @@ -23,7 +23,6 @@ using namespace o2; using namespace o2::framework; using namespace o2::framework::expressions; -using namespace o2::aod::hf_cand; using namespace o2::aod::hf_cand_prong2; using namespace o2::analysis::hf_cuts_d0_topik; @@ -201,7 +200,7 @@ struct TaskD0 { registry.fill(HIST("hPtvsYRecSig_RecoPID"), ptRec, yRec); } - if (candidate.originMCRec() == OriginType::Prompt) { + if (candidate.originMCRec() == RecoDecay::OriginType::Prompt) { registry.fill(HIST("hPtRecSigPrompt"), ptRec); // rec. level pT, prompt if (candidate.isRecoHFFlag() >= d_selectionHFFlag) { registry.fill(HIST("hPtvsYRecSigPrompt_RecoHFFlag"), ptRec, yRec); @@ -314,7 +313,7 @@ struct TaskD0 { auto yGen = RecoDecay::y(array{particle.px(), particle.py(), particle.pz()}, RecoDecay::getMassPDG(particle.pdgCode())); registry.fill(HIST("hPtGen"), ptGen); registry.fill(HIST("hPtvsYGen"), ptGen, yGen); - if (particle.originMCGen() == OriginType::Prompt) { + if (particle.originMCGen() == RecoDecay::OriginType::Prompt) { registry.fill(HIST("hPtGenPrompt"), ptGen); registry.fill(HIST("hPtvsYGenPrompt"), ptGen, yGen); } else { diff --git a/PWGHF/Tasks/taskDPlus.cxx b/PWGHF/Tasks/taskDPlus.cxx index b0ec4bfddd9..e20f900dd1d 100644 --- a/PWGHF/Tasks/taskDPlus.cxx +++ b/PWGHF/Tasks/taskDPlus.cxx @@ -24,7 +24,6 @@ using namespace o2; using namespace o2::framework; using namespace o2::framework::expressions; -using namespace o2::aod::hf_cand; using namespace o2::aod::hf_cand_prong3; #include "Framework/runDataProcessing.h" @@ -159,7 +158,7 @@ struct TaskDPlus { if (candidate.isSelDplusToPiKPi() >= d_selectionFlagDPlus) { registry.fill(HIST("hPtRecSig"), ptRec); // rec. level pT } - if (candidate.originMCRec() == OriginType::Prompt) { + if (candidate.originMCRec() == RecoDecay::OriginType::Prompt) { registry.fill(HIST("hPtvsYRecSigPrompt_RecoSkim"), ptRec, yRec); if (TESTBIT(candidate.isSelDplusToPiKPi(), aod::SelectionStep::RecoTopol)) { registry.fill(HIST("hPtvsYRecSigPrompt_RecoTopol"), ptRec, yRec); @@ -201,7 +200,7 @@ struct TaskDPlus { } registry.fill(HIST("hPtGen"), ptGen); registry.fill(HIST("hPtvsYGen"), ptGen, yGen); - if (particle.originMCGen() == OriginType::Prompt) { + if (particle.originMCGen() == RecoDecay::OriginType::Prompt) { registry.fill(HIST("hPtGenPrompt"), ptGen); registry.fill(HIST("hPtvsYGenPrompt"), ptGen, yGen); } else { diff --git a/PWGHF/Tasks/taskLc.cxx b/PWGHF/Tasks/taskLc.cxx index 549384cc195..835ad4207e5 100644 --- a/PWGHF/Tasks/taskLc.cxx +++ b/PWGHF/Tasks/taskLc.cxx @@ -26,7 +26,6 @@ using namespace o2; using namespace o2::framework; using namespace o2::framework::expressions; -using namespace o2::aod::hf_cand; using namespace o2::aod::hf_cand_prong3; using namespace o2::analysis::hf_cuts_lc_topkpi; @@ -151,7 +150,7 @@ struct TaskLc { registry.fill(HIST("hPtGenSig"), particleMother.pt()); // gen. level pT auto ptRec = candidate.pt(); registry.fill(HIST("hPtRecSig"), ptRec); // rec. level pT - if (candidate.originMCRec() == OriginType::Prompt) { + if (candidate.originMCRec() == RecoDecay::OriginType::Prompt) { registry.fill(HIST("hPtRecSigPrompt"), ptRec); // rec. level pT, prompt } else { registry.fill(HIST("hPtRecSigNonPrompt"), ptRec); // rec. level pT, non-prompt @@ -173,7 +172,7 @@ struct TaskLc { } auto ptGen = particle.pt(); registry.fill(HIST("hPtGen"), ptGen); - if (particle.originMCGen() == OriginType::Prompt) { + if (particle.originMCGen() == RecoDecay::OriginType::Prompt) { registry.fill(HIST("hPtGenPrompt"), ptGen); } else { registry.fill(HIST("hPtGenNonPrompt"), ptGen); diff --git a/PWGHF/Tasks/taskLcCentrality.cxx b/PWGHF/Tasks/taskLcCentrality.cxx index 98c19f3ff61..7bc4d2040b0 100644 --- a/PWGHF/Tasks/taskLcCentrality.cxx +++ b/PWGHF/Tasks/taskLcCentrality.cxx @@ -25,7 +25,6 @@ using namespace o2; using namespace o2::framework; using namespace o2::framework::expressions; -using namespace o2::aod::hf_cand; using namespace o2::aod::hf_cand_prong3; using namespace o2::analysis::hf_cuts_lc_topkpi; @@ -160,7 +159,7 @@ struct TaskLcCentralityMC { registry.fill(HIST("hPtGenSig"), particleMother.pt()); // gen. level pT auto ptRec = candidate.pt(); registry.fill(HIST("hPtRecSig"), ptRec); // rec. level pT - if (candidate.originMCRec() == OriginType::Prompt) { + if (candidate.originMCRec() == RecoDecay::OriginType::Prompt) { registry.fill(HIST("hPtRecSigPrompt"), ptRec); // rec. level pT, prompt } else { registry.fill(HIST("hPtRecSigNonPrompt"), ptRec); // rec. level pT, non-prompt @@ -182,7 +181,7 @@ struct TaskLcCentralityMC { } auto ptGen = particle.pt(); registry.fill(HIST("hPtGen"), ptGen); - if (particle.originMCGen() == OriginType::Prompt) { + if (particle.originMCGen() == RecoDecay::OriginType::Prompt) { registry.fill(HIST("hPtGenPrompt"), ptGen); } else { registry.fill(HIST("hPtGenNonPrompt"), ptGen);