RecoDecay: getDaughters: Skip PDG check for the original particle. (#5303) (#5263)

Do not label the original particle as final, when its species matches the expected species of its daughter.
This could happen when calling this function for a quark that radiates a gluon and becomes its own daughter.
In general, there is no reason to check the PDG code of the original particle against the codes of the expected daughters, so skipping it speeds up this function.

Co-authored-by: Vít Kučera <vit.kucera@cern.ch>

Author: DelloStritto

Analysis/DataModel/include/AnalysisDataModel/HFSecondaryVertex.h

@@ -314,6 +314,8 @@ DECLARE_SOA_DYNAMIC_COLUMN(MaxNormalisedDeltaIP, maxNormalisedDeltaIP, [](float
 // - ±LcToPKPi: Λc± → p± K∓ π±
 DECLARE_SOA_COLUMN(FlagMCMatchRec, flagMCMatchRec, int8_t); // reconstruction level
 DECLARE_SOA_COLUMN(FlagMCMatchGen, flagMCMatchGen, int8_t); // generator level
+DECLARE_SOA_COLUMN(FlagMCDecayChanRec, flagMCDecayChanRec, int8_t); // Resonant decay channel flag, reconstruction level
+DECLARE_SOA_COLUMN(FlagMCDecayChanGen, flagMCDecayChanGen, int8_t); // Resonant decay channel flag, generator level
 
 // mapping of decay types
 enum DecayType { DPlusToPiKPi = 0,
@@ -428,11 +430,13 @@ using HfCandProng3 = HfCandProng3Ext;
 
 // table with results of reconstruction level MC matching
 DECLARE_SOA_TABLE(HfCandProng3MCRec, "AOD", "HFCANDP3MCREC",
-                  hf_cand_prong3::FlagMCMatchRec);
+                  hf_cand_prong3::FlagMCMatchRec,
+                  hf_cand_prong3::FlagMCDecayChanRec);
 
 // table with results of generator level MC matching
 DECLARE_SOA_TABLE(HfCandProng3MCGen, "AOD", "HFCANDP3MCGEN",
-                  hf_cand_prong3::FlagMCMatchGen);
+                  hf_cand_prong3::FlagMCMatchGen,
+                  hf_cand_prong3::FlagMCDecayChanGen);
 
 } // namespace o2::aod
 

Analysis/Tasks/PWGHF/HFCandidateCreator3Prong.cxx

@@ -152,11 +152,19 @@ struct HFCandidateCreator3ProngMC {
   {
     int8_t sign = 0;
     int8_t flag = 0;
+    int8_t DecayChannel = 0;
+    std::vector<int> arrDaughIndex;
+    std::array<int, 2> arrPDGDaugh;
+    std::array<int, 2> arrPDGResonant1 = {2212, 313}; // Λc± → p± K*
+    std::array<int, 2> arrPDGResonant2 = {2224, 321}; // Λc± → Δ(1232)±± K∓
+    std::array<int, 2> arrPDGResonant3 = {3124, 211}; // Λc± → Λ(1520) π±
 
     // Match reconstructed candidates.
     for (auto& candidate : candidates) {
       //Printf("New rec. candidate");
       flag = 0;
+      DecayChannel = 0;
+      arrDaughIndex.clear();
       auto arrayDaughters = array{candidate.index0_as<aod::BigTracksMC>(), candidate.index1_as<aod::BigTracksMC>(), candidate.index2_as<aod::BigTracksMC>()};
 
       // D± → π± K∓ π±
@@ -168,18 +176,37 @@ struct HFCandidateCreator3ProngMC {
       // Λc± → p± K∓ π±
       if (flag == 0) {
         //Printf("Checking Λc± → p± K∓ π±");
-        if (RecoDecay::getMatchedMCRec(particlesMC, std::move(arrayDaughters), 4122, array{+kProton, -kKPlus, +kPiPlus}, true, &sign) > -1) {
+        auto indexRecLc = RecoDecay::getMatchedMCRec(particlesMC, std::move(arrayDaughters), 4122, array{+kProton, -kKPlus, +kPiPlus}, true, &sign, 2);
+        if (indexRecLc > -1) {
           flag = sign * (1 << LcToPKPi);
+
+          //Printf("Flagging the different Λc± → p± K∓ π± decay channels");
+          RecoDecay::getDaughters(particlesMC, indexRecLc, &arrDaughIndex, array{0}, 1);
+          if (arrDaughIndex.size() == 2) {
+            for (auto iProng = 0; iProng < arrDaughIndex.size(); ++iProng) {
+              auto daughI = particlesMC.iteratorAt(arrDaughIndex[iProng]);
+              arrPDGDaugh[iProng] = std::abs(daughI.pdgCode());
+            }
+            if (arrPDGDaugh[0] == arrPDGResonant1[0] && arrPDGDaugh[1] == arrPDGResonant1[1]) {
+              DecayChannel = 1;
+            } else if (arrPDGDaugh[0] == arrPDGResonant2[0] && arrPDGDaugh[1] == arrPDGResonant2[1]) {
+              DecayChannel = 2;
+            } else if (arrPDGDaugh[0] == arrPDGResonant3[0] && arrPDGDaugh[1] == arrPDGResonant3[1]) {
+              DecayChannel = 3;
+            }
+          }
         }
       }
 
-      rowMCMatchRec(flag);
+      rowMCMatchRec(flag, DecayChannel);
     }
 
     // Match generated particles.
     for (auto& particle : particlesMC) {
       //Printf("New gen. candidate");
       flag = 0;
+      DecayChannel = 0;
+      arrDaughIndex.clear();
 
       // D± → π± K∓ π±
       //Printf("Checking D± → π± K∓ π±");
@@ -190,12 +217,29 @@ struct HFCandidateCreator3ProngMC {
       // Λc± → p± K∓ π±
       if (flag == 0) {
         //Printf("Checking Λc± → p± K∓ π±");
-        if (RecoDecay::isMatchedMCGen(particlesMC, particle, 4122, array{+kProton, -kKPlus, +kPiPlus}, true, &sign)) {
+        auto isMatchedGenLc = RecoDecay::isMatchedMCGen(particlesMC, particle, 4122, array{+kProton, -kKPlus, +kPiPlus}, true, &sign, 2);
+        if (isMatchedGenLc) {
           flag = sign * (1 << LcToPKPi);
+
+          //Printf("Flagging the different Λc± → p± K∓ π± decay channels");
+          RecoDecay::getDaughters(particlesMC, particle.globalIndex(), &arrDaughIndex, array{0}, 1);
+          if (arrDaughIndex.size() == 2) {
+            for (auto jProng = 0; jProng < arrDaughIndex.size(); ++jProng) {
+              auto daughJ = particlesMC.iteratorAt(arrDaughIndex[jProng]);
+              arrPDGDaugh[jProng] = std::abs(daughJ.pdgCode());
+            }
+            if (arrPDGDaugh[0] == arrPDGResonant1[0] && arrPDGDaugh[1] == arrPDGResonant1[1]) {
+              DecayChannel = 1;
+            } else if (arrPDGDaugh[0] == arrPDGResonant2[0] && arrPDGDaugh[1] == arrPDGResonant2[1]) {
+              DecayChannel = 2;
+            } else if (arrPDGDaugh[0] == arrPDGResonant3[0] && arrPDGDaugh[1] == arrPDGResonant3[1]) {
+              DecayChannel = 3;
+            }
+          }
         }
       }
 
-      rowMCMatchGen(flag);
+      rowMCMatchGen(flag, DecayChannel);
     }
   }
 };