diff --git a/PWGHF/D2H/Core/SelectorCutsRedDataFormat.h b/PWGHF/D2H/Core/SelectorCutsRedDataFormat.h index a698b8152a6..daa642f1601 100644 --- a/PWGHF/D2H/Core/SelectorCutsRedDataFormat.h +++ b/PWGHF/D2H/Core/SelectorCutsRedDataFormat.h @@ -22,9 +22,11 @@ #include // std::string #include // 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,7 +47,14 @@ 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 labelsPt{}; +static const std::vector 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 labelsCutVar = {"invMassSignalLow", "invMassSignalHigh", "invMassLeftSBLow", "invMassLeftSBHigh", "invMassRightSBLow", "invMassRightSBHigh"}; } // namespace hf_cuts_d_daughter @@ -53,7 +62,7 @@ static const std::vector labelsCutVar = {"invMassSignalLow", "invMa // 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 labelsPt{}; // column labels static const std::vector labelsCutVar = {"invMassLow", "invMassHigh", "cpaMin", "dcaMax", "radiusMin"}; } // namespace hf_cuts_v0_daughter +} // namespace o2::analysis #endif // PWGHF_D2H_CORE_SELECTORCUTSREDDATAFORMAT_H_ diff --git a/PWGHF/D2H/DataModel/ReducedDataModel.h b/PWGHF/D2H/DataModel/ReducedDataModel.h index fc84b17a70f..906859de48e 100644 --- a/PWGHF/D2H/DataModel/ReducedDataModel.h +++ b/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); + +// 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 diff --git a/PWGHF/D2H/TableProducer/candidateCreatorCharmResoReduced.cxx b/PWGHF/D2H/TableProducer/candidateCreatorCharmResoReduced.cxx index 935153adb83..0e39a651880 100644 --- a/PWGHF/D2H/TableProducer/candidateCreatorCharmResoReduced.cxx +++ b/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> 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> binsPtD{"binsPtD", std::vector{hf_cuts_d_daughter::vecBinsPt}, "pT bin limits for D daughter cuts"}; Configurable> 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> binsPtV0{"binsPtV0", std::vector{hf_cuts_v0_daughter::vecBinsPt}, "pT bin limits for V0 daughter cuts"}; using reducedDWithMl = soa::Join; @@ -114,7 +123,7 @@ struct HfCandidateCreatorCharmResoReduced { } // histograms const AxisSpec axisPt{(std::vector)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)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)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)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)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)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 rowResoMcRec; + + HistogramRegistry registry{"registry"}; + + void init(InitContext const&) + { + const AxisSpec axisPt{(std::vector)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 + 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(cfgc)}; + return WorkflowSpec{adaptAnalysisTask(cfgc), + adaptAnalysisTask(cfgc)}; } diff --git a/PWGHF/D2H/TableProducer/dataCreatorCharmResoReduced.cxx b/PWGHF/D2H/TableProducer/dataCreatorCharmResoReduced.cxx index 5e57bbd8598..f1cce927da1 100644 --- a/PWGHF/D2H/TableProducer/dataCreatorCharmResoReduced.cxx +++ b/PWGHF/D2H/TableProducer/dataCreatorCharmResoReduced.cxx @@ -77,7 +77,14 @@ enum DecayTypeMc : uint8_t { Ds1ToDStarK0ToD0PiK0s = 0, Ds2StarToDplusK0s, Ds1ToDStarK0ToDPlusPi0K0s, - Ds1ToDStarK0ToDPlusGammaK0s + Ds1ToDStarK0ToD0PiK0sPart +}; + +enum PartialMatchMc : uint8_t { + K0Matched = 0, + D0Matched, + DStarMatched, + DPlusMatched }; /// Creation of D-V0 pairs @@ -231,12 +238,13 @@ struct HfDataCreatorCharmResoReduced { registry.add("hMassDstarProton", "D^{*}-proton candidates;m_{D^{*}p} - m_{D^{*}} (GeV/#it{c}^{2});entries", {HistType::kTH1F, {{500, 0.9, 1.4}}}); registry.add("hDType", "D selection flag", {HistType::kTH1F, {{5, -2.5, 2.5}}}); - registry.add("hMCRecCounter", "Number of Reconstructed MC Matched candidates per channel", {HistType::kTH1F, {{11, -5.5, 5.5}}}); - registry.add("hMCRecDebug", "Debug of MC Reco", {HistType::kTH1F, {{5, -0.5, 4.5}}}); + registry.add("hMCRecCounter", "Number of Reconstructed MC Matched candidates per channel", {HistType::kTH1F, {{17, -8.5, 8.5}}}); + registry.add("hMCRecDebug", "Debug of MC Reco", {HistType::kTH1F, {{16, -0.5, 15.5}}}); + registry.add("hMCRecOrigin", "Origin of Matched particles", {HistType::kTH1F, {{3, -0.5, 2.5}}}); - registry.add("hMCGenCounter", "Number of Generated particles; Decay Channel Flag; pT [GeV/c]", {HistType::kTH2F, {{11, -5.5, 5.5}, {100, 0, 50}}}); + registry.add("hMCGenCounter", "Number of Generated particles; Decay Channel Flag; pT [GeV/c]", {HistType::kTH2F, {{17, -8.5, 8.5}, {100, 0, 50}}}); registry.add("hMCSignCounter", "Sign of Generated particles", {HistType::kTH1F, {{3, -1.5, 1.5}}}); - registry.add("hMCOriginCounter", "Origin of Generated particles", {HistType::kTH1F, {{3, -0.5, 2.5}}}); + registry.add("hMCGenOrigin", "Origin of Generated particles", {HistType::kTH1F, {{3, -0.5, 2.5}}}); registry.add("hMCOriginCounterWrongDecay", "Origin of Generated particles in Wrong decay", {HistType::kTH1F, {{3, -0.5, 2.5}}}); // Configure CCDB access @@ -501,85 +509,88 @@ struct HfDataCreatorCharmResoReduced { int8_t signV0{0}; int8_t flag{0}; int8_t debug{0}; + int8_t origin{0}; + std::vector idxBhadMothers{}; float motherPt{-1.f}; - bool fillHisto{false}; if constexpr (decChannel == DecayChannel::DstarV0) { // Ds1 → D* K0 → (D0 π+) K0s → ((K-π+) π+)(π+π-) - auto indexRec = RecoDecay::getMatchedMCRec(particlesMc, std::array{vecDaughtersReso[0], vecDaughtersReso[1], vecDaughtersReso[2], vecDaughtersReso[3], vecDaughtersReso[4]}, Pdg::kDS1, std::array{+kPiPlus, -kKPlus, +kPiPlus, +kPiPlus, -kPiPlus}, true, &sign, 3); - if (indexRec > -1) { + auto indexRecReso = RecoDecay::getMatchedMCRec(particlesMc, std::array{vecDaughtersReso[0], vecDaughtersReso[1], vecDaughtersReso[2], vecDaughtersReso[3], vecDaughtersReso[4]}, Pdg::kDS1, std::array{+kPiPlus, -kKPlus, +kPiPlus, +kPiPlus, -kPiPlus}, true, &sign, 3); + if (indexRecReso > -1) { // D* → π- K+ π- - indexRec = RecoDecay::getMatchedMCRec(particlesMc, std::array{vecDaughtersReso[0], vecDaughtersReso[1], vecDaughtersReso[2]}, Pdg::kDStar, std::array{-kKPlus, +kPiPlus, +kPiPlus}, true, &signDStar, 2); + auto indexRecDstar = RecoDecay::getMatchedMCRec(particlesMc, std::array{vecDaughtersReso[0], vecDaughtersReso[1], vecDaughtersReso[2]}, Pdg::kDStar, std::array{-kKPlus, +kPiPlus, +kPiPlus}, true, &signDStar, 2); auto indexRecD0 = RecoDecay::getMatchedMCRec(particlesMc, std::array{vecDaughtersReso[0], vecDaughtersReso[1]}, Pdg::kD0, std::array{+kPiPlus, -kKPlus}, true, &signD0, 1); auto indexRecK0 = RecoDecay::getMatchedMCRec(particlesMc, std::array{vecDaughtersReso[3], vecDaughtersReso[4]}, kK0, std::array{+kPiPlus, -kPiPlus}, true, &signV0, 2); - if (indexRec > -1 && indexRecD0 > -1 && indexRecK0 > -1) { + if (indexRecDstar > -1 && indexRecD0 > -1 && indexRecK0 > -1) { flag = sign * BIT(DecayTypeMc::Ds1ToDStarK0ToD0PiK0s); - } else { - if (indexRec <= -1) { - debug = 1; - LOGF(debug, "DS1 decays in the expected final state but Dstar does not!"); - } - if (indexRec > -1 && indexRecD0 <= -1) { - debug = 2; - LOGF(debug, "DS1 and Dstar decays in the expected final state but D0 does not!"); - } - if (indexRec > -1 && indexRecK0 <= -1) { - debug = 3; - LOGF(debug, "DS1 and Dstar decays in the expected final state but K0 does not!"); - } } - auto indexMother = RecoDecay::getMother(particlesMc, vecDaughtersReso.back().template mcParticle_as(), Pdg::kDS1, true); - if (indexMother >= 0) { - auto particleMother = particlesMc.rawIteratorAt(indexMother); - motherPt = particleMother.pt(); + } else if (RecoDecay::getMatchedMCRec(particlesMc, std::array{vecDaughtersReso[0], vecDaughtersReso[1], vecDaughtersReso[2], vecDaughtersReso[3], vecDaughtersReso[4]}, Pdg::kDS1, std::array{+kPiPlus, -kKPlus, +kPiPlus, +kPiPlus, -kPiPlus}, true, &sign, 3) > -1) { // Verify partly reconstructed decay Ds1 -> D* K0s -> D+ π0 K0s + auto indexRecDstar = RecoDecay::getMatchedMCRec(particlesMc, std::array{vecDaughtersReso[0], vecDaughtersReso[1], vecDaughtersReso[2]}, Pdg::kDStar, std::array{-kKPlus, +kPiPlus, +kPiPlus}, true, &signDStar, 2); + auto indexRecK0 = RecoDecay::getMatchedMCRec(particlesMc, std::array{vecDaughtersReso[3], vecDaughtersReso[4]}, kK0, std::array{+kPiPlus, -kPiPlus}, true, &signV0, 2); + if (indexRecDstar > -1 && indexRecK0 > -1) { + if (RecoDecay::getMatchedMCRec(particlesMc, std::array{vecDaughtersReso[0], vecDaughtersReso[1], vecDaughtersReso[2]}, Pdg::kDPlus, std::array{+kPiPlus, -kKPlus, +kPiPlus}, true, &signDPlus, 2) > -1) { + flag = sign * BIT(DecayTypeMc::Ds1ToDStarK0ToDPlusPi0K0s); + } else if (RecoDecay::getMatchedMCRec(particlesMc, std::array{vecDaughtersReso[0], vecDaughtersReso[1]}, Pdg::kD0, std::array{+kPiPlus, -kKPlus}, true, &signD0, 2) > -1) { + flag = sign * BIT(DecayTypeMc::Ds1ToDStarK0ToD0PiK0sPart); + } } - fillHisto = true; + } + if (flag != 0) { + auto particleReso = particlesMc.iteratorAt(indexRecReso); + origin = RecoDecay::getCharmHadronOrigin(particlesMc, particleReso, false, &idxBhadMothers); } } else if constexpr (decChannel == DecayChannel::DplusV0) { // Ds2Star → D+ K0 → (π+K-π+) K0s → (π+K-π+)(π+π-) - auto indexRec = RecoDecay::getMatchedMCRec(particlesMc, std::array{vecDaughtersReso[0], vecDaughtersReso[1], vecDaughtersReso[2], vecDaughtersReso[3], vecDaughtersReso[4]}, Pdg::kDS2Star, std::array{+kPiPlus, -kKPlus, +kPiPlus, +kPiPlus, -kPiPlus}, true, &sign, 3); - if (indexRec > -1) { + auto indexRecReso = RecoDecay::getMatchedMCRec(particlesMc, std::array{vecDaughtersReso[0], vecDaughtersReso[1], vecDaughtersReso[2], vecDaughtersReso[3], vecDaughtersReso[4]}, Pdg::kDS2Star, std::array{+kPiPlus, -kKPlus, +kPiPlus, +kPiPlus, -kPiPlus}, true, &sign, 3); + if (indexRecReso > -1) { // D+ → π+ K- π+ - indexRec = RecoDecay::getMatchedMCRec(particlesMc, std::array{vecDaughtersReso[0], vecDaughtersReso[1], vecDaughtersReso[2]}, Pdg::kDPlus, std::array{+kPiPlus, -kKPlus, +kPiPlus}, true, &signDPlus, 2); + auto indexRecDplus = RecoDecay::getMatchedMCRec(particlesMc, std::array{vecDaughtersReso[0], vecDaughtersReso[1], vecDaughtersReso[2]}, Pdg::kDPlus, std::array{+kPiPlus, -kKPlus, +kPiPlus}, true, &signDPlus, 2); auto indexRecK0 = RecoDecay::getMatchedMCRec(particlesMc, std::array{vecDaughtersReso[3], vecDaughtersReso[4]}, kK0, std::array{+kPiPlus, -kPiPlus}, true, &signV0, 2); - if (indexRec > -1 && indexRecK0 > -1) { + if (indexRecDplus > -1 && indexRecK0 > -1) { flag = sign * BIT(DecayTypeMc::Ds2StarToDplusK0s); - } else { - if (indexRec <= -1) { - debug = 1; - LOGF(debug, "DS2Star decays in the expected final state but DPlus does not!"); - } - if (indexRec > -1 && indexRecK0 <= -1) { - debug = 3; - LOGF(debug, "DS2Star and DPlus decays in the expected final state but K0 does not!"); - } } auto indexMother = RecoDecay::getMother(particlesMc, vecDaughtersReso.back().template mcParticle_as(), Pdg::kDS2Star, true); if (indexMother >= 0) { auto particleMother = particlesMc.rawIteratorAt(indexMother); motherPt = particleMother.pt(); } - fillHisto = true; } else { // Verify partly reconstructed decay Ds1 -> D* K0s -> D+ π0 K0s - indexRec = RecoDecay::getMatchedMCRec(particlesMc, std::array{vecDaughtersReso[0], vecDaughtersReso[1], vecDaughtersReso[2], vecDaughtersReso[3], vecDaughtersReso[4]}, Pdg::kDS1, std::array{+kPiPlus, -kKPlus, +kPiPlus, +kPiPlus, -kPiPlus}, true, &sign, 3); - if (indexRec > -1) { + indexRecReso = RecoDecay::getMatchedMCRec(particlesMc, std::array{vecDaughtersReso[0], vecDaughtersReso[1], vecDaughtersReso[2], vecDaughtersReso[3], vecDaughtersReso[4]}, Pdg::kDS1, std::array{+kPiPlus, -kKPlus, +kPiPlus, +kPiPlus, -kPiPlus}, true, &sign, 3); + if (indexRecReso > -1) { auto indexRecDstar = RecoDecay::getMatchedMCRec(particlesMc, std::array{vecDaughtersReso[0], vecDaughtersReso[1], vecDaughtersReso[2]}, Pdg::kDStar, std::array{-kKPlus, +kPiPlus, +kPiPlus}, true, &signDStar, 2); auto indexRecK0 = RecoDecay::getMatchedMCRec(particlesMc, std::array{vecDaughtersReso[3], vecDaughtersReso[4]}, kK0, std::array{+kPiPlus, -kPiPlus}, true, &signV0, 2); if (indexRecDstar > -1 && indexRecK0 > -1) { auto indexRecDplus = RecoDecay::getMatchedMCRec(particlesMc, std::array{vecDaughtersReso[0], vecDaughtersReso[1], vecDaughtersReso[2]}, Pdg::kDPlus, std::array{+kPiPlus, -kKPlus, +kPiPlus}, true, &sign, 2); if (indexRecDplus > -1) { flag = sign * BIT(DecayTypeMc::Ds1ToDStarK0ToDPlusPi0K0s); - fillHisto = true; } } } } + if (flag != 0) { + auto particleReso = particlesMc.iteratorAt(indexRecReso); + origin = RecoDecay::getCharmHadronOrigin(particlesMc, particleReso, false, &idxBhadMothers); + } } - if (fillHisto) { + if (flag != 0) { registry.fill(HIST("hMCRecCounter"), flag); + registry.fill(HIST("hMCRecOrigin"), origin); + } else { + if (RecoDecay::getMatchedMCRec(particlesMc, std::array{vecDaughtersReso[3], vecDaughtersReso[4]}, kK0, std::array{+kPiPlus, -kPiPlus}, true, &signV0, 2) > -1) { + SETBIT(debug, PartialMatchMc::K0Matched); + } + if (RecoDecay::getMatchedMCRec(particlesMc, std::array{vecDaughtersReso[0], vecDaughtersReso[1]}, Pdg::kD0, std::array{+kPiPlus, -kKPlus}, true, &signD0, 1) > -1) { + SETBIT(debug, PartialMatchMc::D0Matched); + } + if (RecoDecay::getMatchedMCRec(particlesMc, std::array{vecDaughtersReso[0], vecDaughtersReso[1], vecDaughtersReso[2]}, Pdg::kDStar, std::array{-kKPlus, +kPiPlus, +kPiPlus}, true, &signDStar, 2) > -1) { + SETBIT(debug, PartialMatchMc::DStarMatched); + } + if (RecoDecay::getMatchedMCRec(particlesMc, std::array{vecDaughtersReso[0], vecDaughtersReso[1], vecDaughtersReso[2]}, Pdg::kDPlus, std::array{+kPiPlus, -kKPlus, +kPiPlus}, true, &signDPlus, 2) > -1) { + SETBIT(debug, PartialMatchMc::DPlusMatched); + } registry.fill(HIST("hMCRecDebug"), debug); } - rowHfDV0McRecReduced(indexHfCandCharm, indexCandV0, flag, debug, motherPt); + rowHfDV0McRecReduced(indexHfCandCharm, indexCandV0, flag, debug, origin, motherPt); } template @@ -831,14 +842,27 @@ struct HfDataCreatorCharmResoReduced { if (RecoDecay::isMatchedMCGen(particlesMc, particle, Pdg::kDS1, std::array{static_cast(Pdg::kDStar), +kK0}, true, &sign, 1)) { registry.fill(HIST("hMCSignCounter"), sign); origin = RecoDecay::getCharmHadronOrigin(particlesMc, particle, false, &idxBhadMothers); - registry.fill(HIST("hMCOriginCounter"), origin); + registry.fill(HIST("hMCGenOrigin"), origin); auto candV0MC = particlesMc.rawIteratorAt(particle.daughtersIds().back()); // K0 -> K0s -> π+π- if (RecoDecay::isMatchedMCGen(particlesMc, candV0MC, kK0, std::array{+kPiPlus, -kPiPlus}, true, &signV0, 2)) { auto candDStarMC = particlesMc.rawIteratorAt(particle.daughtersIds().front()); // D* -> D0 π+ -> K-π+π+ - if (RecoDecay::isMatchedMCGen(particlesMc, candDStarMC, Pdg::kDStar, std::array{-kKPlus, +kPiPlus, +kPiPlus}, true, &signDStar, 2)) { - flag = signDStar * BIT(DecayTypeMc::Ds1ToDStarK0ToD0PiK0s); + if (RecoDecay::isMatchedMCGen(particlesMc, candDStarMC, Pdg::kDStar, std::array{static_cast(Pdg::kD0), +static_cast(kPiPlus)}, true, &signDStar, 1)) { + auto candD0MC = particlesMc.rawIteratorAt(candDStarMC.daughtersIds().front()); + if (RecoDecay::isMatchedMCGen(particlesMc, candDStarMC, Pdg::kDStar, std::array{-kKPlus, +kPiPlus, +kPiPlus}, true, &signDStar, 2)) { + flag = signDStar * BIT(DecayTypeMc::Ds1ToDStarK0ToD0PiK0s); + } else if (RecoDecay::isMatchedMCGen(particlesMc, candD0MC, Pdg::kD0, std::array{-kKPlus, +kPiPlus, +kPiPlus, +kPi0}, true, &signDStar, 2) || + RecoDecay::isMatchedMCGen(particlesMc, candD0MC, Pdg::kD0, std::array{-kKPlus, +kPiPlus, +kPiPlus, -kPi0}, true, &signDStar, 2)) { + flag = signDStar * BIT(DecayTypeMc::Ds1ToDStarK0ToD0PiK0sPart); + } + } else if (RecoDecay::isMatchedMCGen(particlesMc, candDStarMC, Pdg::kDStar, std::array{static_cast(Pdg::kDPlus), static_cast(kGamma)}, true, &signDStar, 1) || + RecoDecay::isMatchedMCGen(particlesMc, candDStarMC, Pdg::kDStar, std::array{static_cast(Pdg::kDPlus), -static_cast(kGamma)}, true, &signDStar, 1) || + RecoDecay::isMatchedMCGen(particlesMc, candDStarMC, Pdg::kDStar, std::array{static_cast(Pdg::kDPlus), static_cast(kPi0)}, true, &signDStar, 1) || + RecoDecay::isMatchedMCGen(particlesMc, candDStarMC, Pdg::kDStar, std::array{static_cast(Pdg::kDPlus), -static_cast(kPi0)}, true, &signDStar, 1)) { + auto candDPlusMC = particlesMc.rawIteratorAt(candDStarMC.daughtersIds().front()); + if (RecoDecay::isMatchedMCGen(particlesMc, candDPlusMC, Pdg::kDPlus, std::array{+kPiPlus, -kKPlus, +kPiPlus}, true, &signDPlus, 2)) + flag = sign * BIT(DecayTypeMc::Ds1ToDStarK0ToDPlusPi0K0s); } } } else { @@ -848,7 +872,7 @@ struct HfDataCreatorCharmResoReduced { } } // save information for task - if (!TESTBIT(std::abs(flag), DecayTypeMc::Ds1ToDStarK0ToD0PiK0s)) { + if (!TESTBIT(std::abs(flag), DecayTypeMc::Ds1ToDStarK0ToD0PiK0s) && !TESTBIT(std::abs(flag), DecayTypeMc::Ds1ToDStarK0ToD0PiK0sPart) && !TESTBIT(std::abs(flag), DecayTypeMc::Ds1ToDStarK0ToDPlusPi0K0s)) { continue; } @@ -867,14 +891,14 @@ struct HfDataCreatorCharmResoReduced { counter++; } registry.fill(HIST("hMCGenCounter"), flag, ptParticle); - rowHfResoMcGenReduced(flag, ptParticle, yParticle, etaParticle, + rowHfResoMcGenReduced(flag, origin, ptParticle, yParticle, etaParticle, ptProngs[0], yProngs[0], etaProngs[0], ptProngs[1], yProngs[1], etaProngs[1]); } else if constexpr (decayChannel == DecayChannel::DplusV0) { // Ds2Star → D+ K0 if (RecoDecay::isMatchedMCGen(particlesMc, particle, Pdg::kDS2Star, std::array{static_cast(Pdg::kDPlus), +kK0}, true, &sign, 1)) { registry.fill(HIST("hMCSignCounter"), sign); origin = RecoDecay::getCharmHadronOrigin(particlesMc, particle, false, &idxBhadMothers); - registry.fill(HIST("hMCOriginCounter"), origin); + registry.fill(HIST("hMCGenOrigin"), origin); auto candV0MC = particlesMc.rawIteratorAt(particle.daughtersIds().back()); auto candDPlusMC = particlesMc.rawIteratorAt(particle.daughtersIds().front()); // K0 -> K0s -> π+π- @@ -889,15 +913,14 @@ struct HfDataCreatorCharmResoReduced { // K0 -> K0s -> π+π- if (RecoDecay::isMatchedMCGen(particlesMc, candV0MC, kK0, std::array{+kPiPlus, -kPiPlus}, true, &signV0, 2)) { auto candDStarMC = particlesMc.rawIteratorAt(particle.daughtersIds().front()); - // D* -> D+ γ - if (RecoDecay::isMatchedMCGen(particlesMc, candDStarMC, Pdg::kDStar, std::array{static_cast(Pdg::kDPlus), static_cast(kGamma)}, true, &signDPlus, 1)) { - flag = sign * BIT(DecayTypeMc::Ds1ToDStarK0ToDPlusGammaK0s); - } else if (RecoDecay::isMatchedMCGen(particlesMc, candDStarMC, Pdg::kDStar, std::array{static_cast(Pdg::kDPlus), -static_cast(kGamma)}, true, &signDPlus, 1)) { - flag = sign * BIT(DecayTypeMc::Ds1ToDStarK0ToDPlusGammaK0s); - } else if (RecoDecay::isMatchedMCGen(particlesMc, candDStarMC, Pdg::kDStar, std::array{static_cast(Pdg::kDPlus), static_cast(kPi0)}, true, &signDPlus, 1)) { - flag = sign * BIT(DecayTypeMc::Ds1ToDStarK0ToDPlusPi0K0s); - } else if (RecoDecay::isMatchedMCGen(particlesMc, candDStarMC, Pdg::kDStar, std::array{static_cast(Pdg::kDPlus), -static_cast(kPi0)}, true, &signDPlus, 1)) { - flag = sign * BIT(DecayTypeMc::Ds1ToDStarK0ToDPlusPi0K0s); + // D* -> D+ π0/γ ->π+K-π+ π0/γ + if (RecoDecay::isMatchedMCGen(particlesMc, candDStarMC, Pdg::kDStar, std::array{static_cast(Pdg::kDPlus), static_cast(kGamma)}, true, &signDStar, 1) || + RecoDecay::isMatchedMCGen(particlesMc, candDStarMC, Pdg::kDStar, std::array{static_cast(Pdg::kDPlus), -static_cast(kGamma)}, true, &signDStar, 1) || + RecoDecay::isMatchedMCGen(particlesMc, candDStarMC, Pdg::kDStar, std::array{static_cast(Pdg::kDPlus), static_cast(kPi0)}, true, &signDStar, 1) || + RecoDecay::isMatchedMCGen(particlesMc, candDStarMC, Pdg::kDStar, std::array{static_cast(Pdg::kDPlus), -static_cast(kPi0)}, true, &signDStar, 1)) { + auto candDPlusMC = particlesMc.rawIteratorAt(candDStarMC.daughtersIds().front()); + if (RecoDecay::isMatchedMCGen(particlesMc, candDPlusMC, Pdg::kDPlus, std::array{+kPiPlus, -kKPlus, +kPiPlus}, true, &signDPlus, 2)) + flag = sign * BIT(DecayTypeMc::Ds1ToDStarK0ToDPlusPi0K0s); } } } else { @@ -908,7 +931,7 @@ struct HfDataCreatorCharmResoReduced { } } // save information for task - if (!TESTBIT(std::abs(flag), DecayTypeMc::Ds2StarToDplusK0s) && !TESTBIT(std::abs(flag), DecayTypeMc::Ds1ToDStarK0ToDPlusGammaK0s) && !TESTBIT(std::abs(flag), DecayTypeMc::Ds1ToDStarK0ToDPlusPi0K0s)) { + if (!TESTBIT(std::abs(flag), DecayTypeMc::Ds2StarToDplusK0s) && !TESTBIT(std::abs(flag), DecayTypeMc::Ds1ToDStarK0ToDPlusPi0K0s)) { continue; } @@ -927,7 +950,7 @@ struct HfDataCreatorCharmResoReduced { counter++; } registry.fill(HIST("hMCGenCounter"), flag, ptParticle); - rowHfResoMcGenReduced(flag, ptParticle, yParticle, etaParticle, + rowHfResoMcGenReduced(flag, origin, ptParticle, yParticle, etaParticle, ptProngs[0], yProngs[0], etaProngs[0], ptProngs[1], yProngs[1], etaProngs[1]); } // Dplus V0