[PWGHF] align the table between reduced and not-reduced workflow for Lb

PWGHF/D2H/DataModel/ReducedDataModel.h

@@ -506,15 +506,15 @@ using HfRedCandBs = soa::Join<HfCandBsExt, HfRedBsProngs>;
 
 namespace hf_cand_lb_reduced
 {
-DECLARE_SOA_INDEX_COLUMN_FULL(Prong0Lc, prong0Lc, int, HfRed3Prongs, "_0");          //! Prong0 index
-DECLARE_SOA_INDEX_COLUMN_FULL(Prong1Track, prong1Track, int, HfRedTrackBases, "_1"); //! Prong1 index
+DECLARE_SOA_INDEX_COLUMN_FULL(Prong0, prong0, int, HfRed3Prongs, "_0");              //! Prong0 index
+DECLARE_SOA_INDEX_COLUMN_FULL(Prong1, prong1, int, HfRedTrackBases, "_1");           //! Prong1 index
 DECLARE_SOA_COLUMN(Prong0MlScoreBkg, prong0MlScoreBkg, float);                       //! Bkg ML score of the Lc daughter
 DECLARE_SOA_COLUMN(Prong0MlScorePrompt, prong0MlScorePrompt, float);                 //! Prompt ML score of the Lc daughter
 DECLARE_SOA_COLUMN(Prong0MlScoreNonprompt, prong0MlScoreNonprompt, float);           //! Nonprompt ML score of the Lc daughter
 } // namespace hf_cand_lb_reduced
 
 DECLARE_SOA_TABLE(HfRedLbProngs, "AOD", "HFREDLBPRONG", //! Table with Lb daughter indices
-                  hf_cand_lb_reduced::Prong0LcId, hf_cand_lb_reduced::Prong1TrackId);
+                  hf_cand_lb_reduced::Prong0Id, hf_cand_lb_reduced::Prong1Id);
 
 DECLARE_SOA_TABLE(HfRedLbLcMls, "AOD", "HFREDLBLCML", //! Table with ML scores for the Lc daughter
                   hf_cand_lb_reduced::Prong0MlScoreBkg,
@@ -796,8 +796,8 @@ DECLARE_SOA_COLUMN(PdgCodeProng3, pdgCodeProng3, int);             //! Pdg code
 
 // table with results of reconstruction level MC matching
 DECLARE_SOA_TABLE(HfMcRecRedLcPis, "AOD", "HFMCRECREDLCPI", //! Table with reconstructed MC information on LcPi(<-Lb) pairs for reduced workflow
-                  hf_cand_lb_reduced::Prong0LcId,
-                  hf_cand_lb_reduced::Prong1TrackId,
+                  hf_cand_lb_reduced::Prong0Id,
+                  hf_cand_lb_reduced::Prong1Id,
                   hf_cand_lb::FlagMcMatchRec,
                   hf_cand_lb::FlagWrongCollision,
                   hf_cand_lb::DebugMcRec,

PWGHF/D2H/TableProducer/candidateCreatorLbReduced.cxx

@@ -207,7 +207,7 @@ struct HfCandidateCreatorLbReduced {
                          pVecLc[0], pVecLc[1], pVecLc[2],
                          pVecPion[0], pVecPion[1], pVecPion[2],
                          dcaLc.getY(), dcaPion.getY(),
-                         std::sqrt(dcaLc.getSigmaY2()), std::sqrt(dcaPion.getSigmaY2()), candLc.globalIndex(), trackPion.globalIndex(), -999.);
+                         std::sqrt(dcaLc.getSigmaY2()), std::sqrt(dcaPion.getSigmaY2()));
 
         rowCandidateProngs(candLc.globalIndex(), trackPion.globalIndex());
 
@@ -302,7 +302,7 @@ struct HfCandidateCreatorLbReducedExpressions {
     for (const auto& candLb : candsLb) {
       bool filledMcInfo{false};
       for (const auto& rowLcPiMcRec : rowsLcPiMcRec) {
-        if ((rowLcPiMcRec.prong0LcId() != candLb.prong0LcId()) || (rowLcPiMcRec.prong1TrackId() != candLb.prong1TrackId())) {
+        if ((rowLcPiMcRec.prong0Id() != candLb.prong0Id()) || (rowLcPiMcRec.prong1Id() != candLb.prong1Id())) {
           continue;
         }
         rowLbMcRec(rowLcPiMcRec.flagMcMatchRec(), rowLcPiMcRec.flagWrongCollision(), rowLcPiMcRec.debugMcRec(), rowLcPiMcRec.ptMother());

PWGHF/D2H/TableProducer/candidateSelectorLbToLcPiReduced.cxx

@@ -182,7 +182,7 @@ struct HfCandidateSelectorLbToLcPiReduced {
       }
 
       // track-level PID selection
-      auto trackPi = hfCandLb.template prong1Track_as<TracksPion>();
+      auto trackPi = hfCandLb.template prong1_as<TracksPion>();
       if (pionPidMethod == PidMethod::TpcOrTof || pionPidMethod == PidMethod::TpcAndTof) {
         int pidTrackPi{TrackSelectorPID::Status::NotApplicable};
         if (pionPidMethod == PidMethod::TpcOrTof) {

PWGHF/D2H/Tasks/taskLb.cxx

@@ -16,6 +16,8 @@
 /// \author Martin Voelkl <martin.andreas.volkl@cern.ch>, University of Birmingham
 
 #include <vector>
+#include <array>
+#include <algorithm>
 
 #include "CommonConstants/PhysicsConstants.h"
 #include "Framework/AnalysisTask.h"
@@ -42,7 +44,7 @@ struct HfTaskLb {
   Configurable<int> selectionFlagLb{"selectionFlagLb", 0, "Selection Flag for Lb"};
   Configurable<float> yCandGenMax{"yCandGenMax", 0.5, "max. gen particle rapidity"};
   Configurable<float> yCandRecoMax{"yCandRecoMax", 0.8, "max. cand. rapidity"};
-  Configurable<float> DCALengthParameter{"DCALengthParameter", 0.02, "decay length for DCA"};
+  Configurable<float> lengthDCAParameter{"lengthDCAParameter", 0.02, "decay length for DCA"};
   Configurable<float> minLikelihoodRatio{"minLikelihoodRatio", 10., "min. likelihood ratio for combined DCAs"};
   Configurable<float> minLikelihoodRatioLc{"minLikelihoodRatioLc", 10., "min. likelihood ratio for Lc cross check"};
   Configurable<float> mDiffKStar892Max{"mDiffKStar892Max", 0.0473, "Accepted range around KStar mass peak"};
@@ -57,24 +59,24 @@ struct HfTaskLb {
   HfHelper hfHelper;
   Service<o2::framework::O2DatabasePDG> pdg;
 
-  Filter filterSelectCandidates = (aod::hf_sel_candidate_lb::isSelLbToLcPi >= selectionFlagLb);
-
   using TracksWExt = soa::Join<o2::aod::Tracks, o2::aod::TracksExtra, aod::TrackSelection, o2::aod::TrackSelectionExtension, aod::TracksPidPi, aod::PidTpcTofFullPi, aod::TracksPidKa, aod::PidTpcTofFullKa>;
   using TracksWExtMc = soa::Join<o2::aod::Tracks, o2::aod::TracksExtra, aod::TrackSelection, o2::aod::TrackSelectionExtension, aod::TracksPidPi, aod::PidTpcTofFullPi, aod::TracksPidKa, aod::PidTpcTofFullKa, McTrackLabels>;
 
-  PresliceUnsorted<aod::TracksWMc> McPartID = aod::mctracklabel::mcParticleId;
+  Filter filterSelectCandidates = (aod::hf_sel_candidate_lb::isSelLbToLcPi >= selectionFlagLb);
+
+  PresliceUnsorted<aod::TracksWMc> mcPartID = aod::mctracklabel::mcParticleId;
 
   bool passesImpactParameterResolution(float pT, float d0Resolution)
   {
     float expectedResolution(0.001 + 0.0052 * std::exp(-0.655 * pT));
     return (d0Resolution <= expectedResolution * 1.5);
   } // Compares to pT dependent cut on impact parameter resolution
 
-  float logLikelihoodRatioSingleTrackDCA(float DCA, float reso, float lengthParameter)
+  float logLikelihoodRatioSingleTrackDCA(float dca, float reso, float lengthParameter)
   {
     reso *= resoCorrectionFactor; // In case real resolution is worse
-    float numerator = 1. / lengthParameter * std::exp(-DCA / lengthParameter);
-    float denominator = (1. - largeLifetimeBG) * TMath::Gaus(DCA, 0., reso, true) + largeLifetimeBG / 0.2; // flat distribution to 2 mm
+    float numerator = 1. / lengthParameter * std::exp(-dca / lengthParameter);
+    float denominator = (1. - largeLifetimeBG) * TMath::Gaus(dca, 0., reso, true) + largeLifetimeBG / 0.2; // flat distribution to 2 mm
     return std::log(numerator / denominator);
   } // Creates the single track log likelihood assuming an exonential law for the secondaries
 
@@ -172,6 +174,8 @@ struct HfTaskLb {
   {
     float massKStar892 = 0.892;
     float massDelta1232 = 1.232;
+    std::array<float, 3> dca = {0.f, 0.f, 0.f};
+    std::array<float, 3> dcaResolution = {0.f, 0.f, 0.f};
 
     for (const auto& candidateLc : candidatesLc) {
       if (!candidateLc.isSelLcToPKPi() && !candidateLc.isSelLcToPiKP())
@@ -194,12 +198,26 @@ struct HfTaskLb {
         continue;
       if (!passesImpactParameterResolution(track2.pt(), reso2))
         continue;
-      float DCA0 = candidateLc.impactParameter0();
-      float DCA1 = candidateLc.impactParameter1();
-      float DCA2 = candidateLc.impactParameter2();
-      if (DCA0 > maximumImpactParameterForLambdaCCrossChecks || DCA1 > maximumImpactParameterForLambdaCCrossChecks || DCA2 > maximumImpactParameterForLambdaCCrossChecks)
+
+      dca = {
+        candidateLc.impactParameter0(),
+        candidateLc.impactParameter1(),
+        candidateLc.impactParameter2()};
+
+      bool exceedsMaxDca = std::any_of(dca.begin(), dca.end(), [&](float val) {
+        return val > maximumImpactParameterForLambdaCCrossChecks;
+      });
+
+      if (exceedsMaxDca) {
         continue;
-      float likelihoodRatio = logLikelihoodRatioSingleTrackDCA(DCA0, reso0, DCALengthParameter) + logLikelihoodRatioSingleTrackDCA(DCA1, reso1, DCALengthParameter) + logLikelihoodRatioSingleTrackDCA(DCA2, reso2, DCALengthParameter);
+      }
+      dcaResolution = {reso0, reso1, reso2};
+
+      float likelihoodRatio = 0.0f;
+      for (size_t i = 0; i < dca.size(); ++i) {
+        likelihoodRatio += logLikelihoodRatioSingleTrackDCA(dca[i], dcaResolution[i], lengthDCAParameter);
+      }
+
       registry.get<TH2>(HIST("hPtlogLikelihood"))->Fill(candidateLc.pt(), likelihoodRatio);
       if (likelihoodRatio < minLikelihoodRatioLc)
         continue;
@@ -255,24 +273,31 @@ struct HfTaskLb {
     } // Lambda_c candidates loop for cross checks
 
     for (const auto& candidate : candidates) {
-      if (!(candidate.hfflag() & 1 << hf_cand_lb::DecayType::LbToLcPi)) { // This should never be true as the loop is over Lb candidates
-        continue;
-      }
+
       if (yCandRecoMax >= 0. && std::abs(hfHelper.yLb(candidate)) > yCandRecoMax) {
         continue;
       }
       registry.get<TH1>(HIST("hZVertex"))->Fill(collision.posZ());
 
       auto candLc = candidate.prong0_as<soa::Join<aod::HfCand3Prong, aod::HfSelLc>>();
-      float d0resolution0 = candLc.errorImpactParameter0();
-      float d0resolution1 = candLc.errorImpactParameter1();
-      float d0resolution2 = candLc.errorImpactParameter2();
-      float DCA0 = candLc.impactParameter0();
-      float DCA1 = candLc.impactParameter1();
-      float DCA2 = candLc.impactParameter2();
-      float likelihoodRatio = logLikelihoodRatioSingleTrackDCA(DCA0, d0resolution0, DCALengthParameter) + logLikelihoodRatioSingleTrackDCA(DCA1, d0resolution1, DCALengthParameter) + logLikelihoodRatioSingleTrackDCA(DCA2, d0resolution2, DCALengthParameter);
-      if (likelihoodRatio < minLikelihoodRatio)
+      dca = {
+        candLc.impactParameter0(),
+        candLc.impactParameter1(),
+        candLc.impactParameter2()};
+
+      dcaResolution = {
+        candLc.errorImpactParameter0(),
+        candLc.errorImpactParameter1(),
+        candLc.errorImpactParameter2()};
+
+      float likelihoodRatio = 0.0f;
+      for (size_t i = 0; i < dca.size(); ++i) {
+        likelihoodRatio += logLikelihoodRatioSingleTrackDCA(dca[i], dcaResolution[i], lengthDCAParameter);
+      }
+
+      if (likelihoodRatio < minLikelihoodRatio) {
         continue; // Larger likelihood means more likely to be signal
+      }
       float lbMass = hfHelper.invMassLbToLcPi(candidate);
       registry.get<TH2>(HIST("hPtinvMassLb"))->Fill(candidate.pt(), lbMass);
 
@@ -304,15 +329,14 @@ struct HfTaskLb {
   {
     // MC rec
     for (const auto& candidate : candidates) {
-      if (!(candidate.hfflag() & 1 << hf_cand_lb::DecayType::LbToLcPi)) {
-        continue;
-      }
+
       if (yCandRecoMax >= 0. && std::abs(hfHelper.yLb(candidate)) > yCandRecoMax) {
         continue;
       }
-      auto candLc = candidate.prong0_as<aod::HfCand3Prong>();
+      auto candLc = candidate.prong0_as<soa::Join<aod::HfCand3Prong, aod::HfCand3ProngMcRec>>();
+      int flagMcMatchRecLb = std::abs(candidate.flagMcMatchRec());
 
-      if (std::abs(candidate.flagMcMatchRec()) == 1 << hf_cand_lb::DecayType::LbToLcPi) {
+      if (TESTBIT(flagMcMatchRecLb, hf_cand_lb::DecayType::LbToLcPi)) {
 
         auto indexMother = RecoDecay::getMother(mcParticles, candidate.prong1_as<TracksWExtMc>().mcParticle_as<soa::Join<aod::McParticles, aod::HfCandLbMcGen>>(), o2::constants::physics::Pdg::kLambdaB0, true);
         auto particleMother = mcParticles.rawIteratorAt(indexMother);

PWGHF/D2H/Tasks/taskLbReduced.cxx

@@ -379,7 +379,7 @@ struct HfTaskLbReduced {
   {
     auto ptCandLb = candidate.pt();
     auto invMassLb = hfHelper.invMassLbToLcPi(candidate);
-    auto candLc = candidate.template prong0Lc_as<CandsLc>();
+    auto candLc = candidate.template prong0_as<CandsLc>();
     auto ptLc = candidate.ptProng0();
     auto invMassLc = candLc.invMassHypo0() > 0 ? candLc.invMassHypo0() : candLc.invMassHypo1();
     // TODO: here we are assuming that only one of the two hypotheses is filled, to be checked
@@ -538,7 +538,7 @@ struct HfTaskLbReduced {
         if constexpr (withLbMl) {
           candidateMlScoreSig = candidate.mlProbLbToLcPi();
         }
-        auto prong1 = candidate.template prong1Track_as<TracksPion>();
+        auto prong1 = candidate.template prong1_as<TracksPion>();
 
         float ptMother = -1.;
         if constexpr (doMc) {

PWGHF/DataModel/CandidateReconstructionTables.h

@@ -1931,8 +1931,6 @@ DECLARE_SOA_TABLE(HfCandLbBase, "AOD", "HFCANDLBBASE",
                   hf_cand::PxProng1, hf_cand::PyProng1, hf_cand::PzProng1,
                   hf_cand::ImpactParameter0, hf_cand::ImpactParameter1,
                   hf_cand::ErrorImpactParameter0, hf_cand::ErrorImpactParameter1,
-                  hf_cand_lb::Prong0Id, hf_track_index::Prong1Id,
-                  hf_track_index::HFflag,
                   /* dynamic columns */
                   hf_cand_2prong::M<hf_cand::PxProng0, hf_cand::PyProng0, hf_cand::PzProng0, hf_cand::PxProng1, hf_cand::PyProng1, hf_cand::PzProng1>,
                   hf_cand_2prong::M2<hf_cand::PxProng0, hf_cand::PyProng0, hf_cand::PzProng0, hf_cand::PxProng1, hf_cand::PyProng1, hf_cand::PzProng1>,
@@ -1958,7 +1956,10 @@ DECLARE_SOA_TABLE(HfCandLbBase, "AOD", "HFCANDLBBASE",
 DECLARE_SOA_EXTENDED_TABLE_USER(HfCandLbExt, HfCandLbBase, "HFCANDLBEXT",
                                 hf_cand_2prong::Px, hf_cand_2prong::Py, hf_cand_2prong::Pz);
 
-using HfCandLb = HfCandLbExt;
+DECLARE_SOA_TABLE(HfCandLbProngs, "AOD", "HFCANDLBPRONGS",
+                  hf_cand_lb::Prong0Id, hf_track_index::Prong1Id);
+
+using HfCandLb = soa::Join<HfCandLbExt, HfCandLbProngs>;
 
 // table with results of reconstruction level MC matching
 DECLARE_SOA_TABLE(HfCandLbMcRec, "AOD", "HFCANDLBMCREC", //!

PWGHF/TableProducer/candidateCreatorLb.cxx

@@ -43,7 +43,7 @@ using namespace o2::hf_trkcandsel;
 /// Reconstruction of Λb candidates
 struct HfCandidateCreatorLb {
   Produces<aod::HfCandLbBase> rowCandidateBase;
-
+  Produces<aod::HfCandLbProngs> rowCandidateProngs;
   // vertexing
   Configurable<float> bz{"bz", 20., "magnetic field"};
   Configurable<bool> propagateToPCA{"propagateToPCA", true, "create tracks version propagated to PCA"};
@@ -217,8 +217,6 @@ struct HfCandidateCreatorLb {
         auto errorDecayLength = std::sqrt(getRotatedCovMatrixXX(covMatrixPV, phi, theta) + getRotatedCovMatrixXX(covMatrixPCA, phi, theta));
         auto errorDecayLengthXY = std::sqrt(getRotatedCovMatrixXX(covMatrixPV, phi, 0.) + getRotatedCovMatrixXX(covMatrixPCA, phi, 0.));
 
-        int hfFlag = 1 << hf_cand_lb::DecayType::LbToLcPi;
-
         // fill the candidate table for the Lb here:
         rowCandidateBase(collision.globalIndex(),
                          collision.posX(), collision.posY(), collision.posZ(),
@@ -228,10 +226,8 @@ struct HfCandidateCreatorLb {
                          pvecLc[0], pvecLc[1], pvecLc[2],
                          pvecPion[0], pvecPion[1], pvecPion[2],
                          impactParameter0.getY(), impactParameter1.getY(),
-                         std::sqrt(impactParameter0.getSigmaY2()), std::sqrt(impactParameter1.getSigmaY2()),
-                         lcCand.globalIndex(), trackPion.globalIndex(),
-                         hfFlag);
-
+                         std::sqrt(impactParameter0.getSigmaY2()), std::sqrt(impactParameter1.getSigmaY2()));
+        rowCandidateProngs(lcCand.globalIndex(), trackPion.globalIndex());
         // calculate invariant mass
         auto arrayMomenta = std::array{pvecLc, pvecPion};
         massLcPi = RecoDecay::m(std::move(arrayMomenta), std::array{massLc, massPi});
@@ -257,21 +253,19 @@ struct HfCandidateCreatorLbExpressions {
   /// @brief dummy process function, to be run on data
   void process(aod::Tracks const&) {}
 
-  void processMc(aod::HfCand3Prong const& lcCandidates,
-                 aod::TracksWMc const& tracks,
-                 aod::McParticles const& mcParticles)
+  void processMc(aod::HfCand3Prong const&,
+                 aod::TracksWMc const&,
+                 aod::McParticles const& mcParticles,
+                 aod::HfCandLbProngs const& candsLb)
   {
     int indexRec = -1;
     int8_t sign = 0;
     int8_t flag = 0;
     int8_t origin = 0;
     int8_t debug = 0;
 
-    rowCandidateLb->bindExternalIndices(&tracks);
-    rowCandidateLb->bindExternalIndices(&lcCandidates);
-
     // Match reconstructed candidates.
-    for (const auto& candidate : *rowCandidateLb) {
+    for (const auto& candidate : candsLb) {
       flag = 0;
       origin = 0;
       debug = 0;

PWGHF/TableProducer/candidateSelectorLbToLcPi.cxx

@@ -164,14 +164,6 @@ struct HfCandidateSelectorLbToLcPi {
     for (const auto& hfCandLb : hfCandLbs) { // looping over Lb candidates
 
       int statusLb = 0;
-
-      // check if flagged as Λb --> Λc+ π-
-      if (!(hfCandLb.hfflag() & 1 << hf_cand_lb::DecayType::LbToLcPi)) {
-        hfSelLbToLcPiCandidate(statusLb);
-        // LOGF(debug, "Lb candidate selection failed at hfflag check");
-        continue;
-      }
-
       // Lc is always index0 and pi is index1 by default
       // auto candLc = hfCandLb.prong0();
       auto candLc = hfCandLb.prong0_as<soa::Join<aod::HfCand3Prong, aod::HfSelLc>>();