[PWGCF] Femto: add Lc support to the femto framework

PWGCF/Femto/Core/closePairRejection.h

@@ -84,6 +84,9 @@ struct ConfCpr : o2::framework::ConfigurableGroup {
 constexpr const char PrefixCprTrackTrack[] = "CprTrackTrack";
 constexpr const char PrefixCprTrackV0Daughter[] = "CprTrackV0Daughter";
 constexpr const char PrefixCprTrackD0Daughter[] = "CprTrackD0Daughter";
+constexpr const char PrefixCprTrackLcProton[] = "CprTrackLcProton";
+constexpr const char PrefixCprTrackLcKaon[] = "CprTrackLcKaon";
+constexpr const char PrefixCprTrackLcPion[] = "CprTrackLcPion";
 constexpr const char PrefixCprD0DaughterD0DaughterPos[] = "CprD0DaughterD0DaughterPos";
 constexpr const char PrefixCprD0DaughterD0DaughterNeg[] = "CprD0DaughterD0DaughterNeg";
 constexpr const char PrefixCprTrackResonanceDaughter[] = "CprTrackResonanceDaughter";
@@ -98,6 +101,9 @@ constexpr const char PrefixCprTrackCascadeBachelor[] = "CprTrackCascadeBachelor"
 using ConfCprTrackTrack = ConfCpr<PrefixCprTrackTrack>;
 using ConfCprTrackV0Daughter = ConfCpr<PrefixCprTrackV0Daughter>;
 using ConfCprTrackD0Daughter = ConfCpr<PrefixCprTrackD0Daughter>;
+using ConfCprTrackLcProton = ConfCpr<PrefixCprTrackLcProton>;
+using ConfCprTrackLcKaon = ConfCpr<PrefixCprTrackLcKaon>;
+using ConfCprTrackLcPion = ConfCpr<PrefixCprTrackLcPion>;
 using ConfCprD0DaugherD0DaughterPos = ConfCpr<PrefixCprD0DaughterD0DaughterPos>;
 using ConfCprD0DaugherD0DaughterNeg = ConfCpr<PrefixCprD0DaughterD0DaughterNeg>;
 using ConfCprTrackResonanceDaughter = ConfCpr<PrefixCprTrackResonanceDaughter>;
@@ -123,6 +129,12 @@ constexpr char PrefixD0D0PosSe[] = "CPR_D0D0_PosDau/SE/";
 constexpr char PrefixD0D0NegSe[] = "CPR_D0D0_NegDau/SE/";
 constexpr char PrefixD0D0PosMe[] = "CPR_D0D0_PosDau/ME/";
 constexpr char PrefixD0D0NegMe[] = "CPR_D0D0_NegDau/ME/";
+constexpr char PrefixTrackLcProtonSe[] = "CPR_TrackLc_Proton/SE/";
+constexpr char PrefixTrackLcKaonSe[] = "CPR_TrackLc_Kaon/SE/";
+constexpr char PrefixTrackLcPionSe[] = "CPR_TrackLc_Pion/SE/";
+constexpr char PrefixTrackLcProtonMe[] = "CPR_TrackLc_Proton/ME/";
+constexpr char PrefixTrackLcKaonMe[] = "CPR_TrackLc_Kaon/ME/";
+constexpr char PrefixTrackLcPionMe[] = "CPR_TrackLc_Pion/ME/";
 constexpr char PrefixV0V0PosSe[] = "CPR_V0V0_PosDau/SE/";
 constexpr char PrefixV0V0NegSe[] = "CPR_V0V0_NegDau/SE/";
 constexpr char PrefixV0V0PosMe[] = "CPR_V0V0_PosDau/ME/";
@@ -539,6 +551,58 @@ class ClosePairRejectionTrackV0 // can also be used for any particle type that h
   CloseTrackRejection<prefixTrackV0> mCtr;
 };
 
