Differentiate between b2c2e and b2cc2e, add function hasNCharmHadronsInBeautyDecay()

Author: feisenhu

PWGEM/Dilepton/Utils/MCUtilities.h

@@ -32,6 +32,8 @@ enum class EM_HFeeType : int {
   kBCe_BCe = 2,      // ULS
   kBCe_Be_SameB = 3, // ULS
   kBCe_Be_DiffB = 4, // LS
+  kBCCe_BCe = 5,     // ULS
+  kBCCe_BCCe = 6,    // ULS
 };
 
 //_______________________________________________________________________
@@ -461,6 +463,55 @@ int IsFromCharm(TMCParticle const& p, TMCParticles const& mcparticles)
   return -999;
 }
 //_______________________________________________________________________
+template <typename T, typename U>
+int hasNCharmHadronsInBeautyDecay(T const& mcParticle, U const& mcParticles)
+{
+  // require that the direct mother is beauty hadron via semileptonice decay. e.g. hb->e, not hb->X->pi0->eegamma
+  if (!mcParticle.has_mothers()) {
+    return -999;
+  }
+  if (!IsFromBeauty(mcParticle, mcParticles)) {
+    return -999;
+  }
+
+  auto mp = mcParticles.iteratorAt(mcParticle.mothersIds()[0]);
+  if (!isCharmMeson(mp) && !isCharmBaryon(mp)) {
+    return -999;
+  }
+
+  int motherid = mcParticle.mothersIds()[0]; // first mother index
+  while (motherid > -1) {
+    if (motherid < mcParticles.size()) { // protect against bad mother indices. why is this needed?
+      mp = mcParticles.iteratorAt(motherid);
+      if (std::abs(mp.pdgCode()) < 1e+9 && (std::to_string(std::abs(mp.pdgCode()))[std::to_string(std::abs(mp.pdgCode())).length() - 3] == '5' || std::to_string(std::abs(mp.pdgCode()))[std::to_string(std::abs(mp.pdgCode())).length() - 4] == '5')) {
+        // check if mp has two charm hadrons as daughters
+        if (mp.has_daughters()) {
+          const auto& daughtersIds = mp.daughtersIds();
+          int count_charm_hadron = 0;
+          for (const auto& daughterId : daughtersIds) {
+            if (daughterId >= 0 && daughterId < mcParticles.size()) {
+              auto daughter = mcParticles.iteratorAt(daughterId);
+              if (isCharmMeson(daughter) || isCharmBaryon(daughter)) {
+                count_charm_hadron++;
+              }
+            }
+          }
+          return count_charm_hadron;
+        } else {
+          LOGF(debug, "Something went wrong: Did not find any daughter for the current mother! Can't be a mother if there are no daughters\n");
+        }
+      }
+      if (mp.has_mothers()) {
+        motherid = mp.mothersIds()[0];
+      } else {
+        return -999;
+      }
+    } else {
+      LOGF(info, "Mother label(%d) exceeds the McParticles size(%d)", motherid, mcParticles.size());
+    }
+  }
+}
+//_______________________________________________________________________
 template <typename TMCParticle>
 bool isFlavorOscillationB(TMCParticle const& mcParticle)
 {
@@ -661,7 +712,18 @@ int IsHF(TMCParticle1 const& p1, TMCParticle2 const& p2, TMCParticles const& mcp
     mothers_pdg1.shrink_to_fit();
     mothers_id2.shrink_to_fit();
     mothers_pdg2.shrink_to_fit();
-    return static_cast<int>(EM_HFeeType::kBCe_BCe); // b->c->e and b->c->e, decay type = 1
+    int n_c_from_b1 = hasNCharmHadronsInBeautyDecay(p1, mcparticles);
+    int n_c_from_b2 = hasNCharmHadronsInBeautyDecay(p2, mcparticles);
+    if (n_c_from_b1 == 1 && n_c_from_b2 == 1) {
+      return static_cast<int>(EM_HFeeType::kBCe_BCe); // b->c->e and b->c->e, decay type = 1
+    } else if (n_c_from_b1 == 2 && n_c_from_b2 == 2) {
+      return static_cast<int>(EM_HFeeType::kBCCe_BCCe); // b->cc->e and b->cc->e, decay type = 6
+    } else if ((n_c_from_b1 == 1 && n_c_from_b2 == 2) || (n_c_from_b1 == 2 && n_c_from_b2 == 1)) {
+      return static_cast<int>(EM_HFeeType::kBCCe_BCe); // b->cc->e and b->c->e, decay type = 5
+    } else {
+      LOGF(debug, "Unexpected number of charm hadrons from beauty decay: n_c_from_b1 = %d, n_c_from_b2 = %d. Return kBCe_BCe as default.", n_c_from_b1, n_c_from_b2);
+      return static_cast<int>(EM_HFeeType::kBCe_BCe); // default to b->c->e and b->c->e, decay type = 1
+    }
   }
 
   if ((is_direct_from_b1 && is_c_from_b2) || (is_direct_from_b2 && is_c_from_b1)) {