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
85 changes: 49 additions & 36 deletions Common/Core/RecoDecay.h
Original file line number Diff line number Diff line change
Expand Up @@ -562,62 +562,75 @@ class RecoDecay
}

/// Finds the mother of an MC particle by looking for the expected PDG code in the mother chain.
/// \param particlesMC table with MC particles
/// \param particle MC particle
/// \param PDGMother expected mother PDG code
/// \param acceptAntiParticles switch to accept the antiparticle of the expected mother
/// \param sign antiparticle indicator of the found mother w.r.t. PDGMother; 1 if particle, -1 if antiparticle, 0 if mother not found
/// \param depthMax maximum decay tree level to check; Mothers up to this level will be considered. If -1, all levels are considered.
/// \return index of the mother particle if found, -1 otherwise
template <typename T>
static int getMother(const T& particle,
static int getMother(const T& particlesMC,
const typename T::iterator& particle,
int PDGMother,
bool acceptAntiParticles = false,
int8_t* sign = nullptr,
int8_t depthMax = -1)
{
int8_t sgn = 0; // 1 if the expected mother is particle, -1 if antiparticle (w.r.t. PDGMother)
int indexMother = -1; // index of the final matched mother, if found
auto particleMother = particle; // Initialise loop over mothers.
int stage = 0; // mother tree level (just for debugging)
int8_t sgn = 0; // 1 if the expected mother is particle, -1 if antiparticle (w.r.t. PDGMother)
int indexMother = -1; // index of the final matched mother, if found
int stage = 0; // mother tree level (just for debugging)
bool motherFound = false; // true when the desired mother particle is found in the kine tree
if (sign) {
*sign = sgn;
}
std::vector<int> listId;
listId.push_back(particleMother.globalIndex());
while (particleMother.has_mothers()) {
if (depthMax > -1 && -stage >= depthMax) { // Maximum depth has been reached.
return -1;
}
if (particleMother.mothersIds()[0] == -1) { // Additional check for converted Run 2 MC
return -1;
}
particleMother = particleMother.template mothers_first_as<typename std::decay_t<T>::parent_t>(); // get mother 0
auto indexMotherTmp = particleMother.globalIndex();
// Check mother's PDG code.
auto PDGParticleIMother = particleMother.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, indexMotherTmp);
if (std::find(listId.begin(), listId.end(), indexMotherTmp) != listId.end()) {
LOGF(fatal, "Circular mothership at particle index %d", indexMotherTmp);
return -1;
}
listId.push_back(indexMotherTmp);
if (PDGParticleIMother == PDGMother) { // exact PDG match
sgn = 1;
indexMother = indexMotherTmp;
break;
} else if (acceptAntiParticles && PDGParticleIMother == -PDGMother) { // antiparticle PDG match
sgn = -1;
indexMother = indexMotherTmp;
break;

// 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 (!motherFound && arrayIds[-stage].size() > 0 && (depthMax < 0 || -stage < depthMax)) {
// 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 = 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 (PDGParticleIMother == PDGMother) { // exact PDG match
sgn = 1;
indexMother = iMother;
motherFound = true;
break;
} else if (acceptAntiParticles && PDGParticleIMother == -PDGMother) { // antiparticle PDG match
sgn = -1;
indexMother = iMother;
motherFound = true;
break;
}
// 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 (sign) {
*sign = sgn;
}

return indexMother;
}

Expand Down Expand Up @@ -723,7 +736,7 @@ class RecoDecay
if (iProng == 0) {
// Get the mother index and its sign.
// PDG code of the first daughter's mother determines whether the expected mother is a particle or antiparticle.
indexMother = getMother(particleI, PDGMother, acceptAntiParticles, &sgn, depthMax);
indexMother = getMother(particlesMC, particleI, PDGMother, acceptAntiParticles, &sgn, depthMax);
// Check whether mother was found.
if (indexMother <= -1) {
//Printf("MC Rec: Rejected: bad mother index or PDG");
Expand Down
4 changes: 2 additions & 2 deletions EventFiltering/PWGHF/HFFilter.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -979,7 +979,7 @@ 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(particle, kBottom, true) > -1 ? OriginType::NonPrompt : OriginType::Prompt);
origin = (RecoDecay::getMother(particlesMC, particle, kBottom, true) > -1 ? OriginType::NonPrompt : OriginType::Prompt);
if (origin == OriginType::NonPrompt) {
flag = kNonPrompt;
} else {
Expand Down Expand Up @@ -1038,7 +1038,7 @@ struct HfFilter { // Main struct for HF triggers

if (indexRec > -1) {
auto particle = particlesMC.rawIteratorAt(indexRec);
origin = (RecoDecay::getMother(particle, kBottom, true) > -1 ? OriginType::NonPrompt : OriginType::Prompt);
origin = (RecoDecay::getMother(particlesMC, particle, kBottom, true) > -1 ? OriginType::NonPrompt : OriginType::Prompt);
if (origin == OriginType::NonPrompt) {
flag = kNonPrompt;
} else {
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(particle, kBottom, true) > -1 ? OriginType::NonPrompt : OriginType::Prompt);
origin = (RecoDecay::getMother(particlesMC, particle, kBottom, true) > -1 ? OriginType::NonPrompt : OriginType::Prompt);
}

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(particle, kBottom, true) > -1 ? OriginType::NonPrompt : OriginType::Prompt);
origin = (RecoDecay::getMother(particlesMC, particle, kBottom, true) > -1 ? OriginType::NonPrompt : OriginType::Prompt);
}

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(particle, kBottom, true) > -1 ? OriginType::NonPrompt : OriginType::Prompt);
origin = (RecoDecay::getMother(particlesMC, particle, kBottom, true) > -1 ? OriginType::NonPrompt : OriginType::Prompt);
}

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(particle, kBottom, true) > -1 ? OriginType::NonPrompt : OriginType::Prompt);
origin = (RecoDecay::getMother(particlesMC, particle, kBottom, true) > -1 ? OriginType::NonPrompt : OriginType::Prompt);
}

