[PWGHF] adding DstarTrack resonance (AliceO2Group#8014)

* [PWGHF] adding DstarTrack resonance

* [PWGHF] removing index and renaming candV0Tr

* [PWGHF] fixing indeces

* [PWGHF] fixing bug candV0Tr

* Please consider the following formatting changes

---------

Co-authored-by: ALICE Action Bot <[email protected]>

Author: cterrevo

PWGHF/D2H/DataModel/ReducedDataModel.h

@@ -756,6 +756,15 @@ DECLARE_SOA_TABLE(HfMcGenRedResos, "AOD", "HFMCGENREDRESO", //! Generation-level
                   hf_b0_mc::YProng1,
                   hf_b0_mc::EtaProng1,
                   o2::soa::Marker<1>);
+
+DECLARE_SOA_TABLE(HfCandChaResTr, "AOD", "HFCANDCHARESTR", //! Table with Resonance candidate information for resonances plus tracks reduced workflow
+                                                           // Static
+                  hf_cand::PxProng0, hf_cand::PyProng0, hf_cand::PzProng0,
+                  hf_cand::PxProng1, hf_cand::PyProng1, hf_cand::PzProng1,
+                  hf_reso_cand_reduced::InvMass,
+                  hf_reso_cand_reduced::InvMassProng0,
+                  // Dynamic
+                  hf_reso_cand_reduced::PtProng0<hf_cand::PxProng0, hf_cand::PyProng0>);
 } // namespace aod
 
 namespace soa

PWGHF/D2H/TableProducer/candidateCreatorCharmResoReduced.cxx

