PWGHF: Developement of D Resonances MC analysis

PWGHF/D2H/Core/SelectorCutsRedDataFormat.h

@@ -22,9 +22,11 @@
 #include <string> // std::string
 #include <vector> // std::vector
 
+namespace o2::analysis
+{
 namespace hf_cuts_d_daughter
 {
-const int nBinsPt = 7;
+static constexpr int nBinsPt = 7;
 static constexpr int nCutVars = 6;
 constexpr double binsPt[nBinsPt + 1] = {
   1.,
@@ -45,15 +47,22 @@ constexpr double cuts[nBinsPt][nCutVars] = {{1.84, 1.89, 1.77, 1.81, 1.92, 1.96}
                                             {1.84, 1.89, 1.77, 1.81, 1.92, 1.96},  /* 12   < pt < 24 */
                                             {1.84, 1.89, 1.77, 1.81, 1.92, 1.96}}; /* 24   < pt < 1000 */
 // row labels
-static const std::vector<std::string> labelsPt{};
+static const std::vector<std::string> labelsPt{
+  "pT bin 0",
+  "pT bin 1",
+  "pT bin 2",
+  "pT bin 3",
+  "pT bin 4",
+  "pT bin 5",
+  "pT bin 6"};
 // column labels
 static const std::vector<std::string> labelsCutVar = {"invMassSignalLow", "invMassSignalHigh", "invMassLeftSBLow", "invMassLeftSBHigh", "invMassRightSBLow", "invMassRightSBHigh"};
 } // namespace hf_cuts_d_daughter
 
 // namespace with v0 selections for reduced charmed-resonances analysis
 namespace hf_cuts_v0_daughter
 {
-const int nBinsPt = 7;
+static constexpr int nBinsPt = 7;
 static constexpr int nCutVars = 5;
 constexpr double binsPt[nBinsPt + 1] = {
   0.,
@@ -78,4 +87,5 @@ static const std::vector<std::string> labelsPt{};
 // column labels
 static const std::vector<std::string> labelsCutVar = {"invMassLow", "invMassHigh", "cpaMin", "dcaMax", "radiusMin"};
 } // namespace hf_cuts_v0_daughter
+} // namespace o2::analysis
 #endif // PWGHF_D2H_CORE_SELECTORCUTSREDDATAFORMAT_H_

PWGHF/D2H/DataModel/ReducedDataModel.h

@@ -687,7 +687,8 @@ DECLARE_SOA_INDEX_COLUMN_FULL(Prong0, prong0, int, HfRed3PrNoTrks, "_0"); //! Pr
 DECLARE_SOA_INDEX_COLUMN_FULL(Prong1, prong1, int, HfRedVzeros, "_1");    //! Prong1 index (V0 daughter)
 DECLARE_SOA_COLUMN(FlagMcMatchRec, flagMcMatchRec, int8_t);               // flag for decay channel classification reconstruction level
 DECLARE_SOA_COLUMN(FlagMcMatchGen, flagMcMatchGen, int8_t);               // flag for decay channel classification generator level
-DECLARE_SOA_COLUMN(DebugMcRec, debugMcRec, int8_t);                       // debug flag for mis-association reconstruction level
+DECLARE_SOA_COLUMN(DebugMcRec, debugMcRec, int8_t);                       // debug flag for mis-association at reconstruction level
+DECLARE_SOA_COLUMN(Origin, origin, int8_t);                               // Flag for origin of MC particle 1=promt, 2=FD
 
 DECLARE_SOA_DYNAMIC_COLUMN(Pt, pt, //!
                            [](float pxProng0, float pxProng1, float pyProng0, float pyProng1) -> float { return RecoDecay::pt((1.f * pxProng0 + 1.f * pxProng1), (1.f * pyProng0 + 1.f * pyProng1)); });
@@ -741,11 +742,13 @@ DECLARE_SOA_TABLE(HfMcRecRedDV0s, "AOD", "HFMCRECREDDV0", //! Table with reconst
                   hf_reso_cand_reduced::Prong1Id,
                   hf_reso_cand_reduced::FlagMcMatchRec,
                   hf_reso_cand_reduced::DebugMcRec,
+                  hf_reso_cand_reduced::Origin,
                   hf_b0_mc::PtMother,
                   o2::soa::Marker<1>);
 
-DECLARE_SOA_TABLE(HfMcGenRedResos, "AOD", "HFMCGENREDRESO", //! Generation-level MC information on B0 candidates for reduced workflow
+DECLARE_SOA_TABLE(HfMcGenRedResos, "AOD", "HFMCGENREDRESO", //! Generation-level MC information on Ds-Resonances candidates for reduced workflow
                   hf_cand_b0::FlagMcMatchGen,
+                  hf_reso_cand_reduced::Origin,
                   hf_b0_mc::PtTrack,
                   hf_b0_mc::YTrack,
                   hf_b0_mc::EtaTrack,
@@ -765,6 +768,14 @@ DECLARE_SOA_TABLE(HfCandChaResTr, "AOD", "HFCANDCHARESTR", //! Table with Resona
                   hf_reso_cand_reduced::InvMassProng0,
                   // Dynamic
                   hf_reso_cand_reduced::PtProng0<hf_cand::PxProng0, hf_cand::PyProng0>);
+
+// Table with same size as HfCandCharmReso
+DECLARE_SOA_TABLE(HfMcRecRedResos, "AOD", "HFMCRECREDRESO", //! Reconstruction-level MC information on Ds-Resonances candidates for reduced workflow
+                  hf_reso_cand_reduced::FlagMcMatchRec,
+                  hf_reso_cand_reduced::DebugMcRec,
+                  hf_reso_cand_reduced::Origin,
+                  hf_b0_mc::PtMother,
+                  o2::soa::Marker<1>);
 } // namespace aod
 
 namespace soa

PWGHF/D2H/TableProducer/candidateCreatorCharmResoReduced.cxx

@@ -48,12 +48,20 @@ enum DecayChannel : uint8_t {
   LambdaDminus,
   DstarTrack
 };
+
 enum V0Type : uint8_t {
   K0s = 0,
   Lambda,
   AntiLambda
 };
 
+enum DecayTypeMc : uint8_t {
+  Ds1ToDStarK0ToD0PiK0s = 0,
+  Ds2StarToDplusK0sToPiKaPiPiPi,
+  Ds1ToDStarK0ToDPlusPi0K0s,
+  Ds1ToDStarK0ToDPlusGammaK0s
+};
+
 const int nBinsPt = 7;
 constexpr double binsPt[nBinsPt + 1] = {
   1.,
@@ -83,6 +91,7 @@ struct HfCandidateCreatorCharmResoReduced {
   Configurable<LabeledArray<double>> cutsD{"cutsDdaughter", {hf_cuts_d_daughter::cuts[0], hf_cuts_d_daughter::nBinsPt, hf_cuts_d_daughter::nCutVars, hf_cuts_d_daughter::labelsPt, hf_cuts_d_daughter::labelsCutVar}, "D daughter selections"};
   Configurable<std::vector<double>> binsPtD{"binsPtD", std::vector<double>{hf_cuts_d_daughter::vecBinsPt}, "pT bin limits for D daughter cuts"};
   Configurable<LabeledArray<double>> cutsV0{"cutsV0daughter", {hf_cuts_v0_daughter::cuts[0], hf_cuts_v0_daughter::nBinsPt, hf_cuts_v0_daughter::nCutVars, hf_cuts_v0_daughter::labelsPt, hf_cuts_v0_daughter::labelsCutVar}, "V0 daughter selections"};
+
   Configurable<std::vector<double>> binsPtV0{"binsPtV0", std::vector<double>{hf_cuts_v0_daughter::vecBinsPt}, "pT bin limits for V0 daughter cuts"};
 
   using reducedDWithMl = soa::Join<aod::HfRed3PrNoTrks, aod::HfRed3ProngsMl>;
@@ -114,7 +123,7 @@ struct HfCandidateCreatorCharmResoReduced {
     }
     // histograms
     const AxisSpec axisPt{(std::vector<double>)vecBinsPt, "#it{p}_{T} (GeV/#it{c})"};
-    registry.add("hMassDs1", "Ds1 candidates;m_{Ds1} (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{100, 2.4, 2.7}, {(std::vector<double>)binsPt, "#it{p}_{T} (GeV/#it{c})"}}});
+    registry.add("hMassDs1", "Ds1 candidates;m_{Ds1} (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{200, 2.5, 2.7}, {(std::vector<double>)binsPt, "#it{p}_{T} (GeV/#it{c})"}}});
     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})"}}});
@@ -474,9 +483,62 @@ struct HfCandidateCreatorCharmResoReduced {
   }
   PROCESS_SWITCH(HfCandidateCreatorCharmResoReduced, processDstarTrackWithMl, "Process DStar candidates with Ml info", false);
 
-}; // struct
+}; // struct HfCandidateCreatorCharmResoReduced
+
+struct HfCandidateCreatorCharmResoReducedExpressions {
+
+  Produces<aod::HfMcRecRedResos> rowResoMcRec;
+
+  HistogramRegistry registry{"registry"};
+
+  void init(InitContext const&)
+  {
+    const AxisSpec axisPt{(std::vector<double>)vecBinsPt, "#it{p}_{T} (GeV/#it{c})"};
+    registry.add("hMassMcMatched", "Reso MC candidates Matched with generate particle;m (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{200, 2.5, 2.7}, {axisPt}}});
+    registry.add("hMassMcMatchedIncomplete", "Reso MC candidates Matched with generate particle w. Invcomplete decay;m (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{200, 2.5, 2.7}, {axisPt}}});
+    registry.add("hMassMcUnmatched", "Reso MC candidates NOT Matched with generate particle;m (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{200, 2.5, 2.7}, {axisPt}}});
+    registry.add("hMassMcNoEntry", "Reso MC candidates w.o. entry in MC Reco table;m (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{200, 2.5, 2.7}, {axisPt}}});
+  }
+
+  /// Fill candidate information at MC reconstruction level
+  /// \param rowsDV0McRec MC reco information on DPi pairs
+  /// \param candsReso prong global indices of B0 candidates
+  template <typename McRec>
+  void fillResoMcRec(McRec const& rowsDV0McRec, aod::HfCandCharmReso const& candsReso)
+  {
+    for (const auto& candReso : candsReso) {
+      bool filledMcInfo{false};
+      for (const auto& rowDV0McRec : rowsDV0McRec) {
+        if ((rowDV0McRec.prong0Id() != candReso.prong0Id()) || (rowDV0McRec.prong1Id() != candReso.prong1Id())) {
+          continue;
+        }
+        rowResoMcRec(rowDV0McRec.flagMcMatchRec(), rowDV0McRec.debugMcRec(), rowDV0McRec.origin(), rowDV0McRec.ptMother());
+        filledMcInfo = true;
+        if (TESTBIT(rowDV0McRec.flagMcMatchRec(), DecayTypeMc::Ds1ToDStarK0ToD0PiK0s) || TESTBIT(rowDV0McRec.flagMcMatchRec(), DecayTypeMc::Ds2StarToDplusK0sToPiKaPiPiPi)) {
+          registry.fill(HIST("hMassMcMatched"), candReso.invMass(), candReso.pt());
+        } else if (TESTBIT(rowDV0McRec.flagMcMatchRec(), DecayTypeMc::Ds1ToDStarK0ToDPlusGammaK0s) || TESTBIT(rowDV0McRec.flagMcMatchRec(), DecayTypeMc::Ds1ToDStarK0ToDPlusPi0K0s)) {
+          registry.fill(HIST("hMassMcMatchedIncomplete"), candReso.invMass(), candReso.pt());
+        } else {
+          registry.fill(HIST("hMassMcUnmatched"), candReso.invMass(), candReso.pt());
+        }
 
+        break;
+      }
+      if (!filledMcInfo) { // protection to get same size tables in case something went wrong: we created a candidate that was not preselected in the D-Pi creator
+        rowResoMcRec(0, -1, -1, -1.f);
+        registry.fill(HIST("hMassMcNoEntry"), candReso.invMass(), candReso.pt());
+      }
+    }
+  }
+
+  void processMc(aod::HfMcRecRedDV0s const& rowsDV0McRec, aod::HfCandCharmReso const& candsReso)
+  {
+    fillResoMcRec(rowsDV0McRec, candsReso);
+  }
+  PROCESS_SWITCH(HfCandidateCreatorCharmResoReducedExpressions, processMc, "Process MC", false);
+};
 WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
 {
-  return WorkflowSpec{adaptAnalysisTask<HfCandidateCreatorCharmResoReduced>(cfgc)};
+  return WorkflowSpec{adaptAnalysisTask<HfCandidateCreatorCharmResoReduced>(cfgc),
+                      adaptAnalysisTask<HfCandidateCreatorCharmResoReducedExpressions>(cfgc)};
 }