Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
70 changes: 70 additions & 0 deletions Common/Core/RecoDecay.h
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down Expand Up @@ -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 <typename T>
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<std::vector<long int>> arrayIds{};
std::vector<long int> 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<long int> 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<std::tuple<int, double>> mListMass; ///< list of particle masses in form (PDG code, mass)
};
Expand Down
11 changes: 5 additions & 6 deletions EventFiltering/PWGHF/HFFilter.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -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;

Expand Down Expand Up @@ -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 {
Expand Down Expand Up @@ -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;
Expand Down
4 changes: 0 additions & 4 deletions PWGHF/DataModel/HFSecondaryVertex.h
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
4 changes: 2 additions & 2 deletions PWGHF/TableProducer/HFCandidateCreator2Prong.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down Expand Up @@ -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);
Expand Down
4 changes: 2 additions & 2 deletions PWGHF/TableProducer/HFCandidateCreator3Prong.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down Expand Up @@ -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);
Expand Down
2 changes: 1 addition & 1 deletion PWGHF/TableProducer/HFCandidateCreatorChic.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -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);
}
Expand Down
4 changes: 2 additions & 2 deletions PWGHF/TableProducer/HFCandidateCreatorX.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down Expand Up @@ -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);
Expand Down
15 changes: 7 additions & 8 deletions PWGHF/Tasks/HFMCValidation.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -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;

Expand Down Expand Up @@ -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]++;
}
}
Expand Down Expand Up @@ -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());
}
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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;
Expand Down
9 changes: 4 additions & 5 deletions PWGHF/Tasks/HFSelOptimisation.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,6 @@
using namespace o2;
using namespace o2::framework;
using namespace o2::framework::expressions;
using namespace o2::aod::hf_cand;

namespace
{
Expand Down Expand Up @@ -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:
Expand All @@ -283,7 +282,7 @@ struct HfSelOptimisation {
testSelections2Prong<o2::aod::hf_cand_prong2::DecayType::JpsiToEE, 0>(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:
Expand Down Expand Up @@ -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:
Expand All @@ -345,7 +344,7 @@ struct HfSelOptimisation {
testSelections3Prong<o2::aod::hf_cand_prong3::DecayType::XicToPKPi, 0>(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:
Expand Down
5 changes: 2 additions & 3 deletions PWGHF/Tasks/taskD0.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -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;

Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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 {
Expand Down
5 changes: 2 additions & 3 deletions PWGHF/Tasks/taskDPlus.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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 {
Expand Down
5 changes: 2 additions & 3 deletions PWGHF/Tasks/taskLc.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -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;

Expand Down Expand Up @@ -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
Expand All @@ -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);
Expand Down
5 changes: 2 additions & 3 deletions PWGHF/Tasks/taskLcCentrality.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -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;

Expand Down Expand Up @@ -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
Expand All @@ -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);
Expand Down