Skip to content

Commit

Permalink
PWGHF, RecoDecay: Fix search of mothers to work also at partonic level (
Browse files Browse the repository at this point in the history 
AliceO2Group#939)
  • Loading branch information
@fgrosa
fgrosa authored Jul 6, 2022
1 parent 1acfd89 commit eeab638
Show file tree
Hide file tree
Showing 19 changed files with 75 additions and 62 deletions.
85 changes: 49 additions & 36 deletions Common/Core/RecoDecay.h
View file
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
View file
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
View file
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
View file
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
View file
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
View file
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
View file
Original file line number Diff line number Diff line change
Expand Up @@ -157,7 +157,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 @@ -293,7 +293,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 @@ -358,7 +358,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
View file
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
View file
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
View file
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
View file
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
View file
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
View file
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
View file
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
View file
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
View file
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

0 comments on commit eeab638

Please sign in to comment.