rowMCMatchGen(flag, origin, channel);
Expand Down
6 changes: 3 additions & 3 deletions PWGHF/TableProducer/HFCandidateCreatorChic.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -231,8 +231,8 @@ struct HFCandidateCreatorChicMC {
if (indexRec > -1) {
hMassJpsiToMuMuMatched->Fill(InvMassJpsiToMuMu(candidate.index0()));

int indexMother = RecoDecay::getMother(particlesMC.rawIteratorAt(indexRec), pdg::Code::kChic1);
int indexMotherGamma = RecoDecay::getMother(particlesMC.rawIteratorAt(candidate.index1().mcparticleId()), pdg::Code::kChic1);
int indexMother = RecoDecay::getMother(particlesMC, particlesMC.rawIteratorAt(indexRec), pdg::Code::kChic1);
int indexMotherGamma = RecoDecay::getMother(particlesMC, particlesMC.rawIteratorAt(candidate.index1().mcparticleId()), pdg::Code::kChic1);
if (indexMother > -1 && indexMotherGamma == indexMother && candidate.index1().mcparticle().pdgCode() == kGamma) {
auto particleMother = particlesMC.rawIteratorAt(indexMother);
hEphotonMatched->Fill(candidate.index1().e());
Expand All @@ -249,7 +249,7 @@ struct HFCandidateCreatorChicMC {
}
if (flag != 0) {
auto particle = particlesMC.rawIteratorAt(indexRec);
origin = (RecoDecay::getMother(particle, kBottom, true) > -1 ? NonPrompt : Prompt);
origin = (RecoDecay::getMother(particlesMC, particle, kBottom, true) > -1 ? NonPrompt : Prompt);
}
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(particle, 5, true) > -1 ? NonPrompt : Prompt);
origin = (RecoDecay::getMother(particlesMC, particle, 5, true) > -1 ? NonPrompt : Prompt);
}

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(particle, 5, true) > -1 ? NonPrompt : Prompt);
origin = (RecoDecay::getMother(particlesMC, particle, 5, true) > -1 ? NonPrompt : Prompt);
}

