diff --git a/PWGHF/D2H/Tasks/taskCd.cxx b/PWGHF/D2H/Tasks/taskCd.cxx index 04d586f95b1..49ae006e096 100644 --- a/PWGHF/D2H/Tasks/taskCd.cxx +++ b/PWGHF/D2H/Tasks/taskCd.cxx @@ -69,7 +69,7 @@ namespace full { // Candidate kinematics DECLARE_SOA_COLUMN(MassCd, massCd, float); //! Invariant mass of cd candidate (GeV/c^2) -DECLARE_SOA_COLUMN(MassLc, massLc, float); //! Invariant mass of lc candidate (GeV/c^2) +DECLARE_SOA_COLUMN(MassLc, massLc, float); //! Invariant mass under the p K pi hypothesis (GeV/c^2) DECLARE_SOA_COLUMN(Pt, pt, float); //! Transverse momentum of candidate (GeV/c) DECLARE_SOA_COLUMN(Eta, eta, float); //! eta of candidate (GeV/c) DECLARE_SOA_COLUMN(Phi, phi, float); //! phi of candidate (GeV/c) @@ -114,13 +114,15 @@ DECLARE_SOA_COLUMN(NTpcSignalsKa, nTpcSignalsKa, float); //! Number o DECLARE_SOA_COLUMN(NItsSignalsDe, nItsSignalsDe, float); //! Number of ITS signas DECLARE_SOA_COLUMN(CandidateSelFlag, candidateSelFlag, int8_t); //! Candidates falg DECLARE_SOA_COLUMN(CandidateSign, candidateSign, int8_t); //! Candidates sign -DECLARE_SOA_COLUMN(FlagMc, flagMc, int8_t); //! MC matching flag +DECLARE_SOA_COLUMN(FlagMc, flagMc, int8_t); //! Main MC decay-channel flag; 0 for unmatched candidates +DECLARE_SOA_COLUMN(IsCandidateSwapped, isCandidateSwapped, int8_t); //! MC-matched prong permutation; -1 for data DECLARE_SOA_COLUMN(OriginMcRec, originMcRec, int8_t); //! MC origin for reconstructed candidates DECLARE_SOA_COLUMN(FlagMcDecayChanRec, flagMcDecayChanRec, int8_t); //! Resonant MC decay channel for reconstructed candidates DECLARE_SOA_COLUMN(OriginMcGen, originMcGen, int8_t); //! MC origin for generated particles DECLARE_SOA_COLUMN(FlagMcDecayChanGen, flagMcDecayChanGen, int8_t); //! Resonant MC decay channel for generated candidates DECLARE_SOA_COLUMN(CtGen, ctGen, float); //! Generated ct computed wrt to c-deuteron production vertex, which can be either PV (prompt) or B-hadron decay vertex (non-prompt) DECLARE_SOA_COLUMN(CtRec, ctRec, float); //! Reconstructed ct computed wrt to PV +DECLARE_SOA_COLUMN(NumPvContributors, numPvContributors, uint16_t); //! Number of contributors to the primary vertex DECLARE_SOA_COLUMN(Cent, cent, float); //! Centrality DECLARE_SOA_COLUMN(VtxZ, vtxZ, float); //! Vertex Z DECLARE_SOA_COLUMN(GIndexCol, gIndexCol, int); //! Global index for the collision @@ -192,9 +194,11 @@ DECLARE_SOA_TABLE(HfCandCdFull, "AOD", "HFCANDCDFULL", full::CandidateSelFlag, full::CandidateSign, full::FlagMc, + full::IsCandidateSwapped, full::OriginMcRec, full::FlagMcDecayChanRec, full::CtGen, + full::NumPvContributors, full::Cent, full::VtxZ, full::GIndexCol, @@ -209,6 +213,7 @@ DECLARE_SOA_TABLE(HfCandCdGen, "AOD", "HFCANDCDGEN", full::OriginMcGen, full::FlagMcDecayChanGen, full::CtGen, + full::NumPvContributors, full::Cent, full::VtxZ, full::McCollisionId); @@ -231,6 +236,9 @@ struct HfTaskCd { Configurable cfgCutOnDeuteronDcaOrdering{"cfgCutOnDeuteronDcaOrdering", false, "Require deuteron DCA to be smaller than kaon and pion DCAs"}; Configurable cfgMinDeuteronDcaPreselection{"cfgMinDeuteronDcaPreselection", 0.004, "Minimum deuteron DCA for preselection (cm)"}; Configurable cfgMaxDeuteronTofPidPreselection{"cfgMaxDeuteronTofPidPreselection", 5, "Maximum |nSigma TOF| for deuteron preselection"}; + Configurable fillMcCorrelatedBackgrounds{"fillMcCorrelatedBackgrounds", false, "Store selected MC candidates matched to a supported decay channel other than c-deuteron"}; + Configurable fillMcCombinatorialBackground{"fillMcCombinatorialBackground", false, "Store selected MC candidates with flagMcMatchRec equal to zero"}; + Configurable acceptCandidatesWithoutCdFlag{"acceptCandidatesWithoutCdFlag", false, "Accept selected candidates without the Cd skim bit; intended for MC reflection studies"}; SliceCache cache; @@ -263,9 +271,9 @@ struct HfTaskCd { ConfigurableAxis thnConfigAxisNumPvContr{"thnConfigAxisNumPvContr", {200, -0.5, 199.5}, "Number of PV contributors"}; ConfigurableAxis thnConfigAxisCt{"thnConfigAxisCt", {500, 0., 5000.}, ""}; - constexpr static std::string_view SignalFolders[] = {"signal", "prompt", "nonprompt"}; - constexpr static std::string_view SignalSuffixes[] = {"", "Prompt", "NonPrompt"}; - const float cmToMum = 1.e4; + static constexpr std::array SignalFolders = {"signal", "prompt", "nonprompt"}; + static constexpr std::array SignalSuffixes = {"", "Prompt", "NonPrompt"}; + static constexpr float CmToMum = 1.e4f; enum SignalClasses : int { Signal = 0, @@ -332,6 +340,9 @@ struct HfTaskCd { addHistogramsGen("hPt", "#it{p}_{T}^{gen.} (GeV/#it{c})", "entries", {HistType::kTH1D, {{360, 0., 36.}}}); if (!isData) { registry.add("MC/generated/signal/hPtGenSig", "3-prong candidates (matched);#it{p}_{T}^{gen.} (GeV/#it{c});entries", {HistType::kTH1D, {{360, 0., 36.}}}); + registry.add("MC/reconstructed/allCandidates/hFlagMcMatchRec", "Selected 3-prong candidates;MC decay-channel flag;entries", {HistType::kTH1D, {{49, -24.5, 24.5}}}); + registry.add("MC/reconstructed/allCandidates/hMassCdVsFlagMcMatchRec", "Selected 3-prong candidates;inv. mass (d K #pi) (GeV/#it{c}^{2});MC decay-channel flag", {HistType::kTH2F, {{400, 2.4, 4.4}, {49, -24.5, 24.5}}}); + registry.add("MC/reconstructed/allCandidates/hMassPKPiVsFlagMcMatchRec", "Selected 3-prong candidates;inv. mass (p K #pi) (GeV/#it{c}^{2});MC decay-channel flag", {HistType::kTH2F, {{400, 1.6, 3.2}, {49, -24.5, 24.5}}}); } addHistogramsRec("hPtProng0", "prong 0 #it{p}_{T} (GeV/#it{c})", "entries", {HistType::kTH1D, {{360, 0., 36.}}}); addHistogramsRec("hPtProng1", "prong 1 #it{p}_{T} (GeV/#it{c})", "entries", {HistType::kTH1D, {{360, 0., 36.}}}); @@ -496,10 +507,14 @@ struct HfTaskCd { const int64_t timeStamp = bc.timestamp(); for (const auto& candidate : groupedCdCandidates) { - if (candidate.flagMcMatchRec() == 0) { // we skip combinatorial background + const int absFlagMc = std::abs(candidate.flagMcMatchRec()); + const bool isTrueCd = absFlagMc == hf_decay::hf_cand_3prong::DecayChannelMain::CDeuteronToDeKPi; + const bool keepCorrelatedBackground = fillMcCorrelatedBackgrounds && absFlagMc != 0 && !isTrueCd; + const bool keepCombinatorialBackground = fillMcCombinatorialBackground && absFlagMc == 0; + if (!isTrueCd && !keepCorrelatedBackground && !keepCombinatorialBackground) { continue; } - if (!TESTBIT(candidate.hfflag(), aod::hf_cand_3prong::DecayType::CdToDeKPi)) { + if (!acceptCandidatesWithoutCdFlag && !TESTBIT(candidate.hfflag(), aod::hf_cand_3prong::DecayType::CdToDeKPi)) { continue; } const auto yCd = RecoDecay::y(candidate.pVector(), o2::constants::physics::MassCDeuteron); @@ -507,14 +522,24 @@ struct HfTaskCd { continue; } + registry.fill(HIST("MC/reconstructed/allCandidates/hFlagMcMatchRec"), candidate.flagMcMatchRec()); + if (candidate.isSelCdToDeKPi() >= selectionFlagCd) { + registry.fill(HIST("MC/reconstructed/allCandidates/hMassCdVsFlagMcMatchRec"), HfHelper::invMassCdToDeKPi(candidate), candidate.flagMcMatchRec()); + registry.fill(HIST("MC/reconstructed/allCandidates/hMassPKPiVsFlagMcMatchRec"), HfHelper::invMassLcToPKPi(candidate), candidate.flagMcMatchRec()); + } + if (candidate.isSelCdToPiKDe() >= selectionFlagCd) { + registry.fill(HIST("MC/reconstructed/allCandidates/hMassCdVsFlagMcMatchRec"), HfHelper::invMassCdToPiKDe(candidate), candidate.flagMcMatchRec()); + registry.fill(HIST("MC/reconstructed/allCandidates/hMassPKPiVsFlagMcMatchRec"), HfHelper::invMassLcToPiKP(candidate), candidate.flagMcMatchRec()); + } + float ctGen{-1.f}, ptGen{-1.f}; int pdgCodeProng0{0}; - if (std::abs(candidate.flagMcMatchRec()) == hf_decay::hf_cand_3prong::DecayChannelMain::CDeuteronToDeKPi) { + if (isTrueCd) { const auto& mcParticleProng0 = candidate.template prong0_as().template mcParticle_as(); pdgCodeProng0 = std::abs(mcParticleProng0.pdgCode()); const auto indexMother = RecoDecay::getMother(mcParticles, mcParticleProng0, o2::constants::physics::Pdg::kCDeuteron, true); const auto particleMother = mcParticles.rawIteratorAt(indexMother); - ctGen = RecoDecay::ct(std::array{particleMother.px(), particleMother.py(), particleMother.pz()}, RecoDecay::distance(std::array{particleMother.vx(), particleMother.vy(), particleMother.vz()}, std::array{mcParticleProng0.vx(), mcParticleProng0.vy(), mcParticleProng0.vz()}), o2::constants::physics::MassCDeuteron) * cmToMum; + ctGen = RecoDecay::ct(std::array{particleMother.px(), particleMother.py(), particleMother.pz()}, RecoDecay::distance(std::array{particleMother.vx(), particleMother.vy(), particleMother.vz()}, std::array{mcParticleProng0.vx(), mcParticleProng0.vy(), mcParticleProng0.vz()}), o2::constants::physics::MassCDeuteron) * CmToMum; ptGen = particleMother.pt(); } @@ -566,7 +591,7 @@ struct HfTaskCd { candidate.decayLength(), candidate.cpa(), candidate.chi2PCA(), nSigmaTpcDe, nSigmaTpcPr, nSigmaItsDe, nSigmaTofDe, tofBetaDe, tpcInnerParamDe, tofExpMomDe, - candidate.ct(o2::constants::physics::MassCDeuteron) * cmToMum, + candidate.ct(o2::constants::physics::MassCDeuteron) * CmToMum, candFlag, candSign, candidate.flagMcMatchRec(), candidate.originMcRec(), candidate.flagMcDecayChanRec(), ctGen, o2::hf_centrality::getCentralityColl(collision)); } @@ -581,23 +606,34 @@ struct HfTaskCd { nSigmaTpcDe, nSigmaTpcPr, nSigmaItsDe, nSigmaTofDe, tofBetaDe, tpcInnerParamDe, tofExpMomDe, nSigmaTpcPi, nSigmaTofPi, nSigmaTpcKa, nSigmaTofKa, - candidate.ct(o2::constants::physics::MassCDeuteron) * cmToMum, - candFlag, candSign, candidate.flagMcMatchRec(), candidate.originMcRec(), - candidate.flagMcDecayChanRec(), ctGen, o2::hf_centrality::getCentralityColl(collision), + candidate.ct(o2::constants::physics::MassCDeuteron) * CmToMum, + candFlag, candSign, candidate.flagMcMatchRec(), candidate.isCandidateSwapped(), candidate.originMcRec(), + candidate.flagMcDecayChanRec(), ctGen, collision.numContrib(), o2::hf_centrality::getCentralityColl(collision), collision.posZ(), collision.globalIndex(), timeStamp); } }; - // For MC, retain only the hypothesis that matches the true deuteron prong. - if (selDeKPi && pdgCodeProng0 == o2::constants::physics::Pdg::kDeuteron) { - writeMcHypothesis(true); - } - if (selPiKDe && pdgCodeProng0 == kPiPlus) { - writeMcHypothesis(false); + if (isTrueCd) { + // For signal MC, retain only the hypothesis that matches the true deuteron prong. + if (selDeKPi && pdgCodeProng0 == o2::constants::physics::Pdg::kDeuteron) { + writeMcHypothesis(true); + } + if (selPiKDe && pdgCodeProng0 == kPiPlus) { + writeMcHypothesis(false); + } + } else { + // Reflection and combinatorial candidates have no true deuteron prong. + // Retain every Cd hypothesis accepted by the selector and distinguish them with CandidateSelFlag. + if (selDeKPi) { + writeMcHypothesis(true); + } + if (selPiKDe) { + writeMcHypothesis(false); + } } } - if (std::abs(candidate.flagMcMatchRec()) != hf_decay::hf_cand_3prong::DecayChannelMain::CDeuteronToDeKPi) { + if (!isTrueCd) { continue; } @@ -661,7 +697,7 @@ struct HfTaskCd { const bool isPrompt = particle.originMcGen() == 0; const float ptGenB = isPrompt ? -1.f : mcParticles.rawIteratorAt(particle.idxBhadMotherPart()).pt(); const auto firstDau = particle.template daughters_as().begin(); - const float ctGen = RecoDecay::ct(std::array{particle.px(), particle.py(), particle.pz()}, RecoDecay::distance(std::array{particle.vx(), particle.vy(), particle.vz()}, std::array{firstDau.vx(), firstDau.vy(), firstDau.vz()}), o2::constants::physics::MassCDeuteron) * cmToMum; + const float ctGen = RecoDecay::ct(std::array{particle.px(), particle.py(), particle.pz()}, RecoDecay::distance(std::array{particle.vx(), particle.vy(), particle.vz()}, std::array{firstDau.vx(), firstDau.vy(), firstDau.vz()}), o2::constants::physics::MassCDeuteron) * CmToMum; fillHistogramsGen(particle, yGen); if (isPrompt) { @@ -684,6 +720,7 @@ struct HfTaskCd { particle.originMcGen(), particle.flagMcDecayChanGen(), ctGen, + numPvContributors, cent, vtxZ, particle.mcCollision().globalIndex()); @@ -849,7 +886,7 @@ struct HfTaskCd { tofBetaDe, tpcInnerParamDe, tofExpMomDe, nSigmaTpcPi, nSigmaTofPi, nSigmaTpcKa, nSigmaTofKa, candidate.ct(o2::constants::physics::MassCDeuteron), - candFlag, candSign, 0, 0, -1, -1.f, cent, + candFlag, candSign, 0, -1, 0, 0, -1.f, collision.numContrib(), cent, collision.posZ(), collision.globalIndex(), timeStamp); } }; diff --git a/PWGHF/TableProducer/candidateSelectorCd.cxx b/PWGHF/TableProducer/candidateSelectorCd.cxx index 61c662ada3e..7658c95cf11 100644 --- a/PWGHF/TableProducer/candidateSelectorCd.cxx +++ b/PWGHF/TableProducer/candidateSelectorCd.cxx @@ -71,10 +71,11 @@ struct HfCandidateSelectorCd { Configurable itsChi2PerClusterMax{"itsChi2PerClusterMax", 1e10f, "max its fit chi2 per ITS cluster"}; // DCA track cuts Configurable> binsPtTrack{"binsPtTrack", std::vector{hf_cuts_single_track::vecBinsPtTrack}, "track pT bin limits for DCA XY/Z pT-dependent cut"}; - Configurable> cutsSingleTrack{"cutsSingleTrack", {hf_cuts_single_track::CutsTrack[0], hf_cuts_single_track::NBinsPtTrack, hf_cuts_single_track::NCutVarsTrack, hf_cuts_single_track::labelsPtTrack, hf_cuts_single_track::labelsCutVarTrack}, "Single-track selections"}; + Configurable> cutsSingleTrack{"cutsSingleTrack", {&hf_cuts_single_track::CutsTrack[0][0], hf_cuts_single_track::NBinsPtTrack, hf_cuts_single_track::NCutVarsTrack, hf_cuts_single_track::labelsPtTrack, hf_cuts_single_track::labelsCutVarTrack}, "Single-track selections"}; // topological cuts Configurable> binsPt{"binsPt", std::vector{hf_cuts_cd_to_de_k_pi::vecBinsPt}, "pT bin limits"}; - Configurable> cuts{"cuts", {hf_cuts_cd_to_de_k_pi::Cuts[0], hf_cuts_cd_to_de_k_pi::NBinsPt, hf_cuts_cd_to_de_k_pi::NCutVars, hf_cuts_cd_to_de_k_pi::labelsPt, hf_cuts_cd_to_de_k_pi::labelsCutVar}, "Cd candidate selection per pT bin"}; + Configurable> cuts{"cuts", {&hf_cuts_cd_to_de_k_pi::Cuts[0][0], hf_cuts_cd_to_de_k_pi::NBinsPt, hf_cuts_cd_to_de_k_pi::NCutVars, hf_cuts_cd_to_de_k_pi::labelsPt, hf_cuts_cd_to_de_k_pi::labelsCutVar}, "Cd candidate selection per pT bin"}; + Configurable acceptCandidatesWithoutCdFlag{"acceptCandidatesWithoutCdFlag", false, "Apply the Cd mass hypotheses also to 3-prong candidates without the Cd skim bit; intended for MC reflection studies"}; // QA switch Configurable activateQA{"activateQA", false, "Flag to enable QA histogram"}; @@ -102,15 +103,15 @@ struct HfCandidateSelectorCd { selectorDeuteron = selectorPion; if (activateQA) { - constexpr int kNBinsSelections = aod::SelectionStep::NSelectionSteps; - std::string labels[kNBinsSelections]; + constexpr int NBinsSelections = aod::SelectionStep::NSelectionSteps; + std::string labels[NBinsSelections]; labels[0] = "No selection"; labels[1 + aod::SelectionStep::RecoSkims] = "Skims selection"; labels[1 + aod::SelectionStep::RecoTopol] = "Skims & Topological selections"; labels[1 + aod::SelectionStep::RecoPID] = "Skims & Topological & PID selections"; - static const AxisSpec axisSelections = {kNBinsSelections, 0.5, kNBinsSelections + 0.5, ""}; + static const AxisSpec axisSelections = {NBinsSelections, 0.5, NBinsSelections + 0.5, ""}; registry.add("hSelections", "Selections;;#it{p}_{T} (GeV/#it{c})", {HistType::kTH2F, {axisSelections, {(std::vector)binsPt, "#it{p}_{T} (GeV/#it{c})"}}}); - for (int iBin = 0; iBin < kNBinsSelections; ++iBin) { + for (int iBin = 0; iBin < NBinsSelections; ++iBin) { registry.get(HIST("hSelections"))->GetXaxis()->SetBinLabel(iBin + 1, labels[iBin].data()); } } @@ -264,7 +265,7 @@ struct HfCandidateSelectorCd { auto ptCand = candidate.pt(); - if (!(candidate.hfflag() & 1 << aod::hf_cand_3prong::DecayType::CdToDeKPi)) { + if (!acceptCandidatesWithoutCdFlag && !(candidate.hfflag() & 1 << aod::hf_cand_3prong::DecayType::CdToDeKPi)) { hfSelCdCandidate(statusCdToDeKPi, statusCdToPiKDe); if (activateQA) { registry.fill(HIST("hSelections"), 1, ptCand); @@ -314,24 +315,11 @@ struct HfCandidateSelectorCd { if (usePid) { // track-level PID selection - TrackSelectorPID::Status pidTrackPos1Deuteron; - TrackSelectorPID::Status pidTrackPos2Deuteron; - TrackSelectorPID::Status pidTrackPos1Pion; - TrackSelectorPID::Status pidTrackPos2Pion; - TrackSelectorPID::Status pidTrackNegKaon; - if (usePidTpcAndTof) { - pidTrackPos1Deuteron = selectorDeuteron.statusTpcAndTof(trackPos1, candidate.nSigTpcDe0(), candidate.nSigTofDe0()); - pidTrackPos2Deuteron = selectorDeuteron.statusTpcAndTof(trackPos2, candidate.nSigTpcDe2(), candidate.nSigTofDe2()); - pidTrackPos1Pion = selectorPion.statusTpcAndTof(trackPos1, candidate.nSigTpcPi0(), candidate.nSigTofPi0()); - pidTrackPos2Pion = selectorPion.statusTpcAndTof(trackPos2, candidate.nSigTpcPi2(), candidate.nSigTofPi2()); - pidTrackNegKaon = selectorKaon.statusTpcAndTof(trackNeg, candidate.nSigTpcKa1(), candidate.nSigTofKa1()); - } else { - pidTrackPos1Deuteron = selectorDeuteron.statusTpcOrTof(trackPos1, candidate.nSigTpcDe0(), candidate.nSigTofDe0()); - pidTrackPos2Deuteron = selectorDeuteron.statusTpcOrTof(trackPos2, candidate.nSigTpcDe2(), candidate.nSigTofDe2()); - pidTrackPos1Pion = selectorPion.statusTpcOrTof(trackPos1, candidate.nSigTpcPi0(), candidate.nSigTofPi0()); - pidTrackPos2Pion = selectorPion.statusTpcOrTof(trackPos2, candidate.nSigTpcPi2(), candidate.nSigTofPi2()); - pidTrackNegKaon = selectorKaon.statusTpcOrTof(trackNeg, candidate.nSigTpcKa1(), candidate.nSigTofKa1()); - } + const auto pidTrackPos1Deuteron = usePidTpcAndTof ? selectorDeuteron.statusTpcAndTof(trackPos1, candidate.nSigTpcDe0(), candidate.nSigTofDe0()) : selectorDeuteron.statusTpcOrTof(trackPos1, candidate.nSigTpcDe0(), candidate.nSigTofDe0()); + const auto pidTrackPos2Deuteron = usePidTpcAndTof ? selectorDeuteron.statusTpcAndTof(trackPos2, candidate.nSigTpcDe2(), candidate.nSigTofDe2()) : selectorDeuteron.statusTpcOrTof(trackPos2, candidate.nSigTpcDe2(), candidate.nSigTofDe2()); + const auto pidTrackPos1Pion = usePidTpcAndTof ? selectorPion.statusTpcAndTof(trackPos1, candidate.nSigTpcPi0(), candidate.nSigTofPi0()) : selectorPion.statusTpcOrTof(trackPos1, candidate.nSigTpcPi0(), candidate.nSigTofPi0()); + const auto pidTrackPos2Pion = usePidTpcAndTof ? selectorPion.statusTpcAndTof(trackPos2, candidate.nSigTpcPi2(), candidate.nSigTofPi2()) : selectorPion.statusTpcOrTof(trackPos2, candidate.nSigTpcPi2(), candidate.nSigTofPi2()); + const auto pidTrackNegKaon = usePidTpcAndTof ? selectorKaon.statusTpcAndTof(trackNeg, candidate.nSigTpcKa1(), candidate.nSigTofKa1()) : selectorKaon.statusTpcOrTof(trackNeg, candidate.nSigTpcKa1(), candidate.nSigTofKa1()); if (!isSelectedPID(pidTrackPos1Deuteron, pidTrackNegKaon, pidTrackPos2Pion)) { pidCdToDeKPi = 0; // reject CdToDeKPi