+template <auto& prefixProton, auto& prefixKaon, auto& prefixPion>
+class ClosePairRejectionTrackLc
+{
+ public:
+  template <typename T1, typename T2, typename T3>
+  void init(o2::framework::HistogramRegistry* registry,
+            std::map<CprHist, std::vector<o2::framework::AxisSpec>> const& specsProton,
+            std::map<CprHist, std::vector<o2::framework::AxisSpec>> const& specsKaon,
+            std::map<CprHist, std::vector<o2::framework::AxisSpec>> const& specsPion,
+            T1 const& confCprProton,
+            T2 const& confCprKaon,
+            T3 const& confCprPion,
+            int absChargeTrack)
+  {
+    mCtrProton.init(registry, specsProton, confCprProton, absChargeTrack, 1);
+    mCtrKaon.init(registry, specsKaon, confCprKaon, absChargeTrack, 1);
+    mCtrPion.init(registry, specsPion, confCprPion, absChargeTrack, 1);
+  }
+
+  void setMagField(float magField)
+  {
+    mCtrProton.setMagField(magField);
+    mCtrKaon.setMagField(magField);
+    mCtrPion.setMagField(magField);
+  }
+
+  template <typename T1, typename T2, typename T3>
+  void setPair(T1 const& track, T2 const& lc, T3 const& trackTable)
+  {
+    auto prong0 = trackTable.rawIteratorAt(lc.prong0DauId() - trackTable.offset());
+    mCtrProton.compute(track, prong0);
+    auto prong1 = trackTable.rawIteratorAt(lc.prong1DauId() - trackTable.offset());
+    mCtrKaon.compute(track, prong1);
+    auto prong2 = trackTable.rawIteratorAt(lc.prong2DauId() - trackTable.offset());
+    mCtrPion.compute(track, prong2);
+  }
+
+  [[nodiscard]] bool isClosePair() const { return mCtrProton.isClosePair() || mCtrKaon.isClosePair() || mCtrPion.isClosePair(); }
+
+  void fill(float kstar)
+  {
+    mCtrProton.fill(kstar);
+    mCtrKaon.fill(kstar);
+    mCtrPion.fill(kstar);
+  }
+
+ private:
+  CloseTrackRejection<prefixProton> mCtrProton;
+  CloseTrackRejection<prefixKaon> mCtrKaon;
+  CloseTrackRejection<prefixPion> mCtrPion;
+};
+
 template <auto& prefixBachelor, auto& prefixV0Daughter>
 class ClosePairRejectionTrackCascade
 {

PWGCF/Femto/Core/mcBuilder.h

@@ -63,6 +63,7 @@ struct McBuilderProducts : o2::framework::ProducesGroup {
   o2::framework::Produces<o2::aod::FLambdaLabels> producedLambdaLabels;
   o2::framework::Produces<o2::aod::FK0shortLabels> producedK0shortLabels;
   o2::framework::Produces<o2::aod::FD0Labels> producedD0Labels;
+  o2::framework::Produces<o2::aod::FLcLabels> producedLcLabels;
   o2::framework::Produces<o2::aod::FSigmaLabels> producedSigmaLabels;
   o2::framework::Produces<o2::aod::FSigmaPlusLabels> producedSigmaPlusLabels;
   o2::framework::Produces<o2::aod::FXiLabels> producedXiLabels;
@@ -82,6 +83,7 @@ struct ConfMcTables : o2::framework::ConfigurableGroup {
   o2::framework::Configurable<int> producedLambdaLabels{"producedLambdaLabels", -1, "Produce lambda labels (-1: auto; 0 off; 1 on)"};
   o2::framework::Configurable<int> producedK0shortLabels{"producedK0shortLabels", -1, "Produce k0short labels (-1: auto; 0 off; 1 on)"};
   o2::framework::Configurable<int> producedD0Labels{"producedD0Labels", -1, "Produce D0 labels (-1: auto; 0 off; 1 on)"};
+  o2::framework::Configurable<int> producedLcLabels{"producedLcLabels", -1, "Produce Lc labels (-1: auto; 0 off; 1 on)"};
   o2::framework::Configurable<int> producedSigmaLabels{"producedSigmaLabels", -1, "Produce k0short labels (-1: auto; 0 off; 1 on)"};
   o2::framework::Configurable<int> producedSigmaPlusLabels{"producedSigmaPlusLabels", -1, "Produce k0short labels (-1: auto; 0 off; 1 on)"};
   o2::framework::Configurable<int> producedXiLabels{"producedXiLabels", -1, "Produce xi labels (-1: auto; 0 off; 1 on)"};
@@ -148,6 +150,7 @@ class McBuilder
     mProduceLambdaLabels = utils::enableTable("FLambdaLabels", table.producedLambdaLabels.value, initContext);
     mProduceK0shortLabels = utils::enableTable("FK0shortLabels", table.producedK0shortLabels.value, initContext);
     mProduceD0Labels = utils::enableTable("FD0Labels", table.producedD0Labels.value, initContext);
+    mProduceLcLabels = utils::enableTable("FLcLabels", table.producedLcLabels.value, initContext);
     mProduceSigmaLabels = utils::enableTable("FSigmaLabels", table.producedSigmaLabels.value, initContext);
     mProduceSigmaPlusLabels = utils::enableTable("FSigmaPlusLabels", table.producedSigmaPlusLabels.value, initContext);
     mProduceXiLabels = utils::enableTable("FXiLabels", table.producedXiLabels.value, initContext);
@@ -160,7 +163,7 @@ class McBuilder
         mProduceLambdaLabels || mProduceK0shortLabels ||
         mProduceSigmaLabels || mProduceSigmaPlusLabels ||
         mProduceXiLabels || mProduceOmegaLabels ||
-        mProduceD0Labels) {
+        mProduceD0Labels || mProduceLcLabels) {
       mFillAnyTable = true;
     } else {
       LOG(info) << "No tables configured...";
@@ -257,7 +260,8 @@ class McBuilder
   {
     // charm hadrons get a prompt/non-prompt origin, consistent with the reco-matched path;
     // all other particles use the generic getOrigin inside getOrCreateMcParticleRow
-    if (std::abs(mcParticle.pdgCode()) == o2::constants::physics::Pdg::kD0) {
+    const int pdgAbs = std::abs(mcParticle.pdgCode());
+    if (pdgAbs == o2::constants::physics::Pdg::kD0) {
       // truth-level acceptance for the efficiency denominator: keep only
       // generated D0 -> K pi decays inside the rapidity acceptance
       int8_t sign = 0;
@@ -268,6 +272,15 @@ class McBuilder
         return;
       }
     }
+    if (pdgAbs == o2::constants::physics::Pdg::kLambdaCPlus) {
+      int8_t sign = 0;
+      if (!RecoDecay::isMatchedMCGen(mcParticles, mcParticle, o2::constants::physics::Pdg::kLambdaCPlus, std::array{+kProton, -kKPlus, +kPiPlus}, true, &sign, 2)) {
+        return;
+      }
+      if (std::abs(mcParticle.y()) > mCharmYGenMax) {
+        return;
+      }
+    }
     this->getOrCreateMcParticleRow<system>(mcParticle, mcParticles, mcCol, mcProducts);
   }
 
@@ -357,6 +370,35 @@ class McBuilder
     mcProducts.producedD0Labels(mcParticleRow);
   }
 
+  // Lc has no direct MC label. The three prongs are matched to a generated
+  // Lc -> p K pi decay;depthMax = 2 because the decay also proceeds through intermediate resonances.
+  template <modes::System system, typename T1, typename T2, typename T3, typename T4, typename T5>
+  void fillMcLcWithLabel(T1 const& lcCandidate, T2 const& /*tracks*/, T3 const& mcParticles, T4 const& /*mcCols*/, T5& mcProducts)
+  {
+    if (!mProduceLcLabels) {
+      mcProducts.producedLcLabels(-1);
+      return;
+    }
+
+    auto prong0 = lcCandidate.template prong0_as<T2>();
+    auto prong1 = lcCandidate.template prong1_as<T2>();
+    auto prong2 = lcCandidate.template prong2_as<T2>();
+    auto arrayDaughters = std::array{prong0, prong1, prong2};
+    int8_t sign = 0;
+    const int indexMcRec = RecoDecay::getMatchedMCRec(mcParticles, arrayDaughters, o2::constants::physics::Pdg::kLambdaCPlus, std::array{+kProton, -kKPlus, +kPiPlus}, true, &sign, 2);
+
+    if (indexMcRec < 0) {
+      mcProducts.producedLcLabels(-1);
+      return;
+    }
+
+    auto mcParticle = mcParticles.rawIteratorAt(indexMcRec);
+    auto mcCol = mcParticle.template mcCollision_as<T4>();
+    int64_t mcParticleRow = this->getOrCreateMcParticleRow<system>(mcParticle, mcParticles, mcCol, mcProducts);
+
+    mcProducts.producedLcLabels(mcParticleRow);
+  }
+
   template <modes::System system, typename T1, typename T2, typename T3, typename T4>
   void fillMcSigmaWithLabel(T1 const& sigmaDaughter, T2 const& mcParticles, T3 const& mcCols, T4& mcProducts)
   {
@@ -467,7 +509,10 @@ class McBuilder
   template <modes::System system, typename T1, typename T2, typename T3, typename T4>
   int64_t getOrCreateMcParticleRow(T1 const& mcParticle, T2 const& mcParticles, T3 const& mcCol, T4& mcProducts)
   {
-    auto origin = std::abs(mcParticle.pdgCode()) == o2::constants::physics::Pdg::kD0
+    const int pdgAbs = std::abs(mcParticle.pdgCode());
+    const bool isCharmHadron = (pdgAbs == o2::constants::physics::Pdg::kD0 ||
+                                pdgAbs == o2::constants::physics::Pdg::kLambdaCPlus);
+    auto origin = isCharmHadron
                     ? this->getHeavyFlavourOrigin(mcParticle, mcParticles)
                     : this->getOrigin(mcParticle);
     return this->buildMcParticleRow<system>(mcParticle, mcParticles, mcCol, origin, mcProducts);
@@ -686,6 +731,7 @@ class McBuilder
   bool mProduceLambdaLabels = false;
   bool mProduceK0shortLabels = false;
   bool mProduceD0Labels = false;
+  bool mProduceLcLabels = false;
   bool mProduceSigmaLabels = false;
   bool mProduceSigmaPlusLabels = false;
   bool mProduceXiLabels = false;

PWGCF/Femto/Core/modes.h

@@ -174,8 +174,8 @@ enum class TwoTrackResonance : o2::analysis::femto::datatypes::TwoTrackResonance
 enum class CharmHadron : o2::analysis::femto::datatypes::CharmHadronType {
   kD0,
   kD0Bar,
-  kDplus,
-  kLc
+  kLc,
+  kLcBar
 };
 
 enum class QvecDetector : o2::analysis::femto::datatypes::QvecDetectorType {