rowMCMatchGen(flag, origin, channel);
Expand Down
6 changes: 3 additions & 3 deletions PWGHF/Tasks/HFMCValidation.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -153,7 +153,7 @@ 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(particle, kBottom, true) > -1 ? OriginType::NonPrompt : OriginType::Prompt);
origin = (RecoDecay::getMother(particlesMC, particle, kBottom, true) > -1 ? OriginType::NonPrompt : OriginType::Prompt);
if (origin == OriginType::Prompt) {
counterPrompt[iD]++;
} else if (origin == OriginType::NonPrompt) {
Expand Down Expand Up @@ -282,7 +282,7 @@ struct ValidationRecLevel {
if (whichHad >= 0 && whichOrigin >= 0) {
int indexParticle = 0;
if (cand2Prong.index0_as<aod::BigTracksMC>().has_mcParticle()) {
indexParticle = RecoDecay::getMother(cand2Prong.index0_as<aod::BigTracksMC>().mcParticle(), PDGArrayParticle[whichHad], true);
indexParticle = RecoDecay::getMother(particlesMC, cand2Prong.index0_as<aod::BigTracksMC>().mcParticle(), PDGArrayParticle[whichHad], true);
}
auto mother = particlesMC.rawIteratorAt(indexParticle);
histDeltaPt[whichHad]->Fill(cand2Prong.pt() - mother.pt());
Expand Down Expand Up @@ -347,7 +347,7 @@ struct ValidationRecLevel {
if (whichHad >= 0) {
int indexParticle = 0;
if (cand3Prong.index0_as<aod::BigTracksMC>().has_mcParticle()) {
indexParticle = RecoDecay::getMother(cand3Prong.index0_as<aod::BigTracksMC>().mcParticle(), PDGArrayParticle[whichHad], true);
indexParticle = RecoDecay::getMother(particlesMC, cand3Prong.index0_as<aod::BigTracksMC>().mcParticle(), PDGArrayParticle[whichHad], true);
}
auto mother = particlesMC.rawIteratorAt(indexParticle);
histDeltaPt[whichHad]->Fill(cand3Prong.pt() - mother.pt());
Expand Down
2 changes: 1 addition & 1 deletion PWGHF/Tasks/taskBPlus.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -182,7 +182,7 @@ struct HfTaskBplusMc {
}
if (std::abs(candidate.flagMCMatchRec()) == 1 << hf_cand_bplus::DecayType::BPlusToD0Pi) {

auto indexMother = RecoDecay::getMother(candidate.index1_as<aod::BigTracksMC>().mcParticle_as<soa::Join<aod::McParticles, aod::HfCandBPMCGen>>(), pdg::Code::kBPlus, true);
auto indexMother = RecoDecay::getMother(particlesMC, candidate.index1_as<aod::BigTracksMC>().mcParticle_as<soa::Join<aod::McParticles, aod::HfCandBPMCGen>>(), pdg::Code::kBPlus, true);
auto particleMother = particlesMC.rawIteratorAt(indexMother);
registry.fill(HIST("hPtGenSig"), particleMother.pt());
registry.fill(HIST("hPtRecSig"), candidate.pt());
Expand Down
2 changes: 1 addition & 1 deletion PWGHF/Tasks/taskChic.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -165,7 +165,7 @@ struct TaskChicMC {
}
if (candidate.flagMCMatchRec() == 1 << decayMode) {
//FIXME the access to the MC particle gen not yet functional
//int indexMother = RecoDecay::getMother(particlesMC.rawIteratorAt(candidate.index1().mcParticle_as<aod::McParticles_000>().globalIndex()), 20443);
//int indexMother = RecoDecay::getMother(particlesMC, particlesMC.rawIteratorAt(candidate.index1().mcParticle_as<aod::McParticles_000>().globalIndex()), 20443);
//auto particleMother = particlesMC.rawIteratorAt(indexMother);
//registry.fill(HIST("hPtGenSig"), particleMother.pt());
registry.fill(HIST("hPtRecSig"), candidate.pt());
Expand Down
2 changes: 1 addition & 1 deletion PWGHF/Tasks/taskD0.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -182,7 +182,7 @@ struct TaskD0 {
}
if (std::abs(candidate.flagMCMatchRec()) == 1 << DecayType::D0ToPiK) {
// Get the corresponding MC particle.
auto indexMother = RecoDecay::getMother(candidate.index0_as<aod::BigTracksMC>().mcParticle_as<soa::Join<aod::McParticles, aod::HfCandProng2MCGen>>(), pdg::Code::kD0, true);
auto indexMother = RecoDecay::getMother(particlesMC, candidate.index0_as<aod::BigTracksMC>().mcParticle_as<soa::Join<aod::McParticles, aod::HfCandProng2MCGen>>(), pdg::Code::kD0, true);
auto particleMother = particlesMC.rawIteratorAt(indexMother);
registry.fill(HIST("hPtGenSig"), particleMother.pt()); // gen. level pT
auto ptRec = candidate.pt();
Expand Down
2 changes: 1 addition & 1 deletion PWGHF/Tasks/taskDPlus.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -144,7 +144,7 @@ struct TaskDPlus {
}
if (std::abs(candidate.flagMCMatchRec()) == 1 << DecayType::DPlusToPiKPi) {
// Get the corresponding MC particle.
auto indexMother = RecoDecay::getMother(candidate.index0_as<aod::BigTracksMC>().mcParticle_as<soa::Join<aod::McParticles, aod::HfCandProng3MCGen>>(), pdg::Code::kDPlus, true);
auto indexMother = RecoDecay::getMother(particlesMC, candidate.index0_as<aod::BigTracksMC>().mcParticle_as<soa::Join<aod::McParticles, aod::HfCandProng3MCGen>>(), pdg::Code::kDPlus, true);
auto particleMother = particlesMC.rawIteratorAt(indexMother);
registry.fill(HIST("hPtGenSig"), particleMother.pt()); // gen. level pT
auto ptRec = candidate.pt();
Expand Down
2 changes: 1 addition & 1 deletion PWGHF/Tasks/taskJpsi.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -231,7 +231,7 @@ struct TaskJpsiMC {
}
if (candidate.flagMCMatchRec() == 1 << decayMode) {
//Get the corresponding MC particle.
auto indexMother = RecoDecay::getMother(candidate.index0_as<aod::BigTracksMC>().mcParticle_as<McParticlesHf>(), pdg::Code::kJpsi, true);
auto indexMother = RecoDecay::getMother(particlesMC, candidate.index0_as<aod::BigTracksMC>().mcParticle_as<McParticlesHf>(), pdg::Code::kJpsi, true);
auto particleMother = particlesMC.rawIteratorAt(indexMother);
registry.fill(HIST("hPtGenSig"), particleMother.pt()); // gen. level pT
registry.fill(HIST("hPtRecSig"), candidate.pt()); // rec. level pT
Expand Down
2 changes: 1 addition & 1 deletion PWGHF/Tasks/taskLb.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -188,7 +188,7 @@ struct HfTaskLbMc {
auto candLc = candidate.index0_as<aod::HfCandProng3>();
if (std::abs(candidate.flagMCMatchRec()) == 1 << hf_cand_lb::DecayType::LbToLcPi) {

auto indexMother = RecoDecay::getMother(candidate.index1_as<aod::BigTracksMC>().mcParticle_as<soa::Join<aod::McParticles, aod::HfCandLbMCGen>>(), pdg::Code::kLambdaB0, true);
auto indexMother = RecoDecay::getMother(particlesMC, candidate.index1_as<aod::BigTracksMC>().mcParticle_as<soa::Join<aod::McParticles, aod::HfCandLbMCGen>>(), pdg::Code::kLambdaB0, true);
auto particleMother = particlesMC.rawIteratorAt(indexMother);
registry.fill(HIST("hPtGenSig"), particleMother.pt());
registry.fill(HIST("hPtRecSig"), candidate.pt());
Expand Down
2 changes: 1 addition & 1 deletion PWGHF/Tasks/taskLc.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -146,7 +146,7 @@ struct TaskLc {
}
if (std::abs(candidate.flagMCMatchRec()) == 1 << DecayType::LcToPKPi) {
// Get the corresponding MC particle.
auto indexMother = RecoDecay::getMother(candidate.index0_as<aod::BigTracksMC>().mcParticle_as<soa::Join<aod::McParticles, aod::HfCandProng3MCGen>>(), pdg::Code::kLambdaCPlus, true);
auto indexMother = RecoDecay::getMother(particlesMC, candidate.index0_as<aod::BigTracksMC>().mcParticle_as<soa::Join<aod::McParticles, aod::HfCandProng3MCGen>>(), pdg::Code::kLambdaCPlus, true);
auto particleMother = particlesMC.rawIteratorAt(indexMother);
registry.fill(HIST("hPtGenSig"), particleMother.pt()); // gen. level pT
auto ptRec = candidate.pt();
Expand Down
2 changes: 1 addition & 1 deletion PWGHF/Tasks/taskLcCentrality.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -155,7 +155,7 @@ struct TaskLcCentralityMC {
}
if (std::abs(candidate.flagMCMatchRec()) == 1 << DecayType::LcToPKPi) {
// Get the corresponding MC particle.
auto indexMother = RecoDecay::getMother(candidate.index0_as<aod::BigTracksMC>().mcParticle_as<soa::Join<aod::McParticles, aod::HfCandProng3MCGen>>(), pdg::Code::kLambdaCPlus, true);
auto indexMother = RecoDecay::getMother(particlesMC, candidate.index0_as<aod::BigTracksMC>().mcParticle_as<soa::Join<aod::McParticles, aod::HfCandProng3MCGen>>(), pdg::Code::kLambdaCPlus, true);
auto particleMother = particlesMC.rawIteratorAt(indexMother);
registry.fill(HIST("hPtGenSig"), particleMother.pt()); // gen. level pT
auto ptRec = candidate.pt();
Expand Down
2 changes: 1 addition & 1 deletion PWGHF/Tasks/taskLcK0sP.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -115,7 +115,7 @@ struct TaskLcK0SpMC {
}
if (std::abs(candidate.flagMCMatchRec()) == 1) {
// Get the corresponding MC particle.
auto indexMother = RecoDecay::getMother(candidate.index0_as<aod::BigTracksMC>().mcParticle_as<soa::Join<aod::McParticles, aod::HfCandCascadeMCGen>>(), pdg::Code::kLambdaCPlus, true);
auto indexMother = RecoDecay::getMother(particlesMC, candidate.index0_as<aod::BigTracksMC>().mcParticle_as<soa::Join<aod::McParticles, aod::HfCandCascadeMCGen>>(), pdg::Code::kLambdaCPlus, true);
auto particleMother = particlesMC.rawIteratorAt(indexMother);
registry.fill(HIST("hPtGenSig"), particleMother.pt()); // gen. level pT
registry.fill(HIST("hPtRecSig"), candidate.pt()); // rec. level pT
Expand Down
Loading