[PWGEM] Differentiate between b2c2e and b2cc2e, add function hasNCharmHadrons… (#14941)

Co-authored-by: Daiki Sekihata <daiki.sekihata@cern.ch>

Author: feisenhu

PWGEM/Dilepton/Tasks/checkMCPairTemplate.cxx

@@ -309,7 +309,7 @@ struct checkMCPairTemplate {
   HistogramRegistry fRegistry{"output", {}, OutputObjHandlingPolicy::AnalysisObject, false, false};
   static constexpr std::string_view event_cut_types[2] = {"before/", "after/"};
   static constexpr std::string_view pair_sign_types[3] = {"uls/", "lspp/", "lsmm/"};
-  static constexpr std::string_view dilepton_source_types[20] = {
+  static constexpr std::string_view dilepton_source_types[22] = {
     "sm/Photon/",             // 0
     "sm/PromptPi0/",          // 1
     "sm/NonPromptPi0/",       // 2
@@ -329,7 +329,9 @@ struct checkMCPairTemplate {
     "bbbar/b2l_b2l/",         // 16
     "bbbar/b2c2l_b2c2l/",     // 17
     "bbbar/b2c2l_b2l_sameb/", // 18
-    "bbbar/b2c2l_b2l_diffb/"  // 19
+    "bbbar/b2c2l_b2l_diffb/", // 19
+    "bbbar/b2cc2l_b2c2l/",    // 20
+    "bbbar/b2cc2l_b2cc2l/",   // 21
   }; // unordered_map is better, but cannot be constexpr.
   static constexpr std::string_view unfolding_dilepton_source_types[3] = {"sm/", "ccbar/", "bbbar/"};
 
@@ -424,6 +426,8 @@ struct checkMCPairTemplate {
     fRegistry.addClone("Generated/ccbar/c2l_c2l/", "Generated/bbbar/b2c2l_b2c2l/");
     fRegistry.addClone("Generated/ccbar/c2l_c2l/", "Generated/bbbar/b2c2l_b2l_sameb/");
     fRegistry.addClone("Generated/ccbar/c2l_c2l/", "Generated/bbbar/b2c2l_b2l_diffb/"); // LS
+    fRegistry.addClone("Generated/ccbar/c2l_c2l/", "Generated/bbbar/b2cc2l_b2c2l/");
+    fRegistry.addClone("Generated/ccbar/c2l_c2l/", "Generated/bbbar/b2cc2l_b2cc2l/");
 
     // for charmed hadrons // create 28 combinations
     static constexpr std::string_view charmed_mesons[] = {"Dplus", "D0", "Dsplus"}; // 411, 421, 431
@@ -562,6 +566,8 @@ struct checkMCPairTemplate {
     fRegistry.addClone("Pair/ccbar/c2l_c2l/", "Pair/bbbar/b2c2l_b2c2l/");
     fRegistry.addClone("Pair/ccbar/c2l_c2l/", "Pair/bbbar/b2c2l_b2l_sameb/");
     fRegistry.addClone("Pair/ccbar/c2l_c2l/", "Pair/bbbar/b2c2l_b2l_diffb/"); // LS
+    fRegistry.addClone("Pair/ccbar/c2l_c2l/", "Pair/bbbar/b2cc2l_b2c2l/");
+    fRegistry.addClone("Pair/ccbar/c2l_c2l/", "Pair/bbbar/b2cc2l_b2cc2l/");
 
     if (cfgFillSeparateCharmHadronPairs) {
       for (int im = 0; im < nm_c; im++) {
@@ -1949,10 +1955,10 @@ struct checkMCPairTemplate {
     // o2::math_utils::bringToPMPi(phiPol);
     // float quadmom = (3.f * std::pow(cos_thetaPol, 2) - 1.f) / 2.f;
 
-    if ((FindCommonMotherFrom2ProngsWithoutPDG(t1mc, t2mc) > 0 || IsHF(t1mc, t2mc, mcparticles) > 0) && is_pair_from_same_mcevent) { // for bkg study
-      if (std::abs(t1mc.pdgCode()) != pdg_lepton || std::abs(t2mc.pdgCode()) != pdg_lepton) {                                        // hh or lh correlated bkg
-        if (std::abs(t1mc.pdgCode()) != pdg_lepton && std::abs(t2mc.pdgCode()) != pdg_lepton) {                                      // hh correlated bkg
-          if (t1.sign() * t2.sign() < 0) {                                                                                           // ULS
+    if ((FindCommonMotherFrom2ProngsWithoutPDG(t1mc, t2mc) > 0 || IsHF<true>(t1mc, t2mc, mcparticles) > 0) && is_pair_from_same_mcevent) { // for bkg study
+      if (std::abs(t1mc.pdgCode()) != pdg_lepton || std::abs(t2mc.pdgCode()) != pdg_lepton) {                                              // hh or lh correlated bkg
+        if (std::abs(t1mc.pdgCode()) != pdg_lepton && std::abs(t2mc.pdgCode()) != pdg_lepton) {                                            // hh correlated bkg
+          if (t1.sign() * t2.sign() < 0) {                                                                                                 // ULS
             fRegistry.fill(HIST("Pair/corr_bkg_hh/uls/hs"), v12.M(), v12.Pt(), pair_dca, weight);
           } else if (t1.sign() > 0 && t2.sign() > 0) { // LS++
             fRegistry.fill(HIST("Pair/corr_bkg_hh/lspp/hs"), v12.M(), v12.Pt(), pair_dca, weight);
@@ -1990,7 +1996,7 @@ struct checkMCPairTemplate {
       return false;
     }
     int mother_id = std::max({FindSMULS(t1mc, t2mc, mcparticles), FindSMULS(t2mc, t1mc, mcparticles), FindSMLSPP(t1mc, t2mc, mcparticles), FindSMLSMM(t1mc, t2mc, mcparticles)});
-    int hfee_type = IsHF(t1mc, t2mc, mcparticles);
+    int hfee_type = IsHF<true>(t1mc, t2mc, mcparticles);
     if (mother_id < 0 && hfee_type < 0) {
       return false;
     }
@@ -2159,7 +2165,7 @@ struct checkMCPairTemplate {
     }
 
     int mother_id = std::max({FindSMULS(t1, t2, mcparticles), FindSMULS(t2, t1, mcparticles), FindSMLSPP(t1, t2, mcparticles), FindSMLSMM(t1, t2, mcparticles)});
-    int hfee_type = IsHF(t1, t2, mcparticles);
+    int hfee_type = IsHF<true>(t1, t2, mcparticles);
     if (mother_id < 0 && hfee_type < 0) {
       return false;
     }
@@ -2326,6 +2332,12 @@ struct checkMCPairTemplate {
         case static_cast<int>(EM_HFeeType::kBCe_Be_DiffB):
           fillGenHistograms<19>(sign1, sign2, mp1.pdgCode(), mp2.pdgCode(), v12.M(), v12.Pt(), weight); // b2c2l_b2l_diffb
           break;
+        case static_cast<int>(EM_HFeeType::kBCCe_BCe):
+          fillGenHistograms<20>(sign1, sign2, mp1.pdgCode(), mp2.pdgCode(), v12.M(), v12.Pt(), weight); // b2cc2l_b2c2l
+          break;
+        case static_cast<int>(EM_HFeeType::kBCCe_BCCe):
+          fillGenHistograms<21>(sign1, sign2, mp1.pdgCode(), mp2.pdgCode(), v12.M(), v12.Pt(), weight); // b2cc2l_b2cc2l
+          break;
         default:
           break;
       }
@@ -2791,7 +2803,7 @@ struct checkMCPairTemplate {
       return false;
     }
     int mother_id = std::max({FindSMULS(t1mc, t2mc, mcparticles), FindSMULS(t2mc, t1mc, mcparticles), FindSMLSPP(t1mc, t2mc, mcparticles), FindSMLSMM(t1mc, t2mc, mcparticles)});
-    int hfee_type = IsHF(t1mc, t2mc, mcparticles);
+    int hfee_type = IsHF<true>(t1mc, t2mc, mcparticles);
     if (mother_id < 0 && hfee_type < 0) {
       return false;
     }

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,56 @@ 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());
+    }
+  }
+  return -999;
+}
+//_______________________________________________________________________
 template <typename TMCParticle>
 bool isFlavorOscillationB(TMCParticle const& mcParticle)
 {
@@ -551,7 +603,7 @@ int find1stHadron(TMCParticle const& mcParticle, TMCParticles const& mcParticles
   return hadronId;
 }
 //_______________________________________________________________________
-template <typename TMCParticle1, typename TMCParticle2, typename TMCParticles>
+template <bool doMoreDifferentially = false, typename TMCParticle1, typename TMCParticle2, typename TMCParticles>
 int IsHF(TMCParticle1 const& p1, TMCParticle2 const& p2, TMCParticles const& mcparticles)
 {
   if (!p1.has_mothers() || !p2.has_mothers()) {
@@ -661,7 +713,23 @@ 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
+
+    if constexpr (!doMoreDifferentially) {
+      return static_cast<int>(EM_HFeeType::kBCe_BCe); // default to b->c->e and b->c->e, decay type = 1
+    } else {
+      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)) {