Skip to content

Commit 5c16e33

Browse files
mhemmer-cernmhemmer-cern
authored
[PWGEM,PWGEM-36] PM: Fix bug in MCUtilitits isMotherPDG and change … (AliceO2Group#15876)
1 parent c5436b0 commit 5c16e33

4 files changed

Lines changed: 100 additions & 36 deletions

File tree

‎PWGEM/PhotonMeson/Tasks/Pi0EtaToGammaGammaEMCEMC.cxx‎

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -14,6 +14,7 @@
1414
/// \author D. Sekihata, daiki.sekihata@cern.ch
1515

1616
#include "PWGEM/PhotonMeson/Core/Pi0EtaToGammaGamma.h"
17+
#include "PWGEM/PhotonMeson/DataModel/GammaTablesRedux.h"
1718
#include "PWGEM/PhotonMeson/DataModel/gammaTables.h"
1819
#include "PWGEM/PhotonMeson/Utils/PairUtilities.h"
1920

@@ -26,11 +27,11 @@ using namespace o2::aod;
2627
using namespace o2::framework;
2728
using namespace o2::aod::pwgem::photonmeson::photonpair;
2829

29-
using MyEMCClusters = soa::Join<aod::EmEmcClusters, aod::EMCEMEventIds>;
30+
using MyEMCClusters = soa::Join<aod::EMCEMEventIds, aod::MinClusters>;
3031

3132
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
3233
{
3334
return WorkflowSpec{
34-
adaptAnalysisTask<Pi0EtaToGammaGamma<PairType::kEMCEMC, MyEMCClusters, EmEmcMTracks, EmEmcMSTracks>>(cfgc, TaskName{"pi0eta-to-gammagamma-emcemc"}),
35+
adaptAnalysisTask<Pi0EtaToGammaGamma<PairType::kEMCEMC, MyEMCClusters, MinMTracks, MinMSTracks>>(cfgc, TaskName{"pi0eta-to-gammagamma-emcemc"}),
3536
};
3637
}

‎PWGEM/PhotonMeson/Tasks/Pi0EtaToGammaGammaMCEMCEMC.cxx‎

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -14,6 +14,7 @@
1414
/// \author D. Sekihata, daiki.sekihata@cern.ch
1515

1616
#include "PWGEM/PhotonMeson/Core/Pi0EtaToGammaGammaMC.h"
17+
#include "PWGEM/PhotonMeson/DataModel/GammaTablesRedux.h"
1718
#include "PWGEM/PhotonMeson/DataModel/gammaTables.h"
1819
#include "PWGEM/PhotonMeson/Utils/PairUtilities.h"
1920

@@ -26,11 +27,11 @@ using namespace o2::aod;
2627
using namespace o2::framework;
2728
using namespace o2::aod::pwgem::photonmeson::photonpair;
2829

29-
using MyEMCClusters = soa::Join<aod::EmEmcClusters, aod::EMEMCClusterMCLabels, aod::EMCEMEventIds>;
30+
using MyEMCClusters = soa::Join<aod::MinClusters, aod::EMEMCClusterMCLabels, aod::EMCEMEventIds>;
3031

3132
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
3233
{
3334
return WorkflowSpec{
34-
adaptAnalysisTask<Pi0EtaToGammaGammaMC<PairType::kEMCEMC, MyEMCClusters, EmEmcMTracks, EmEmcMSTracks>>(cfgc, TaskName{"pi0eta-to-gammagamma-mc-emcemc"}),
35+
adaptAnalysisTask<Pi0EtaToGammaGammaMC<PairType::kEMCEMC, MyEMCClusters, MinMTracks, MinMSTracks>>(cfgc, TaskName{"pi0eta-to-gammagamma-mc-emcemc"}),
3536
};
3637
}

‎PWGEM/PhotonMeson/Tasks/photonResoTask.cxx‎

Lines changed: 63 additions & 26 deletions
Original file line numberDiff line numberDiff line change
@@ -86,8 +86,10 @@ struct PhotonResoTask {
8686
// configurable axis
8787
ConfigurableAxis thnConfigAxisInvMass{"thnConfigAxisInvMass", {400, 0.0, 0.8}, "invariant mass axis for the neutral meson"};
8888
ConfigurableAxis thnConfigAxisPt{"thnConfigAxisPt", {400, 0., 20.}, "pT axis for the neutral meson"};
89-
ConfigurableAxis thnConfigAxisERelative{"thnConfigAxisERelative", {400, -1., 19.}, "(E rec - E true) / E true axis"};
90-
ConfigurableAxis thnConfigAxisPRelative{"thnConfigAxisPRelative", {60, -1., 2.}, "(P rec - P true) / P true axis"};
89+
ConfigurableAxis thnConfigAxisERelative{"thnConfigAxisERelative", {600, -1., 5.}, "(E rec - E true) / E true axis"};
90+
ConfigurableAxis thnConfigAxisPRelative{"thnConfigAxisPRelative", {300, -1., 2.}, "(P rec - P true) / P true axis"};
91+
ConfigurableAxis thnConfigAxisEtaRelative{"thnConfigAxisEtaRelative", {300, -1., 2.}, "(eta rec - eta true) / eta true axis"};
92+
ConfigurableAxis thnConfigAxisPhiRelative{"thnConfigAxisPhiRelative", {300, -1., 2.}, "(phi rec - phi true) / phi true axis"};
9193
ConfigurableAxis thnConfigAxisCent{"thnConfigAxisCent", {20, 0., 100.}, "centrality axis for the current event"};
9294
ConfigurableAxis thnConfigAxisMult{"thnConfigAxisMult", {60, 0., 60000.}, "multiplicity axis for the current event"};
9395
Configurable<bool> useCent{"useCent", 0, "flag to enable usage of centrality instead of multiplicity as axis."};
@@ -296,11 +298,11 @@ struct PhotonResoTask {
296298
const AxisSpec thnAxisERelative{thnConfigAxisERelative, "#it{E}_{Rec} - #it{E}_{Gen} / #it{E}_{Gen}"};
297299
const AxisSpec thnAxisInvMass{thnConfigAxisInvMass, "#it{M}_{#gamma#gamma} (GeV/#it{c}^{2})"};
298300

299-
const AxisSpec thnAxisEtaGen{280, -0.7, 0.7, "#it{#eta}_{Gen}"};
300-
const AxisSpec thnAxisEtaRec{280, -0.7, 0.7, "#it{#eta}_{Rec}"};
301+
const AxisSpec thnAxisEtaRelative{thnConfigAxisEtaRelative, "#it{#eta}_{Rec} - #it{#eta}_{Gen} / #it{#eta}_{Gen}"};
302+
const AxisSpec thnAxisPhiRelative{thnConfigAxisPhiRelative, "#it{#varphi}_{Rec} - #it{#varphi}_{Gen} / #it{#varphi}_{Gen}"};
301303

304+
const AxisSpec thnAxisEtaGen{280, -0.7, 0.7, "#it{#eta}_{Gen}"};
302305
const AxisSpec thnAxisPhiGen{360, 0., o2::constants::math::TwoPI, "#it{#varphi}_{Gen} (rad)"};
303-
const AxisSpec thnAxisPhiRec{360, 0., o2::constants::math::TwoPI, "#it{#varphi}_{Rec} (rad)"};
304306

305307
AxisSpec thnAxisCentOrMult{1, 0., 1., "Centrality/Multiplicity"}; // placeholder, overwritten in init
306308
if (useCent.value) {
@@ -313,24 +315,28 @@ struct PhotonResoTask {
313315

314316
registry.add("EMCal/hPhotonReso", "EMCal photon rec pT vs true pT vs cent", HistType::kTH3D, {thnAxisPtRec, thnAxisPtGen, thnAxisCentOrMult});
315317
registry.add("EMCal/hConvPhotonReso", "EMCal conversion photon rec pT vs true pT vs cent ", HistType::kTH3D, {thnAxisPtRec, thnAxisPtGen, thnAxisCentOrMult});
318+
registry.add("EMCal/hFullConvPhotonReso", "full EMCal conversion photon rec pT vs true pT vs cent ", HistType::kTH3D, {thnAxisPtRec, thnAxisPtGen, thnAxisCentOrMult});
316319

317320
registry.add("EMCal/hErecEmcPhotons", "EMCal photon rec E - true E vs true E vs cent", HistType::kTH3D, {thnAxisERelative, thnAxisEGen, thnAxisCentOrMult});
318321
registry.add("EMCal/hErecEmcConvPhotons", "EMCal conversion photon rec E - true E vs true E vs cent ", HistType::kTH3D, {thnAxisERelative, thnAxisEGen, thnAxisCentOrMult});
322+
registry.add("EMCal/hErecEmcFullConvPhotons", "full EMCal conversion photon rec E - true E vs true E vs cent ", HistType::kTH3D, {thnAxisERelative, thnAxisEGen, thnAxisCentOrMult});
319323

320324
registry.add("EMCal/hPi0Reso", "EMCal pi0 rec pT vs true pT vs min vs cent ", HistType::kTHnSparseF, {thnAxisPtRec, thnAxisPtGen, thnConfigAxisInvMass, thnAxisCentOrMult});
321325
registry.add("EMCal/hEtaReso", "EMCal eta rec pT vs true pT vs min vs cent ", HistType::kTHnSparseF, {thnAxisPtRec, thnAxisPtGen, thnConfigAxisInvMass, thnAxisCentOrMult});
322326

323-
registry.add("EMCal/hPhotonResoEta", "EMCal photon rec eta vs true eta vs cent", HistType::kTH3D, {thnAxisEtaRec, thnAxisEtaGen, thnAxisCentOrMult});
324-
registry.add("EMCal/hConvPhotonResoEta", "EMCal conversion photon rec eta vs true eta vs cent ", HistType::kTH3D, {thnAxisEtaRec, thnAxisEtaGen, thnAxisCentOrMult});
327+
registry.add("EMCal/hPhotonResoEta", "EMCal photon rec eta vs true eta vs cent", HistType::kTH3D, {thnAxisEtaGen, thnAxisEtaRelative, thnAxisCentOrMult});
328+
registry.add("EMCal/hConvPhotonResoEta", "EMCal conversion photon rec eta vs true eta vs cent ", HistType::kTH3D, {thnAxisEtaGen, thnAxisEtaRelative, thnAxisCentOrMult});
329+
registry.add("EMCal/hFullConvPhotonResoEta", "full EMCal conversion photon rec eta vs true eta vs cent ", HistType::kTH3D, {thnAxisEtaGen, thnAxisEtaRelative, thnAxisCentOrMult});
325330

326-
registry.add("EMCal/hPi0ResoEta", "EMCal pi0 rec eta vs true eta vs min vs cent ", HistType::kTHnSparseF, {thnAxisEtaRec, thnAxisEtaGen, thnConfigAxisInvMass, thnAxisCentOrMult});
327-
registry.add("EMCal/hEtaResoEta", "EMCal eta rec eta vs true eta vs min vs cent ", HistType::kTHnSparseF, {thnAxisEtaRec, thnAxisEtaGen, thnConfigAxisInvMass, thnAxisCentOrMult});
331+
registry.add("EMCal/hPi0ResoEta", "EMCal pi0 rec eta vs true eta vs min vs cent ", HistType::kTHnSparseF, {thnAxisEtaGen, thnAxisEtaRelative, thnConfigAxisInvMass, thnAxisCentOrMult});
332+
registry.add("EMCal/hEtaResoEta", "EMCal eta rec eta vs true eta vs min vs cent ", HistType::kTHnSparseF, {thnAxisEtaGen, thnAxisEtaRelative, thnConfigAxisInvMass, thnAxisCentOrMult});
328333

329-
registry.add("EMCal/hPhotonResoPhi", "EMCal photon rec phi vs true phi vs cent", HistType::kTH3D, {thnAxisPhiRec, thnAxisPhiGen, thnAxisCentOrMult});
330-
registry.add("EMCal/hConvPhotonResoPhi", "EMCal conversion photon rec phi vs true phi vs cent ", HistType::kTH3D, {thnAxisPhiRec, thnAxisPhiGen, thnAxisCentOrMult});
334+
registry.add("EMCal/hPhotonResoPhi", "EMCal photon rec phi vs true phi vs cent", HistType::kTH3D, {thnAxisPhiGen, thnAxisPhiRelative, thnAxisCentOrMult});
335+
registry.add("EMCal/hConvPhotonResoPhi", "EMCal conversion photon rec phi vs true phi vs cent ", HistType::kTH3D, {thnAxisPhiGen, thnAxisPhiRelative, thnAxisCentOrMult});
336+
registry.add("EMCal/hFullConvPhotonResoPhi", "full EMCal conversion photon rec phi vs true phi vs cent ", HistType::kTH3D, {thnAxisPhiGen, thnAxisPhiRelative, thnAxisCentOrMult});
331337

332-
registry.add("EMCal/hPi0ResoPhi", "EMCal pi0 rec phi vs true phi vs min vs cent ", HistType::kTHnSparseF, {thnAxisPhiRec, thnAxisPhiGen, thnConfigAxisInvMass, thnAxisCentOrMult});
333-
registry.add("EMCal/hEtaResoPhi", "EMCal eta rec phi vs true phi vs min vs cent ", HistType::kTHnSparseF, {thnAxisPhiRec, thnAxisPhiGen, thnConfigAxisInvMass, thnAxisCentOrMult});
338+
registry.add("EMCal/hPi0ResoPhi", "EMCal pi0 rec phi vs true phi vs min vs cent ", HistType::kTHnSparseF, {thnAxisPhiGen, thnAxisPhiRelative, thnConfigAxisInvMass, thnAxisCentOrMult});
339+
registry.add("EMCal/hEtaResoPhi", "EMCal eta rec phi vs true phi vs min vs cent ", HistType::kTHnSparseF, {thnAxisPhiGen, thnAxisPhiRelative, thnConfigAxisInvMass, thnAxisCentOrMult});
334340

335341
registry.add("PCM/hPhotonReso", "PCM photon rec pT vs true pT vs cent", HistType::kTH3D, {thnAxisPtRec, thnAxisPtGen, thnAxisCentOrMult});
336342

@@ -460,6 +466,7 @@ struct PhotonResoTask {
460466
// create iterators for photon mc particles
461467
auto mcPhoton1 = mcParticles.begin();
462468
auto mcPhoton2 = mcParticles.begin();
469+
auto mcConvLeg = mcParticles.begin();
463470

464471
// leg iterators for PCM
465472
auto pos1 = legs.begin();
@@ -489,21 +496,49 @@ struct PhotonResoTask {
489496
// we only want to look at the largest contribution
490497
mcPhoton1.setCursor(photonEMC.emmcparticleIds()[0]);
491498

499+
// if the largest contribution is photon, just fill the photon histograms
492500
if (std::abs(mcPhoton1.pdgCode()) == PDG_t::kGamma) {
493501
registry.fill(HIST("EMCal/hPhotonReso"), photonEMC.pt(), mcPhoton1.pt(), centOrMult);
494-
registry.fill(HIST("EMCal/hPhotonResoEta"), photonEMC.eta(), mcPhoton1.eta(), centOrMult);
495-
registry.fill(HIST("EMCal/hPhotonResoPhi"), photonEMC.phi(), mcPhoton1.phi(), centOrMult);
502+
registry.fill(HIST("EMCal/hPhotonResoEta"), mcPhoton1.eta(), (photonEMC.eta() - mcPhoton1.eta()) / mcPhoton1.eta(), centOrMult);
503+
registry.fill(HIST("EMCal/hPhotonResoPhi"), mcPhoton1.phi(), (photonEMC.phi() - mcPhoton1.phi()) / mcPhoton1.phi(), centOrMult);
496504
registry.fill(HIST("EMCal/hErecEmcPhotons"), (photonEMC.e() - mcPhoton1.e()) / mcPhoton1.e(), mcPhoton1.e(), centOrMult);
497505
} else if (std::abs(mcPhoton1.pdgCode()) == PDG_t::kElectron) {
498-
if (!o2::aod::pwgem::photonmeson::utils::mcutil::isMotherPDG(mcPhoton1, PDG_t::kGamma)) {
506+
// if largest contribution is e+ or e-, check if its from a photon conversion
507+
int32_t mcPhoton1MotherId = o2::aod::pwgem::photonmeson::utils::mcutil::getMotherIndexFromChain(mcPhoton1, PDG_t::kGamma);
508+
int32_t mcConvLegMotherId = -1;
509+
bool hasBothLegs = false;
510+
if (mcPhoton1MotherId == -1) {
499511
continue;
500512
}
501-
registry.fill(HIST("EMCal/hConvPhotonReso"), photonEMC.pt(), mcPhoton1.pt(), centOrMult);
502-
registry.fill(HIST("EMCal/hConvPhotonResoEta"), photonEMC.eta(), mcPhoton1.eta(), centOrMult);
503-
registry.fill(HIST("EMCal/hConvPhotonResoPhi"), photonEMC.phi(), mcPhoton1.phi(), centOrMult);
504-
registry.fill(HIST("EMCal/hErecEmcConvPhotons"), (photonEMC.e() - mcPhoton1.e()) / mcPhoton1.e(), mcPhoton1.e(), centOrMult);
505-
}
506-
}
513+
// mcPhoton1 now points to a photon that produced the e+/e- that was the largest contributor in the cluster
514+
// check if both legs of the photon hit the cluster to decide which conversion histogram to fill
515+
for (int32_t emmcparticleId = 1; emmcparticleId < static_cast<int32_t>(photonEMC.emmcparticleIds().size()); ++emmcparticleId) {
516+
mcConvLeg.setCursor(photonEMC.emmcparticleIds()[emmcparticleId]);
517+
const int convLegPDG = mcConvLeg.pdgCode(); // store the pdg value, because getMotherIndexFromChain will alter the cursor!
518+
mcConvLegMotherId = o2::aod::pwgem::photonmeson::utils::mcutil::getMotherIndexFromChain(mcConvLeg, PDG_t::kGamma);
519+
if (!(std::abs(convLegPDG) == PDG_t::kElectron) || mcConvLegMotherId == -1) {
520+
continue;
521+
}
522+
if (mcPhoton1MotherId == mcConvLegMotherId) {
523+
hasBothLegs = true;
524+
break;
525+
}
526+
}
527+
if (hasBothLegs) {
528+
// if its a conversion with both legs, fill the full conversion histograms
529+
registry.fill(HIST("EMCal/hFullConvPhotonReso"), photonEMC.pt(), mcPhoton1.pt(), centOrMult);
530+
registry.fill(HIST("EMCal/hFullConvPhotonResoEta"), mcPhoton1.eta(), (photonEMC.eta() - mcPhoton1.eta()) / mcPhoton1.eta(), centOrMult);
531+
registry.fill(HIST("EMCal/hFullConvPhotonResoPhi"), mcPhoton1.phi(), (photonEMC.phi() - mcPhoton1.phi()) / mcPhoton1.phi(), centOrMult);
532+
registry.fill(HIST("EMCal/hErecEmcFullConvPhotons"), (photonEMC.e() - mcPhoton1.e()) / mcPhoton1.e(), mcPhoton1.e(), centOrMult);
533+
} else {
534+
// if its a conversion with only one leg, fill the conversion histograms
535+
registry.fill(HIST("EMCal/hConvPhotonReso"), photonEMC.pt(), mcPhoton1.pt(), centOrMult);
536+
registry.fill(HIST("EMCal/hConvPhotonResoEta"), mcPhoton1.eta(), (photonEMC.eta() - mcPhoton1.eta()) / mcPhoton1.eta(), centOrMult);
537+
registry.fill(HIST("EMCal/hConvPhotonResoPhi"), mcPhoton1.phi(), (photonEMC.phi() - mcPhoton1.phi()) / mcPhoton1.phi(), centOrMult);
538+
registry.fill(HIST("EMCal/hErecEmcConvPhotons"), (photonEMC.e() - mcPhoton1.e()) / mcPhoton1.e(), mcPhoton1.e(), centOrMult);
539+
}
540+
} // else if (std::abs(mcPhoton1.pdgCode()) == PDG_t::kElectron)
541+
} // for (const auto& photonEMC : photonsEMCPerCollision) {
507542

508543
for (const auto& photonPCM : photonsPCMPerCollision) {
509544
if (!(v0flags.test(photonPCM.globalIndex()))) {
@@ -543,6 +578,8 @@ struct PhotonResoTask {
543578
ROOT::Math::PtEtaPhiMVector v2(g2.pt(), g2.eta(), g2.phi(), 0.);
544579
ROOT::Math::PtEtaPhiMVector vMeson = v1 + v2;
545580

581+
const float mesonPhi = RecoDecay::constrainAngle(vMeson.Phi());
582+
546583
float openingAngle = std::acos(v1.Vect().Dot(v2.Vect()) / (v1.P() * v2.P()));
547584

548585
registry.fill(HIST("hMesonCuts"), 1);
@@ -586,14 +623,14 @@ struct PhotonResoTask {
586623
if (pi0id >= 0) {
587624
const auto pi0mc = mcParticles.iteratorAt(pi0id);
588625
registry.fill(HIST("EMCal/hPi0Reso"), vMeson.Pt(), pi0mc.pt(), vMeson.M(), centOrMult);
589-
registry.fill(HIST("EMCal/hPi0ResoEta"), vMeson.Eta(), pi0mc.eta(), vMeson.M(), centOrMult);
590-
registry.fill(HIST("EMCal/hPi0ResoPhi"), RecoDecay::constrainAngle(vMeson.Phi()), pi0mc.phi(), vMeson.M(), centOrMult);
626+
registry.fill(HIST("EMCal/hPi0ResoEta"), pi0mc.eta(), (vMeson.Eta() - pi0mc.eta()) / pi0mc.eta(), vMeson.M(), centOrMult);
627+
registry.fill(HIST("EMCal/hPi0ResoPhi"), pi0mc.phi(), (mesonPhi - pi0mc.phi()) / pi0mc.phi(), vMeson.M(), centOrMult);
591628
}
592629
if (etaid >= 0) {
593630
const auto etamc = mcParticles.iteratorAt(etaid);
594631
registry.fill(HIST("EMCal/hEtaReso"), vMeson.Pt(), etamc.pt(), vMeson.M(), centOrMult);
595-
registry.fill(HIST("EMCal/hEtaResoEta"), vMeson.Eta(), etamc.eta(), vMeson.M(), centOrMult);
596-
registry.fill(HIST("EMCal/hEtaResoPhi"), RecoDecay::constrainAngle(vMeson.Phi()), etamc.phi(), vMeson.M(), centOrMult);
632+
registry.fill(HIST("EMCal/hEtaResoEta"), etamc.eta(), (vMeson.Eta() - etamc.eta()) / etamc.eta(), vMeson.M(), centOrMult);
633+
registry.fill(HIST("EMCal/hEtaResoPhi"), etamc.phi(), (mesonPhi - etamc.phi()) / etamc.phi(), vMeson.M(), centOrMult);
597634
}
598635
}
599636
}

‎PWGEM/PhotonMeson/Utils/MCUtilities.h‎

Lines changed: 31 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -282,21 +282,46 @@ bool isGammaGammaDecay(TMCParticle const& mcParticle, TMCParticles const& mcPart
282282
}
283283

284284
//_______________________________________________________________________
285-
// Go up the decay chain of a mcparticle looking for a mother with the given pdg codes, if found return true else false
286-
// E.g. if electron cluster is coming from a photon return true, if primary electron return false
285+
/// \brief Go up the decay chain of a mcparticle looking for a mother with the given pdg codes, if found return true else false
286+
/// E.g. if electron cluster is coming from a photon return true, if primary electron return false
287+
/// \param mcparticle iterator of mxparticle, WILL BE CHANGED by this function!
288+
/// \param motherPDG target mother PDG value
289+
/// \param depth how many steps in the chain this check should go maximum before failing
287290
template <o2::soa::is_iterator T>
288-
bool isMotherPDG(T& mcparticle, const int motherPDG, const int Depth = 10) // o2-linter: disable=pdg/explicit-code (false positive)
291+
bool isMotherPDG(T& mcparticle, const int motherPDG, const int depth = 10) // o2-linter: disable=pdg/explicit-code (false positive)
289292
{
290-
if (!mcparticle.has_mothers() || Depth < 1) {
293+
if (!mcparticle.has_mothers() || depth < 1) {
291294
return false;
292295
}
293296

294297
int motherid = mcparticle.mothersIds()[0];
295298
mcparticle.setCursor(motherid);
296-
if (mcparticle.pdgCode() != motherPDG) {
299+
if (mcparticle.pdgCode() == motherPDG) {
297300
return true; // The mother has the required pdg code, so return its daughters global mc particle code.
298301
} else {
299-
return isMotherPDG(mcparticle, motherPDG, Depth - 1);
302+
return isMotherPDG(mcparticle, motherPDG, depth - 1);
303+
}
304+
}
305+
306+
//_______________________________________________________________________
307+
/// \brief Go up the decay chain of a mcparticle looking for a mother with the given pdg codes, if found return id else -1
308+
/// E.g. if electron cluster is coming from a photon return true, if primary electron return false
309+
/// \param mcparticle iterator of mxparticle, WILL BE CHANGED by this function!
310+
/// \param motherPDG target mother PDG value
311+
/// \param depth how many steps in the chain this check should go maximum before failing
312+
template <o2::soa::is_iterator T>
313+
int32_t getMotherIndexFromChain(T& mcparticle, const int motherPDG, const int depth = 10) // o2-linter: disable=pdg/explicit-code (false positive)
314+
{
315+
if (!mcparticle.has_mothers() || depth < 1) {
316+
return -1;
317+
}
318+
319+
int32_t motherid = mcparticle.mothersIds()[0];
320+
mcparticle.setCursor(motherid);
321+
if (mcparticle.pdgCode() == motherPDG) {
322+
return motherid; // The mother has the required pdg code, so return its daughters global mc particle code.
323+
} else {
324+
return getMotherIndexFromChain(mcparticle, motherPDG, depth - 1);
300325
}
301326
}
302327

0 commit comments

Comments
 (0)