Add MC matching for c-deuteron

Author: fgrosa

PWGHF/Core/DecayChannels.h

@@ -92,8 +92,10 @@ enum DecayChannelMain : HfDecayChannel {
   XicToPKPi = 21,  // p K− π+
   XicToPKK = 22,   // p K− K+
   XicToSPiPi = 23, // Σ+ π− π+
+  // cd+
+  CDeuteronToDeKPi = 24, // de K− π+
   //
-  NChannelsMain = XicToSPiPi // last channel
+  NChannelsMain = CDeuteronToDeKPi // last channel
 };
 /// @brief 3-prong candidates: resonant channels
 enum DecayChannelResonant : HfDecayChannel {
@@ -131,8 +133,11 @@ enum DecayChannelResonant : HfDecayChannel {
   // Ξc+
   XicToPKstar0 = 27, // p anti-K*0(892)
   XicToPPhi = 28,    // p φ
-  //
-  NChannelsResonant = XicToPPhi // last channel
+  // cd+
+  CDeuteronToDeKstar0 = 30,
+  CDeuteronToNeDeltaplusK = 31,
+  CDeuteronToNeL1520Pi = 32,
+  NChannelsResonant = CDeuteronToNeL1520Pi // last channel
 };
 } // namespace hf_cand_3prong
 

PWGHF/TableProducer/candidateCreator3Prong.cxx

@@ -1327,6 +1327,25 @@ struct HfCandidateCreator3ProngExpressions {
             }
           }
         }
+
+        // cd± → de± K∓ π±
+        if (flagChannelMain == 0) {
+          auto arrPdgDaughtersCDeuteronToDeKPi{std::array{+Pdg::kDeuteron, -kKPlus, +kPiPlus}};
+          if (matchKinkedDecayTopology && matchInteractionsWithMaterial) {
+            indexRec = RecoDecay::getMatchedMCRec<false, false, false, true, true>(mcParticles, arrayDaughters, Pdg::kCDeuteron, arrPdgDaughtersCDeuteronToDeKPi, true, &sign, 1, &nKinkedTracks, &nInteractionsWithMaterial);
+          } else if (matchKinkedDecayTopology && !matchInteractionsWithMaterial) {
+            indexRec = RecoDecay::getMatchedMCRec<false, false, false, true, false>(mcParticles, arrayDaughters, Pdg::kCDeuteron, arrPdgDaughtersCDeuteronToDeKPi, true, &sign, 1, &nKinkedTracks);
+          } else if (!matchKinkedDecayTopology && matchInteractionsWithMaterial) {
+            indexRec = RecoDecay::getMatchedMCRec<false, false, false, false, true>(mcParticles, arrayDaughters, Pdg::kCDeuteron, arrPdgDaughtersCDeuteronToDeKPi, true, &sign, 1, nullptr, &nInteractionsWithMaterial);
+          } else {
+            indexRec = RecoDecay::getMatchedMCRec(mcParticles, arrayDaughters, Pdg::kCDeuteron, arrPdgDaughtersCDeuteronToDeKPi, true, &sign, 1);
+          }
+          if (indexRec > -1) {
+            flagChannelMain = static_cast<int8_t>(sign * DecayChannelMain::CDeuteronToDeKPi);
+            auto particle = mcParticles.rawIteratorAt(indexRec);
+            flagChannelResonant = hf_decay::getResonantDecayCDeuteron(particle);
+          }
+        }
       }
 
       // Check whether the particle is non-prompt (from a b quark).

PWGHF/Utils/utilsMcGen.h

@@ -301,6 +301,14 @@ void fillMcMatchGen3Prong(TMcParticles const& mcParticles,
           flagChannelMain = sign * DecayChannelMain::XicToPKPi;
         }
       }
+
+      // cd± → de± K∓ π±
+      if (flagChannelMain == 0) {
+        if (RecoDecay::isMatchedMCGen(mcParticles, particle, Pdg::kCDeuteron, std::array{+Pdg::kDeuteron, -kKPlus, +kPiPlus}, true, &sign, 1)) {
+          flagChannelMain = sign * DecayChannelMain::CDeuteronToDeKPi;
+          flagChannelResonant = o2::hf_decay::getResonantDecayCDeuteron(particle);
+        }
+      }
     }
 
     // Check whether the particle is non-prompt (from a b quark).

PWGHF/Utils/utilsMcMatching.h

@@ -154,6 +154,17 @@ static const std::unordered_map<DecayChannelResonant, const std::array<int, 2>>
   {DecayChannelResonant::XicToPPhi, {+PDG_t::kProton, +o2::constants::physics::Pdg::kPhi}},
 };
 
+// cd+
+
+static const std::unordered_map<DecayChannelMain, const std::vector<int>> daughtersCDeuteronMain{
+  {DecayChannelMain::CDeuteronToDeKPi, {+o2::constants::physics::Pdg::kDeuteron, +PDG_t::kKMinus, +PDG_t::kPiPlus}}
+};
+
+/// resonances in c-deuteron decay are not stored in the particle stack for c-deuteron, but tagged with specific status codes
+static constexpr int statusCodeCDeuteronToDeKstar0{95};
+static constexpr int statusCodeCDeuteronToNeDeltaplusK{96};
+static constexpr int statusCodeCDeuteronToNeL1520Pi{97};
+
 /// Returns a map of the possible final states for a specific 3-prong particle specie
 /// \param pdgMother PDG code of the mother particle
 /// \return a map of final states with their corresponding PDG codes
@@ -170,6 +181,8 @@ inline std::unordered_map<DecayChannelMain, const std::vector<int>> getDecayChan
       return daughtersLcMain;
     case o2::constants::physics::Pdg::kXiCPlus:
       return daughtersXicMain;
+    case o2::constants::physics::Pdg::kCDeuteron:
+      return daughtersCDeuteronMain;
     default:
       LOG(fatal) << "Unknown PDG code for 3-prong final states: " << pdgMother;
       return {};
@@ -317,6 +330,27 @@ inline void flipPdgSign(const int pdgMother, const int pdgToFlip, std::array<int
     }
   }
 }
+/// Get resonant channel for c-deuteron
+/// resonances are not stored in the particle stack for c-deuteron, but tagged with specific status codes
+/// \tparam particle is the c-deuteron
+/// \param pdgMother PDG code of the mother particle
+/// \param pdgToFlip PDG code to be flipped
+/// \param arrPdg array of PDG codes to be modified
+template <typename Part>
+inline int getResonantDecayCDeuteron(Part const& particle) 
+{
+  auto statusCode = std::abs(particle.getGenStatusCode());
+  if (statusCode == o2::hf_decay::hf_cand_3prong::statusCodeCDeuteronToDeKstar0) {
+    return o2::hf_decay::hf_cand_3prong::DecayChannelResonant::CDeuteronToDeKstar0;
+  }
+  if (statusCode == o2::hf_decay::hf_cand_3prong::statusCodeCDeuteronToNeDeltaplusK) {
+    return o2::hf_decay::hf_cand_3prong::DecayChannelResonant::CDeuteronToNeDeltaplusK;
+  }
+  if (statusCode == o2::hf_decay::hf_cand_3prong::statusCodeCDeuteronToNeL1520Pi) {
+    return o2::hf_decay::hf_cand_3prong::DecayChannelResonant::CDeuteronToNeL1520Pi;
+  }
+  return 0;
+}
 } // namespace o2::hf_decay
 
 #endif // PWGHF_UTILS_UTILSMCMATCHING_H_