@@ -38,13 +38,15 @@ enum Selections : uint8_t {
   NoSel = 0,
   DSel,
   V0Sel,
+  TrackSel,
   NSelSteps
 };
 enum DecayChannel : uint8_t {
   Ds1ToDstarK0s = 0,
   Ds2StarToDplusK0s,
   XcToDplusLambda,
-  LambdaDminus
+  LambdaDminus,
+  DstarTrack
 };
 enum V0Type : uint8_t {
   K0s = 0,
@@ -67,6 +69,7 @@ auto vecBinsPt = std::vector<double>{binsPt, binsPt + nBinsPt + 1};
 struct HfCandidateCreatorCharmResoReduced {
   // Produces: Tables with resonance info
   Produces<aod::HfCandCharmReso> rowCandidateReso;
+  Produces<aod::HfCandChaResTr> rowCandidateResoTrack;
   // Optional daughter ML scores table
   Produces<aod::HfCharmResoMLs> mlScores;
 
@@ -89,11 +92,13 @@ struct HfCandidateCreatorCharmResoReduced {
   Partition<aod::HfRedVzeros> candidatesLambda = aod::hf_reso_v0::v0Type == (uint8_t)2 || aod::hf_reso_v0::v0Type == (uint8_t)4;
 
   Preslice<aod::HfRedVzeros> candsV0PerCollision = aod::hf_track_index_reduced::hfRedCollisionId;
+  Preslice<aod::HfRedTracks> candsTrackPerCollision = aod::hf_track_index_reduced::hfRedCollisionId;
   Preslice<aod::HfRed3PrNoTrks> candsDPerCollision = hf_track_index_reduced::hfRedCollisionId;
 
   // Useful constants
   double massK0{0.};
   double massLambda{0.};
+  double massProton{0.};
   double massDplus{0.};
   double massDstar{0.};
   double massD0{0.};
@@ -103,7 +108,7 @@ struct HfCandidateCreatorCharmResoReduced {
   void init(InitContext const&)
   {
     // check that only one process function is enabled
-    std::array<bool, 8> doprocess{doprocessDs2StarToDplusK0s, doprocessDs2StarToDplusK0sWithMl, doprocessDs1ToDstarK0s, doprocessDs1ToDstarK0sWithMl, doprocessXcToDplusLambda, doprocessXcToDplusLambdaWithMl, doprocessLambdaDminus, doprocessLambdaDminusWithMl};
+    std::array<bool, 10> doprocess{doprocessDs2StarToDplusK0s, doprocessDs2StarToDplusK0sWithMl, doprocessDs1ToDstarK0s, doprocessDs1ToDstarK0sWithMl, doprocessXcToDplusLambda, doprocessXcToDplusLambdaWithMl, doprocessLambdaDminus, doprocessLambdaDminusWithMl, doprocessDstarTrack, doprocessDstarTrackWithMl};
     if ((std::accumulate(doprocess.begin(), doprocess.end(), 0)) != 1) {
       LOGP(fatal, "Only one process function should be enabled! Please check your configuration!");
     }
@@ -113,12 +118,15 @@ struct HfCandidateCreatorCharmResoReduced {
     registry.add("hMassDs2Star", "Ds^{*}2 candidates; m_Ds^{*}2 (GeV/#it{c}^{2}) ;entries", {HistType::kTH2F, {{100, 2.4, 2.7}, {(std::vector<double>)binsPt, "#it{p}_{T} (GeV/#it{c})"}}});
     registry.add("hMassXcRes", "XcRes candidates; m_XcRes (GeV/#it{c}^{2}) ;entries", {HistType::kTH2F, {{100, 2.9, 3.3}, {(std::vector<double>)binsPt, "#it{p}_{T} (GeV/#it{c})"}}});
     registry.add("hMassLambdaDminus", "LambdaDminus candidates; m_LambdaDminus (GeV/#it{c}^{2}) ;entries", {HistType::kTH2F, {{100, 2.9, 3.3}, {(std::vector<double>)binsPt, "#it{p}_{T} (GeV/#it{c})"}}});
+    registry.add("hMassDstarTrack", "DstarTrack candidates; m_DstarTrack (GeV/#it{c}^{2}) ;entries", {HistType::kTH2F, {{100, 0.9, 1.4}, {(std::vector<double>)binsPt, "#it{p}_{T} (GeV/#it{c})"}}});
+
     if (activateQA) {
       constexpr int kNBinsSelections = Selections::NSelSteps;
       std::string labels[kNBinsSelections];
       labels[Selections::NoSel] = "No selection";
       labels[Selections::DSel] = "D Candidates Selection";
       labels[Selections::V0Sel] = "D & V0 candidate Selection";
+      labels[Selections::TrackSel] = "D & Track candidate Selection";
       static const AxisSpec axisSelections = {kNBinsSelections, 0.5, kNBinsSelections + 0.5, ""};
       registry.add("hSelections", "Selections", {HistType::kTH1F, {axisSelections}});
       for (int iBin = 0; iBin < kNBinsSelections; ++iBin) {
@@ -128,6 +136,7 @@ struct HfCandidateCreatorCharmResoReduced {
     // mass constants
     massK0 = o2::constants::physics::MassK0Short;
     massLambda = o2::constants::physics::MassLambda;
+    massProton = o2::constants::physics::MassProton;
     massDplus = o2::constants::physics::MassDPlus;
     massDstar = o2::constants::physics::MassDStar;
     massD0 = o2::constants::physics::MassD0;
@@ -147,7 +156,7 @@ struct HfCandidateCreatorCharmResoReduced {
     // slection on D candidate mass
     if (channel == DecayChannel::Ds2StarToDplusK0s || channel == DecayChannel::XcToDplusLambda || channel == DecayChannel::LambdaDminus) {
       invMassD = candD.invMassDplus();
-    } else if (channel == DecayChannel::Ds1ToDstarK0s) {
+    } else if (channel == DecayChannel::Ds1ToDstarK0s || channel == DecayChannel::DstarTrack) {
       if (candD.dType() > 0)
         invMassD = candD.invMassDstar();
       else
@@ -169,7 +178,7 @@ struct HfCandidateCreatorCharmResoReduced {
     return true;
   }
 
-  /// Basic selection of V0 candidates
+  /// Basic selection of V0 and track candidates
   /// \param candV0 is the reduced V0 candidate
   /// \param candD is the reduced D meson candidate
   /// \return true if selections are passed
@@ -215,10 +224,10 @@ struct HfCandidateCreatorCharmResoReduced {
     return true;
   }
 
-  template <bool fillMl, DecayChannel channel, typename Coll, typename DRedTable, typename V0RedTable>
+  template <bool fillMl, DecayChannel channel, typename Coll, typename DRedTable, typename V0TrRedTable>
   void runCandidateCreation(Coll const& collision,
                             DRedTable const& candsD,
-                            V0RedTable const& candsV0)
+                            V0TrRedTable const& candsV0Tr)
   {
     // loop on D candidates
     for (const auto& candD : candsD) {
@@ -242,71 +251,92 @@ struct HfCandidateCreatorCharmResoReduced {
       std::array<float, 3> pVecD = {candD.px(), candD.py(), candD.pz()};
       std::array<int, 3> dDaughtersIds = {candD.prong0Id(), candD.prong1Id(), candD.prong2Id()};
 
-      // loop on V0 candidates
+      // loop on V0 or track candidates
       bool alreadyCounted{false};
-      for (const auto& candV0 : candsV0) {
+      for (const auto& candV0Tr : candsV0Tr) {
         if (rejectDV0PairsWithCommonDaughter) {
           const std::array<int, 3> dDaughtersIDs = {candD.prong0Id(), candD.prong1Id(), candD.prong2Id()};
-          if (std::find(dDaughtersIDs.begin(), dDaughtersIDs.end(), candV0.prong0Id()) != dDaughtersIDs.end() || std::find(dDaughtersIDs.begin(), dDaughtersIDs.end(), candV0.prong1Id()) != dDaughtersIDs.end()) {
+          if constexpr (channel == DecayChannel::DstarTrack) {
+            if (std::find(dDaughtersIDs.begin(), dDaughtersIDs.end(), candV0Tr.globalIndex()) != dDaughtersIDs.end()) {
+              continue;
+            }
+          } else {
+            if (std::find(dDaughtersIDs.begin(), dDaughtersIDs.end(), candV0Tr.prong0Id()) != dDaughtersIDs.end() || std::find(dDaughtersIDs.begin(), dDaughtersIDs.end(), candV0Tr.prong1Id()) != dDaughtersIDs.end()) {
+              continue;
+            }
+          }
+        }
+        if constexpr (channel != DecayChannel::DstarTrack) {
+          if (!isV0Selected<channel>(candV0Tr, candD)) {
             continue;
           }
+          if (activateQA && !alreadyCounted) {
+            registry.fill(HIST("hSelections"), 1 + Selections::V0Sel);
+            alreadyCounted = true;
+          }
         }
 
-        if (!isV0Selected<channel>(candV0, candD)) {
-          continue;
-        }
-        if (activateQA && !alreadyCounted) {
-          registry.fill(HIST("hSelections"), 1 + Selections::V0Sel);
-          alreadyCounted = true;
-        }
         float invMassReso{0.};
         float invMassV0{0.};
-        std::array<float, 3> pVecV0 = {candV0.px(), candV0.py(), candV0.pz()};
-        float ptReso = RecoDecay::pt(RecoDecay::sumOfVec(pVecV0, pVecD));
-        switch (channel) {
-          case DecayChannel::Ds1ToDstarK0s:
-            invMassV0 = candV0.invMassK0s();
-            invMassReso = RecoDecay::m(std::array{pVecD, pVecV0}, std::array{massDstar, massK0});
-            registry.fill(HIST("hMassDs1"), invMassReso, ptReso);
-            break;
-          case DecayChannel::Ds2StarToDplusK0s:
-            invMassV0 = candV0.invMassK0s();
-            invMassReso = RecoDecay::m(std::array{pVecD, pVecV0}, std::array{massDplus, massK0});
-            registry.fill(HIST("hMassDs2Star"), invMassReso, ptReso);
-            break;
-          case DecayChannel::XcToDplusLambda:
-            if (candD.dType() > 0) {
-              invMassV0 = candV0.invMassLambda();
-            } else {
-              invMassV0 = candV0.invMassAntiLambda();
-            }
-            invMassReso = RecoDecay::m(std::array{pVecD, pVecV0}, std::array{massDplus, massLambda});
-            registry.fill(HIST("hMassXcRes"), invMassReso, ptReso);
-            break;
-          case DecayChannel::LambdaDminus:
-            if (candD.dType() < 0) {
-              invMassV0 = candV0.invMassLambda();
-            } else {
-              invMassV0 = candV0.invMassAntiLambda();
-            }
-            invMassReso = RecoDecay::m(std::array{pVecD, pVecV0}, std::array{massDplus, massLambda});
-            registry.fill(HIST("hMassLambdaDminus"), invMassReso, ptReso);
-            break;
-          default:
-            break;
+        std::array<float, 3> pVecV0Tr = {candV0Tr.px(), candV0Tr.py(), candV0Tr.pz()};
+        float ptReso = RecoDecay::pt(RecoDecay::sumOfVec(pVecV0Tr, pVecD));
+
+        if constexpr (channel == DecayChannel::DstarTrack) {
+          invMassReso = RecoDecay::m(std::array{pVecD, pVecV0Tr}, std::array{massDstar, massProton});
+          registry.fill(HIST("hMassDstarTrack"), invMassReso, ptReso);
+        } else {
+          switch (channel) {
+            case DecayChannel::Ds1ToDstarK0s:
+              invMassV0 = candV0Tr.invMassK0s();
+              invMassReso = RecoDecay::m(std::array{pVecD, pVecV0Tr}, std::array{massDstar, massK0});
+              registry.fill(HIST("hMassDs1"), invMassReso, ptReso);
+              break;
+            case DecayChannel::Ds2StarToDplusK0s:
+              invMassV0 = candV0Tr.invMassK0s();
+              invMassReso = RecoDecay::m(std::array{pVecD, pVecV0Tr}, std::array{massDplus, massK0});
+              registry.fill(HIST("hMassDs2Star"), invMassReso, ptReso);
+              break;
+            case DecayChannel::XcToDplusLambda:
+              if (candD.dType() > 0) {
+                invMassV0 = candV0Tr.invMassLambda();
+              } else {
+                invMassV0 = candV0Tr.invMassAntiLambda();
+              }
+              invMassReso = RecoDecay::m(std::array{pVecD, pVecV0Tr}, std::array{massDplus, massLambda});
+              registry.fill(HIST("hMassXcRes"), invMassReso, ptReso);
+              break;
+            case DecayChannel::LambdaDminus:
+              if (candD.dType() < 0) {
+                invMassV0 = candV0Tr.invMassLambda();
+              } else {
+                invMassV0 = candV0Tr.invMassAntiLambda();
+              }
+              invMassReso = RecoDecay::m(std::array{pVecD, pVecV0Tr}, std::array{massDplus, massLambda});
+              registry.fill(HIST("hMassLambdaDminus"), invMassReso, ptReso);
+              break;
+            default:
+              break;
+          }
         }
         // Filling Output table
-        rowCandidateReso(collision.globalIndex(),
-                         candD.globalIndex(),
-                         candV0.globalIndex(),
-                         pVecD[0], pVecD[1], pVecD[2],
-                         pVecV0[0], pVecV0[1], pVecV0[2],
-                         invMassReso,
-                         invMassD,
-                         invMassV0,
-                         candV0.cpa(),
-                         candV0.dca(),
-                         candV0.v0Radius());
+        if constexpr (channel == DecayChannel::DstarTrack) {
+          rowCandidateResoTrack(pVecD[0], pVecD[1], pVecD[2],
+                                candV0Tr.px(), candV0Tr.py(), candV0Tr.pz(),
+                                invMassReso,
+                                invMassD);
+        } else {
+          rowCandidateReso(collision.globalIndex(),
+                           candD.globalIndex(),
+                           candV0Tr.globalIndex(),
+                           pVecD[0], pVecD[1], pVecD[2],
+                           pVecV0Tr[0], pVecV0Tr[1], pVecV0Tr[2],
+                           invMassReso,
+                           invMassD,
+                           invMassV0,
+                           candV0Tr.cpa(),
+                           candV0Tr.dca(),
+                           candV0Tr.v0Radius());
+        }
         if constexpr (fillMl) {
           mlScores(candD.mlScoreBkgMassHypo0(), candD.mlScorePromptMassHypo0(), candD.mlScoreNonpromptMassHypo0());
         }
@@ -417,6 +447,32 @@ struct HfCandidateCreatorCharmResoReduced {
     }
   }
   PROCESS_SWITCH(HfCandidateCreatorCharmResoReduced, processLambdaDminusWithMl, "Process LambdaDminus candidates with Ml info", false);
+  void processDstarTrack(aod::HfRedCollisions const& collisions,
+                         aod::HfRed3PrNoTrks const& candsD,
+                         soa::Join<aod::HfRedTrackBases, aod::HfRedTracksCov> const& candidatesTrack)
+  {
+    for (const auto& collision : collisions) {
+      auto thisCollId = collision.globalIndex();
+      auto candsDThisColl = candsD.sliceBy(candsDPerCollision, thisCollId);
+      auto trackThisColl = candidatesTrack.sliceBy(candsTrackPerCollision, thisCollId);
+      runCandidateCreation<false, DecayChannel::DstarTrack>(collision, candsDThisColl, trackThisColl);
+    }
+  }
+  PROCESS_SWITCH(HfCandidateCreatorCharmResoReduced, processDstarTrack, "Process DStar candidates without Ml info", false);
+
+  void processDstarTrackWithMl(aod::HfRedCollisions const& collisions,
+                               soa::Join<aod::HfRed3PrNoTrks, aod::HfRed3ProngsMl> const& candsD,
+                               soa::Join<aod::HfRedTrackBases, aod::HfRedTracksCov> const& candidatesTrack)
+  //                                    aod::HfRedTracks const&)
+  {
+    for (const auto& collision : collisions) {
+      auto thisCollId = collision.globalIndex();
+      auto candsDThisColl = candsD.sliceBy(candsDPerCollision, thisCollId);
+      auto trackThisColl = candidatesTrack.sliceBy(candsTrackPerCollision, thisCollId);
+      runCandidateCreation<true, DecayChannel::DstarTrack>(collision, candsDThisColl, trackThisColl);
+    }
+  }
+  PROCESS_SWITCH(HfCandidateCreatorCharmResoReduced, processDstarTrackWithMl, "Process DStar candidates with Ml info", false);
 
 }; // struct