From 8a4d1e29ea5f3c4c32e958d269a54fdf43360f08 Mon Sep 17 00:00:00 2001 From: fchinu Date: Tue, 9 Jun 2026 18:40:21 +0200 Subject: [PATCH] Implement online workflow for B0->J/Psi K*0 --- PWGHF/D2H/DataModel/ReducedDataModel.h | 39 ++- .../candidateCreatorBToJpsiReduced.cxx | 2 +- .../dataCreatorJpsiHadReduced.cxx | 287 +++++++++++++++++- .../DataModel/CandidateReconstructionTables.h | 67 +++- 4 files changed, 373 insertions(+), 22 deletions(-) diff --git a/PWGHF/D2H/DataModel/ReducedDataModel.h b/PWGHF/D2H/DataModel/ReducedDataModel.h index 09f8bcd2d1b..9499fc76a31 100644 --- a/PWGHF/D2H/DataModel/ReducedDataModel.h +++ b/PWGHF/D2H/DataModel/ReducedDataModel.h @@ -749,14 +749,17 @@ using HfRedPidDau2 = HfRedPidDau2s::iterator; // Beauty candidates prongs namespace hf_cand_b0_reduced { -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_INDEX_COLUMN_FULL(ProngD0, prongD0, int, HfRed2Prongs, "_0"); //! ProngD0 index -DECLARE_SOA_INDEX_COLUMN_FULL(ProngBachPi, prongBachPi, int, HfRedTrackBases, "_1"); //! ProngBachPi index -DECLARE_SOA_INDEX_COLUMN_FULL(ProngSoftPi, prongSoftPi, int, HfRedSoftPiBases, "_2"); //! ProngSoftPi index -DECLARE_SOA_COLUMN(Prong0MlScoreBkg, prong0MlScoreBkg, float); //! Bkg ML score of the D daughter -DECLARE_SOA_COLUMN(Prong0MlScorePrompt, prong0MlScorePrompt, float); //! Prompt ML score of the D daughter -DECLARE_SOA_COLUMN(Prong0MlScoreNonprompt, prong0MlScoreNonprompt, float); //! Nonprompt ML score of the D daughter +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_INDEX_COLUMN_FULL(ProngD0, prongD0, int, HfRed2Prongs, "_0"); //! ProngD0 index +DECLARE_SOA_INDEX_COLUMN_FULL(ProngBachPi, prongBachPi, int, HfRedTrackBases, "_1"); //! ProngBachPi index +DECLARE_SOA_INDEX_COLUMN_FULL(ProngSoftPi, prongSoftPi, int, HfRedSoftPiBases, "_2"); //! ProngSoftPi index +DECLARE_SOA_INDEX_COLUMN_FULL(Jpsi, jpsi, int, HfRedJpsis, "_0"); //! J/Psi index +DECLARE_SOA_INDEX_COLUMN_FULL(Prong0K0Star, prong0K0Star, int, HfRedBach0Bases, "_0"); //! J/Psi index +DECLARE_SOA_INDEX_COLUMN_FULL(Prong1K0Star, prong1K0Star, int, HfRedBach1Bases, "_0"); //! J/Psi index +DECLARE_SOA_COLUMN(Prong0MlScoreBkg, prong0MlScoreBkg, float); //! Bkg ML score of the D daughter +DECLARE_SOA_COLUMN(Prong0MlScorePrompt, prong0MlScorePrompt, float); //! Prompt ML score of the D daughter +DECLARE_SOA_COLUMN(Prong0MlScoreNonprompt, prong0MlScoreNonprompt, float); //! Nonprompt ML score of the D daughter } // namespace hf_cand_b0_reduced DECLARE_SOA_TABLE(HfRedB0Prongs, "AOD", "HFREDB0PRONG", //! Table with B0 daughter indices @@ -904,6 +907,16 @@ DECLARE_SOA_TABLE(HfMcRecRedDStarPis, "AOD", "HFMCRECREDDSTPI", //! Table with r hf_cand_mc_flag::DebugMcRec, hf_b0_mc::PtMother); +DECLARE_SOA_TABLE(HfMcRecRedJPK0ss, "AOD", "HFMCRECREDJPK0S", //! Table with reconstructed MC information on J/Psi/K*0(<-B0) pairs for reduced workflow + hf_cand_b0_reduced::JpsiId, + hf_cand_b0_reduced::Prong0K0StarId, + hf_cand_b0_reduced::Prong1K0StarId, + hf_cand_mc_flag::FlagMcMatchRec, + hf_cand_mc_flag::FlagMcDecayChanRec, + hf_cand_mc_flag::FlagWrongCollision, + hf_cand_mc_flag::DebugMcRec, + hf_b0_mc::PtMother); + // Table with same size as HFCANDB0 DECLARE_SOA_TABLE(HfMcRecRedB0s, "AOD", "HFMCRECREDB0", //! Reconstruction-level MC information on B0 candidates for reduced workflow hf_cand_mc_flag::FlagMcMatchRec, @@ -941,14 +954,18 @@ DECLARE_SOA_TABLE(HfMcGenRedB0s, "AOD", "HFMCGENREDB0", //! Generation-level MC // so we can use them in the B0 part namespace hf_cand_b0_config { -DECLARE_SOA_COLUMN(MySelectionFlagD, mySelectionFlagD, int8_t); //! Flag to filter selected D+ mesons -DECLARE_SOA_COLUMN(MyInvMassWindowDPi, myInvMassWindowDPi, float); //! Half-width of the B0 invariant-mass window in GeV/c2 +DECLARE_SOA_COLUMN(MySelectionFlagD, mySelectionFlagD, int8_t); //! Flag to filter selected D+ mesons +DECLARE_SOA_COLUMN(MyInvMassWindowDPi, myInvMassWindowDPi, float); //! Half-width of the B0 invariant-mass window in GeV/c2 +DECLARE_SOA_COLUMN(MyInvMassWindowJpsiK0Star, myInvMassWindowJpsiK0Star, float); //! Half-width of the B0 invariant-mass window in GeV/c2 } // namespace hf_cand_b0_config DECLARE_SOA_TABLE(HfCandB0Configs, "AOD", "HFCANDB0CONFIG", //! Table with configurables information for reduced workflow hf_cand_b0_config::MySelectionFlagD, hf_cand_b0_config::MyInvMassWindowDPi); +DECLARE_SOA_TABLE(HfCfgB0ToJpsis, "AOD", "HFCFGB0TOJPSI", //! Table with configurables information for reduced workflow + hf_cand_b0_config::MyInvMassWindowJpsiK0Star); + namespace hf_bplus_mc { // MC Rec @@ -1046,7 +1063,7 @@ DECLARE_SOA_TABLE(HfCandBpConfigs, "AOD", "HFCANDBPCONFIG", //! Table with confi hf_cand_bplus_config::MySelectionFlagD0bar, hf_cand_bplus_config::MyInvMassWindowD0Pi); -DECLARE_SOA_TABLE(HfCfgBpToJpsi, "AOD", "HFCFGBPTOJPSI", //! Table with configurables information for reduced workflow +DECLARE_SOA_TABLE(HfCfgBpToJpsis, "AOD", "HFCFGBPTOJPSI", //! Table with configurables information for reduced workflow hf_cand_bplus_config::MyInvMassWindowJpsiK); namespace hf_bs_mc diff --git a/PWGHF/D2H/TableProducer/candidateCreatorBToJpsiReduced.cxx b/PWGHF/D2H/TableProducer/candidateCreatorBToJpsiReduced.cxx index 9a2a6621fd6..b42d84c54fc 100644 --- a/PWGHF/D2H/TableProducer/candidateCreatorBToJpsiReduced.cxx +++ b/PWGHF/D2H/TableProducer/candidateCreatorBToJpsiReduced.cxx @@ -354,7 +354,7 @@ struct HfCandidateCreatorBToJpsiReduced { soa::Join const& candsJpsi, soa::Join const& tracksKaon, aod::HfOrigColCounts const& collisionsCounter, - aod::HfCfgBpToJpsi const& configs) + aod::HfCfgBpToJpsis const& configs) { // Jpsi K invariant-mass window cut for (const auto& config : configs) { diff --git a/PWGHF/D2H/TableProducer/dataCreatorJpsiHadReduced.cxx b/PWGHF/D2H/TableProducer/dataCreatorJpsiHadReduced.cxx index a8b7df7303c..83716354e3b 100644 --- a/PWGHF/D2H/TableProducer/dataCreatorJpsiHadReduced.cxx +++ b/PWGHF/D2H/TableProducer/dataCreatorJpsiHadReduced.cxx @@ -121,11 +121,14 @@ struct HfDataCreatorJpsiHadReduced { Produces hfTrackCovLfDau1; // MC related tables Produces rowHfJpsiKMcRecReduced; + Produces rowHfJpsiK0StarMcRecReduced; Produces rowHfJpsiPhiMcRecReduced; Produces rowHfBpMcGenReduced; + Produces rowHfB0McGenReduced; Produces rowHfBsMcGenReduced; - Produces rowCandidateConfigBplus; + Produces rowCandidateConfigBplus; + Produces rowCandidateConfigB0; Produces rowCandidateConfigBs; Configurable skipRejectedCollisions{"skipRejectedCollisions", true, "skips collisions rejected by the event selection, instead of flagging only"}; @@ -156,6 +159,7 @@ struct HfDataCreatorJpsiHadReduced { Configurable> cuts{"cuts", {hf_cuts_jpsi_to_mu_mu::Cuts[0], hf_cuts_jpsi_to_mu_mu::NBinsPt, hf_cuts_jpsi_to_mu_mu::NCutVars, hf_cuts_jpsi_to_mu_mu::labelsPt, hf_cuts_jpsi_to_mu_mu::labelsCutVar}, "J/Psi candidate selection per pT bin"}; Configurable invMassWindowJpsiHad{"invMassWindowJpsiHad", 0.3, "invariant-mass window for Jpsi-Had pair preselections (GeV/c2)"}; Configurable deltaMPhiMax{"deltaMPhiMax", 0.02, "invariant-mass window for phi preselections (GeV/c2) (only for Bs->J/PsiPhi)"}; + Configurable deltaMK0StarMax{"deltaMK0StarMax", 0.02, "invariant-mass window for K*0 preselections (GeV/c2) (only for B0->J/PsiK0*)"}; Configurable cpaMin{"cpaMin", 0., "Minimum cosine of pointing angle for B candidates"}; Configurable decLenMin{"decLenMin", 0., "Minimum decay length for B candidates"}; @@ -202,7 +206,7 @@ struct HfDataCreatorJpsiHadReduced { selectorElectron.setRangePtTpc(selectionsPid.ptPidTpcMin, selectionsPid.ptPidTpcMax); selectorElectron.setRangeNSigmaTpc(selectionsPid.nSigmaTpcElMinForVeto, selectionsPid.nSigmaTpcElMaxForVeto); - std::array doProcess = {doprocessJpsiKData, doprocessJpsiKMc, doprocessJpsiPhiData, doprocessJpsiPhiMc}; + std::array doProcess = {doprocessJpsiKData, doprocessJpsiKMc, doprocessJpsiK0StarData, doprocessJpsiK0StarMc, doprocessJpsiPhiData, doprocessJpsiPhiMc}; if (std::accumulate(doProcess.begin(), doProcess.end(), 0) != 1) { LOGP(fatal, "One and only one process function can be enabled at a time, please fix your configuration!"); } @@ -236,13 +240,21 @@ struct HfDataCreatorJpsiHadReduced { registry.add("hCpaJpsi", "J/Psi cos#theta_{p};J/Psi cos#theta_{p};Counts", {HistType::kTH1F, {{200, -1., 1, "J/Psi cos#theta_{p}"}}}); std::shared_ptr hFitCandidatesJpsi = registry.add("hFitCandidatesJpsi", "Jpsi candidate counter", {HistType::kTH1D, {axisCands}}); std::shared_ptr hFitCandidatesBPlus = registry.add("hFitCandidatesBPlus", "hFitCandidatesBPlus candidate counter", {HistType::kTH1D, {axisCands}}); + std::shared_ptr hFitCandidatesB0 = registry.add("hFitCandidatesB0", "hFitCandidatesB0 candidate counter", {HistType::kTH1D, {axisCands}}); std::shared_ptr hFitCandidatesBS = registry.add("hFitCandidatesBS", "hFitCandidatesBS candidate counter", {HistType::kTH1D, {axisCands}}); setLabelHistoCands(hFitCandidatesJpsi); setLabelHistoCands(hFitCandidatesBPlus); + setLabelHistoCands(hFitCandidatesB0); setLabelHistoCands(hFitCandidatesBS); if (doprocessJpsiKData || doprocessJpsiKMc) { registry.add("hPtKaon", "Kaon #it{p}_{T};#it{p}_{T} (GeV/#it{c});Counts", {HistType::kTH1F, {{100, 0., 10.}}}); registry.add("hMassJpsiKaon", "J/Psi Kaon mass;#it{M}_{J/#PsiK} (GeV/#it{c}^{2});Counts", {HistType::kTH1F, {{800, 4.9, 5.7}}}); + } else if (doprocessJpsiK0StarData || doprocessJpsiK0StarMc) { + registry.add("hPtK0Star", "K*0 #it{p}_{T};#it{p}_{T} (GeV/#it{c});Counts", {HistType::kTH1F, {{100, 0., 10.}}}); + registry.add("hMassK0Star", "K*0 mass;#it{M}_{#piK} (GeV/#it{c}^{2});Counts", {HistType::kTH1F, {{200, 0.9, 1.2}}}); + registry.add("hMassJpsiK0Star", "J/Psi K*0 mass;#it{M}_{J/#PsiK*0} (GeV/#it{c}^{2});Counts", {HistType::kTH1F, {{800, 4.9, 5.7}}}); + std::shared_ptr hFitCandidatesK0Star = registry.add("hFitCandidatesK0Star", "K*0 candidate counter", {HistType::kTH1D, {axisCands}}); + setLabelHistoCands(hFitCandidatesK0Star); } else if (doprocessJpsiPhiData || doprocessJpsiPhiMc) { registry.add("hPtPhi", "Phi #it{p}_{T};#it{p}_{T} (GeV/#it{c});Counts", {HistType::kTH1F, {{100, 0., 10.}}}); registry.add("hMassPhi", "Phi mass;#it{M}_{KK} (GeV/#it{c}^{2});Counts", {HistType::kTH1F, {{200, 0.9, 1.2}}}); @@ -492,7 +504,36 @@ struct HfDataCreatorJpsiHadReduced { checkWrongCollision(particleMother, collision, indexCollisionMaxNumContrib, flagWrongCollision); } } - rowHfJpsiKMcRecReduced(indexHfCandJpsi, selectedTracksBach[0][vecDaughtersB.back().globalIndex()], flag, channel, flagWrongCollision, debug, motherPt); + rowHfJpsiKMcRecReduced(indexHfCandJpsi, selectedTracksBach[0][vecDaughtersB[2].globalIndex()], flag, channel, flagWrongCollision, debug, motherPt); + } else if constexpr (DecChannel == DecayChannel::B0ToJpsiK0Star) { + // B0 → J/Psi K0* → (µ+µ-) (K+pi-) + int indexRec = -1; + indexRec = RecoDecay::getMatchedMCRec(particlesMc, std::array{vecDaughtersB[0], vecDaughtersB[1], vecDaughtersB[2], vecDaughtersB[3]}, Pdg::kB0, std::array{-kMuonMinus, +kMuonMinus, +kKPlus, -kPiPlus}, true, &sign, 4); + if (indexRec > -1) { + // J/Psi → µ+µ- + indexRec = RecoDecay::getMatchedMCRec(particlesMc, std::array{vecDaughtersB[0], vecDaughtersB[1]}, Pdg::kJPsi, std::array{-kMuonMinus, +kMuonMinus}, true, &sign, 1); + if (indexRec > -1) { + flag = sign * o2::hf_decay::hf_cand_beauty::B0ToJpsiPiK; + indexRec = RecoDecay::getMatchedMCRec(particlesMc, std::array{vecDaughtersB[2], vecDaughtersB[3]}, Pdg::kK0Star892, std::array{+kKPlus, -kPiPlus}, true, &sign, 1); + if (indexRec > -1) { + channel = o2::hf_decay::hf_cand_beauty::B0ToJpsiKstar0; + } else { + debug = 1; + LOGF(debug, "B0 decays in the expected final state but the condition on the K0* intermediate state is not fulfilled"); + } + } else { + debug = 1; + LOGF(debug, "B0 decays in the expected final state but the condition on the J/Psi intermediate state is not fulfilled"); + } + + auto indexMother = RecoDecay::getMother(particlesMc, vecDaughtersB.back().template mcParticle_as(), Pdg::kB0, true); + if (indexMother >= 0) { + auto particleMother = particlesMc.rawIteratorAt(indexMother); + motherPt = particleMother.pt(); + checkWrongCollision(particleMother, collision, indexCollisionMaxNumContrib, flagWrongCollision); + } + } + rowHfJpsiK0StarMcRecReduced(indexHfCandJpsi, selectedTracksBach[0][vecDaughtersB[2].globalIndex()], selectedTracksBach[1][vecDaughtersB[3].globalIndex()], flag, channel, flagWrongCollision, debug, motherPt); } else if constexpr (DecChannel == DecayChannel::BsToJpsiPhi) { // Bs → J/Psi phi → (µ+µ-) (K+K-) int indexRec = -1; @@ -521,7 +562,7 @@ struct HfDataCreatorJpsiHadReduced { checkWrongCollision(particleMother, collision, indexCollisionMaxNumContrib, flagWrongCollision); } } - rowHfJpsiPhiMcRecReduced(indexHfCandJpsi, selectedTracksBach[0][vecDaughtersB.back().globalIndex()], selectedTracksBach[1][vecDaughtersB.back().globalIndex()], flag, channel, flagWrongCollision, debug, motherPt); + rowHfJpsiPhiMcRecReduced(indexHfCandJpsi, selectedTracksBach[0][vecDaughtersB[2].globalIndex()], selectedTracksBach[1][vecDaughtersB[3].globalIndex()], flag, channel, flagWrongCollision, debug, motherPt); } } @@ -594,6 +635,44 @@ struct HfDataCreatorJpsiHadReduced { rowHfBpMcGenReduced(flag, channel, ptParticle, yParticle, etaParticle, ptProngs[0], yProngs[0], etaProngs[0], ptProngs[1], yProngs[1], etaProngs[1], hfRejMap, centFT0C, centFT0M); + } else if constexpr (DecChannel == DecayChannel::B0ToJpsiK0Star) { + // B0 → J/Psi K*0 → (µ+µ-) (K+pi-) + if (RecoDecay::isMatchedMCGen(particlesMc, particle, Pdg::kB0, std::array{static_cast(Pdg::kJPsi), +kKPlus, -kPiPlus}, true, &sign, 2)) { + // Match J/Psi -> µ+µ- and K*0 -> K+pi- + auto candJpsiMC = particlesMc.rawIteratorAt(particle.daughtersIds().front()); + auto candK0StarMC = particlesMc.rawIteratorAt(particle.daughtersIds().back()); + // Printf("Checking J/Psi -> µ+µ- and K*0 -> K+pi-"); + if (RecoDecay::isMatchedMCGen(particlesMc, candJpsiMC, static_cast(Pdg::kJPsi), std::array{-kMuonMinus, +kMuonMinus}, true)) { + flag = sign * o2::hf_decay::hf_cand_beauty::B0ToJpsiPiK; + } + // Check K*0 -> K+pi- + if (RecoDecay::isMatchedMCGen(particlesMc, candK0StarMC, static_cast(Pdg::kK0Star892), std::array{+kKPlus, -kPiPlus}, true)) { + channel = o2::hf_decay::hf_cand_beauty::B0ToJpsiPiK; + } + } + + // save information for B0 task + if (std::abs(flag) != o2::hf_decay::hf_cand_beauty::B0ToJpsiPiK) { + continue; + } + + auto ptParticle = particle.pt(); + auto yParticle = RecoDecay::y(particle.pVector(), MassB0); + auto etaParticle = particle.eta(); + + std::array ptProngs{}; + std::array yProngs{}; + std::array etaProngs{}; + int counter = 0; + for (const auto& daught : particle.daughters_as()) { + ptProngs[counter] = daught.pt(); + etaProngs[counter] = daught.eta(); + yProngs[counter] = RecoDecay::y(daught.pVector(), pdg->Mass(daught.pdgCode())); + counter++; + } + rowHfB0McGenReduced(flag, channel, ptParticle, yParticle, etaParticle, + ptProngs[0], yProngs[0], etaProngs[0], + ptProngs[1], yProngs[1], etaProngs[1], hfRejMap, centFT0C, centFT0M); } else if constexpr (DecChannel == DecayChannel::BsToJpsiPhi) { // Bs → J/Psi phi → (µ+µ-) (K+K-) if (RecoDecay::isMatchedMCGen(particlesMc, particle, Pdg::kBS, std::array{static_cast(Pdg::kJPsi), +kKPlus, -kKPlus}, true, &sign, 2)) { @@ -603,7 +682,6 @@ struct HfDataCreatorJpsiHadReduced { // Printf("Checking J/Psi -> µ+µ- and phi -> K+K-"); if (RecoDecay::isMatchedMCGen(particlesMc, candJpsiMC, static_cast(Pdg::kJPsi), std::array{-kMuonMinus, +kMuonMinus}, true)) { flag = sign * o2::hf_decay::hf_cand_beauty::BsToJpsiKK; - flag = sign * o2::hf_decay::hf_cand_beauty::BsToJpsiKK; } // Check phi -> K+K- if (RecoDecay::isMatchedMCGen(particlesMc, candPhiMC, static_cast(Pdg::kPhi), std::array{-kKPlus, +kKPlus}, true)) { @@ -617,7 +695,7 @@ struct HfDataCreatorJpsiHadReduced { } auto ptParticle = particle.pt(); - auto yParticle = RecoDecay::y(particle.pVector(), MassBPlus); + auto yParticle = RecoDecay::y(particle.pVector(), MassBS); auto etaParticle = particle.eta(); std::array ptProngs{}; @@ -806,7 +884,7 @@ struct HfDataCreatorJpsiHadReduced { continue; } // compute invariant mass square and apply selection - invMass2JpsiHad = RecoDecay::m2(std::array{pVecJpsi, pVecBach}, std::array{MassJPsi, MassKPlus}); + invMass2JpsiHad = RecoDecay::m2(std::array{pVecJpsi, pVec2}, std::array{MassJPsi, MassKPlus}); if ((invMass2JpsiHad < invMass2JpsiHadMin) || (invMass2JpsiHad > invMass2JpsiHadMax)) { continue; } @@ -846,6 +924,133 @@ struct HfDataCreatorJpsiHadReduced { fillMcRecoInfo(collision, particlesMc, beautyHadDauTracks, indexHfCandJpsi, std::array, 2>{selectedTracksBach}, indexCollisionMaxNumContrib); } fillHfCandJpsi = true; + } else if constexpr (DecChannel == DecayChannel::B0ToJpsiK0Star) { + for (auto trackBachId2 = trackId + 1; trackBachId2 != trackIndices.end(); ++trackBachId2) { + auto trackBach2 = trackBachId2.template track_as(); + auto trackBach2ParCov = getTrackParCov(trackBach2); + + std::array dcaBach2{trackBach2.dcaXY(), trackBach2.dcaZ()}; + std::array pVecBach2 = trackBach2.pVector(); + if (trackBach2.collisionId() != thisCollId) { + o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, trackBach2ParCov, 2.f, noMatCorr, &dcaBach2); + getPxPyPz(trackBach2ParCov, pVecBach2); + } + + // apply selections on bachelor tracks + if (!isTrackSelected(trackBach2, trackBach2ParCov, dcaBach2, jPsiDauTracks)) { + continue; + } + std::array pVec2{trackBach.pVector()}, pVec3{trackBach2.pVector()}; + auto invMassPiK = RecoDecay::m(std::array{pVec2, pVec3}, std::array{MassPiPlus, MassKPlus}); + auto invMassKPi = RecoDecay::m(std::array{pVec2, pVec3}, std::array{MassKPlus, MassPiPlus}); + bool isK0StarPiK = std::abs(invMassPiK - MassK0Star892) < deltaMK0StarMax; + bool isK0StarKPi = std::abs(invMassKPi - MassK0Star892) < deltaMK0StarMax; + + if (!isK0StarPiK && !isK0StarKPi) { + continue; + } + + // --------------------------------- + // reconstruct B0 candidate secondary vertex + + registry.fill(HIST("hFitCandidatesB0"), SVFitting::BeforeFit); + try { + if (df4.process(trackPosParCov, trackNegParCov, trackParCovBach, trackBach2ParCov) == 0) { + continue; + } + } catch (const std::runtime_error& error) { + LOG(info) << "Run time error found: " << error.what() << ". DCAFitterN cannot work, skipping the candidate."; + registry.fill(HIST("hFitCandidatesB0"), SVFitting::Fail); + continue; + } + registry.fill(HIST("hFitCandidatesB0"), SVFitting::FitOk); + + o2::track::TrackParCov trackParCovB0{}; + std::array pVecB0{}, pVec0{}, pVec1{}, pVecK0Star{}; + + auto secondaryVertexB0 = df4.getPCACandidate(); + df4.getTrack(0).getPxPyPzGlo(pVec0); + df4.getTrack(1).getPxPyPzGlo(pVec1); + df4.getTrack(2).getPxPyPzGlo(pVec2); + df4.getTrack(3).getPxPyPzGlo(pVec3); + pVecB0 = RecoDecay::pVec(pVec0, pVec1, pVec2, pVec3); + pVecJpsi = RecoDecay::pVec(pVec0, pVec1); + pVecK0Star = RecoDecay::pVec(pVec2, pVec3); + trackParCovB0 = df4.createParentTrackParCov(); + trackParCovB0.setAbsCharge(0); // to be sure + + if (!isBSelected(pVecB0, secondaryVertexB0, collision)) { + continue; + } + + registry.fill(HIST("hPtK0Star"), RecoDecay::pt(pVecBach, pVecBach2)); + registry.fill(HIST("hMassK0Star"), RecoDecay::m(std::array{pVecBach, pVecBach2}, isK0StarPiK ? std::array{MassPiPlus, MassKPlus} : std::array{MassKPlus, MassPiPlus})); + invMass2JpsiHad = RecoDecay::m2(std::array{pVecJpsi, pVecK0Star}, std::array{MassJPsi, MassK0Star892}); + if ((invMass2JpsiHad < invMass2JpsiHadMin) || (invMass2JpsiHad > invMass2JpsiHadMax)) { + continue; + } + registry.fill(HIST("hMassJpsiK0Star"), std::sqrt(invMass2JpsiHad)); + + // fill daughter tracks table + // if information on track already stored, go to next track + if (!selectedTracksBach.count(trackBach.globalIndex())) { + hfTrackLfDau0(trackBach.globalIndex(), indexHfReducedCollision, + trackParCovBach.getX(), trackParCovBach.getAlpha(), + trackParCovBach.getY(), trackParCovBach.getZ(), trackParCovBach.getSnp(), + trackParCovBach.getTgl(), trackParCovBach.getQ2Pt(), + trackBach.itsNCls(), trackBach.tpcNClsCrossedRows(), trackBach.tpcChi2NCl(), trackBach.itsChi2NCl(), + trackBach.hasTPC(), trackBach.hasTOF(), + trackBach.tpcNSigmaPi(), trackBach.tofNSigmaPi(), + trackBach.tpcNSigmaKa(), trackBach.tofNSigmaKa(), + trackBach.tpcNSigmaPr(), trackBach.tofNSigmaPr()); + hfTrackCovLfDau0(trackParCovBach.getSigmaY2(), trackParCovBach.getSigmaZY(), trackParCovBach.getSigmaZ2(), + trackParCovBach.getSigmaSnpY(), trackParCovBach.getSigmaSnpZ(), + trackParCovBach.getSigmaSnp2(), trackParCovBach.getSigmaTglY(), trackParCovBach.getSigmaTglZ(), + trackParCovBach.getSigmaTglSnp(), trackParCovBach.getSigmaTgl2(), + trackParCovBach.getSigma1PtY(), trackParCovBach.getSigma1PtZ(), trackParCovBach.getSigma1PtSnp(), + trackParCovBach.getSigma1PtTgl(), trackParCovBach.getSigma1Pt2()); + // add trackBach.globalIndex() to a list + // to keep memory of the pions filled in the table and avoid refilling them if they are paired to another Jpsi candidate + // and keep track of their index in hfTrackLfDau0 for McRec purposes + selectedTracksBach[trackBach.globalIndex()] = hfTrackLfDau0.lastIndex(); + } + + // fill daughter tracks table + // if information on track already stored, go to next track + if (!selectedTracksBach2.count(trackBach2.globalIndex())) { + hfTrackLfDau1(trackBach2.globalIndex(), indexHfReducedCollision, + trackBach2ParCov.getX(), trackBach2ParCov.getAlpha(), + trackBach2ParCov.getY(), trackBach2ParCov.getZ(), trackBach2ParCov.getSnp(), + trackBach2ParCov.getTgl(), trackBach2ParCov.getQ2Pt(), + trackBach2.itsNCls(), trackBach2.tpcNClsCrossedRows(), trackBach2.tpcChi2NCl(), trackBach2.itsChi2NCl(), + trackBach2.hasTPC(), trackBach2.hasTOF(), + trackBach2.tpcNSigmaPi(), trackBach2.tofNSigmaPi(), + trackBach2.tpcNSigmaKa(), trackBach2.tofNSigmaKa(), + trackBach2.tpcNSigmaPr(), trackBach2.tofNSigmaPr()); + hfTrackCovLfDau1(trackBach2ParCov.getSigmaY2(), trackBach2ParCov.getSigmaZY(), trackBach2ParCov.getSigmaZ2(), + trackBach2ParCov.getSigmaSnpY(), trackBach2ParCov.getSigmaSnpZ(), + trackBach2ParCov.getSigmaSnp2(), trackBach2ParCov.getSigmaTglY(), trackBach2ParCov.getSigmaTglZ(), + trackBach2ParCov.getSigmaTglSnp(), trackBach2ParCov.getSigmaTgl2(), + trackBach2ParCov.getSigma1PtY(), trackBach2ParCov.getSigma1PtZ(), trackBach2ParCov.getSigma1PtSnp(), + trackBach2ParCov.getSigma1PtTgl(), trackBach2ParCov.getSigma1Pt2()); + // add trackBach2.globalIndex() to a list + // to keep memory of the pions filled in the table and avoid refilling them if they are paired to another Jpsi candidate + // and keep track of their index in hfTrackLfDau1 for McRec purposes + selectedTracksBach2[trackBach2.globalIndex()] = hfTrackLfDau1.lastIndex(); + } + + if constexpr (DoMc) { + std::vector beautyHadDauTracks{}; + beautyHadDauTracks.reserve(jPsiDauTracks.size()); + for (const auto& track : jPsiDauTracks) { + beautyHadDauTracks.push_back(track); + } + beautyHadDauTracks.push_back(trackBach); + beautyHadDauTracks.push_back(trackBach2); + fillMcRecoInfo(collision, particlesMc, beautyHadDauTracks, indexHfCandJpsi, std::array, 2>{selectedTracksBach, selectedTracksBach2}, indexCollisionMaxNumContrib); + } + fillHfCandJpsi = true; + } } else if constexpr (DecChannel == DecayChannel::BsToJpsiPhi) { for (auto trackBachId2 = trackId + 1; trackBachId2 != trackIndices.end(); ++trackBachId2) { auto trackBach2 = trackBachId2.template track_as(); @@ -965,7 +1170,8 @@ struct HfDataCreatorJpsiHadReduced { beautyHadDauTracks.push_back(track); } beautyHadDauTracks.push_back(trackBach); - fillMcRecoInfo(collision, particlesMc, beautyHadDauTracks, indexHfCandJpsi, std::array, 2>{selectedTracksBach, selectedTracksBach2}, indexCollisionMaxNumContrib); + beautyHadDauTracks.push_back(trackBach2); + fillMcRecoInfo(collision, particlesMc, beautyHadDauTracks, indexHfCandJpsi, std::array, 2>{selectedTracksBach, selectedTracksBach2}, indexCollisionMaxNumContrib); } fillHfCandJpsi = true; } @@ -1060,7 +1266,7 @@ struct HfDataCreatorJpsiHadReduced { TracksPidWithSel const& tracks, aod::BCsWithTimestamps const& bcs) { - // store configurables needed for B0 workflow + // store configurables needed for Bs workflow if (!isHfCandBhadConfigFilled) { rowCandidateConfigBs(invMassWindowJpsiHad.value); isHfCandBhadConfigFilled = true; @@ -1082,6 +1288,34 @@ struct HfDataCreatorJpsiHadReduced { } PROCESS_SWITCH(HfDataCreatorJpsiHadReduced, processJpsiPhiData, "Process J/Psi phi without MC info", false); + void processJpsiK0StarData(soa::Join const& collisions, + aod::HfCand2ProngWPid const& candsJpsi, + aod::TrackAssoc const& trackIndices, + TracksPidWithSel const& tracks, + aod::BCsWithTimestamps const& bcs) + { + // store configurables needed for B0 workflow + if (!isHfCandBhadConfigFilled) { + rowCandidateConfigB0(invMassWindowJpsiHad.value); + isHfCandBhadConfigFilled = true; + } + + int zvtxColl{0}; + int sel8Coll{0}; + int zvtxAndSel8Coll{0}; + int zvtxAndSel8CollAndSoftTrig{0}; + int allSelColl{0}; + for (const auto& collision : collisions) { + auto thisCollId = collision.globalIndex(); + auto candsJpsiThisColl = candsJpsi.sliceBy(candsJpsiPerCollision, thisCollId); + auto trackIdsThisCollision = trackIndices.sliceBy(trackIndicesPerCollision, thisCollId); + runDataCreation(collision, candsJpsiThisColl, trackIdsThisCollision, tracks, tracks, -1, bcs, zvtxColl, sel8Coll, zvtxAndSel8Coll, zvtxAndSel8CollAndSoftTrig, allSelColl); + } + // handle normalization by the right number of collisions + hfCollisionCounter(collisions.tableSize(), zvtxColl, sel8Coll, zvtxAndSel8Coll, zvtxAndSel8CollAndSoftTrig, allSelColl); + } + PROCESS_SWITCH(HfDataCreatorJpsiHadReduced, processJpsiK0StarData, "Process J/Psi K*0 without MC info", false); + void processJpsiKMc(CollisionsWCMcLabels const& collisions, aod::HfCand2ProngWPid const& candsJpsi, aod::TrackAssoc const& trackIndices, @@ -1117,6 +1351,41 @@ struct HfDataCreatorJpsiHadReduced { } PROCESS_SWITCH(HfDataCreatorJpsiHadReduced, processJpsiKMc, "Process J/Psi K with MC info", false); + void processJpsiK0StarMc(CollisionsWCMcLabels const& collisions, + aod::HfCand2ProngWPid const& candsJpsi, + aod::TrackAssoc const& trackIndices, + TracksPidWithSelAndMc const& tracks, + aod::McParticles const& particlesMc, + BCsInfo const& bcs, + McCollisions const& mcCollisions) + { + // store configurables needed for B+ workflow + if (!isHfCandBhadConfigFilled) { + rowCandidateConfigB0(invMassWindowJpsiHad.value); + isHfCandBhadConfigFilled = true; + } + + int zvtxColl{0}; + int sel8Coll{0}; + int zvtxAndSel8Coll{0}; + int zvtxAndSel8CollAndSoftTrig{0}; + int allSelColl{0}; + for (const auto& collision : collisions) { + auto thisCollId = collision.globalIndex(); + auto candsJpsiThisColl = candsJpsi.sliceBy(candsJpsiPerCollision, thisCollId); + auto trackIdsThisCollision = trackIndices.sliceBy(trackIndicesPerCollision, thisCollId); + auto collsSameMcCollision = collisions.sliceBy(colPerMcCollision, collision.mcCollisionId()); + int64_t const indexCollisionMaxNumContrib = getIndexCollisionMaxNumContrib(collsSameMcCollision); + runDataCreation(collision, candsJpsiThisColl, trackIdsThisCollision, tracks, particlesMc, indexCollisionMaxNumContrib, bcs, zvtxColl, sel8Coll, zvtxAndSel8Coll, zvtxAndSel8CollAndSoftTrig, allSelColl); + } + // handle normalization by the right number of collisions + hfCollisionCounter(collisions.tableSize(), zvtxColl, sel8Coll, zvtxAndSel8Coll, zvtxAndSel8CollAndSoftTrig, allSelColl); + for (const auto& mcCollision : mcCollisions) { + runMcGen(mcCollision, particlesMc, collisions, bcs); + } + } + PROCESS_SWITCH(HfDataCreatorJpsiHadReduced, processJpsiK0StarMc, "Process J/Psi K0* with MC info", false); + void processJpsiPhiMc(CollisionsWCMcLabels const& collisions, aod::HfCand2ProngWPid const& candsJpsi, aod::TrackAssoc const& trackIndices, diff --git a/PWGHF/DataModel/CandidateReconstructionTables.h b/PWGHF/DataModel/CandidateReconstructionTables.h index 46931f96988..04f0cdde0cd 100644 --- a/PWGHF/DataModel/CandidateReconstructionTables.h +++ b/PWGHF/DataModel/CandidateReconstructionTables.h @@ -1860,7 +1860,17 @@ DECLARE_SOA_TABLE(HfCandLbMcGen, "AOD", "HFCANDLBMCGEN", //! // specific B0 candidate properties namespace hf_cand_b0 { -DECLARE_SOA_INDEX_COLUMN_FULL(Prong0, prong0, int, HfCand3Prong, "_0"); // D index +DECLARE_SOA_INDEX_COLUMN_FULL(Prong0, prong0, int, HfCand3Prong, "_0"); // D index +DECLARE_SOA_DYNAMIC_COLUMN(ImpactParameterProduct, impactParameterProduct, // Impact parameter product for B0 -> J/Psi K*0 + [](float pxJpsiDauPos, float pyJpsiDauPos, float pzJpsiDauPos, float pxJpsiDauNeg, float pyJpsiDauNeg, float pzJpsiDauNeg, float pxLfTrack0, float pyLfTrack0, float pzLfTrack0, float pxLfTrack1, float pyLfTrack1, float pzLfTrack1, float xVtxP, float yVtxP, float zVtxP, float xVtxS, float yVtxS, float zVtxS) -> float { + float impParJpsi = RecoDecay::impParXY(std::array{xVtxP, yVtxP, zVtxP}, std::array{xVtxS, yVtxS, zVtxS}, RecoDecay::pVec(std::array{pxJpsiDauPos, pyJpsiDauPos, pzJpsiDauPos}, std::array{pxJpsiDauNeg, pyJpsiDauNeg, pzJpsiDauNeg})); + float impParK0Star = RecoDecay::impParXY(std::array{xVtxP, yVtxP, zVtxP}, std::array{xVtxS, yVtxS, zVtxS}, RecoDecay::pVec(std::array{pxLfTrack0, pyLfTrack0, pzLfTrack0}, std::array{pxLfTrack1, pyLfTrack1, pzLfTrack1})); + return impParJpsi * impParK0Star; + }); +DECLARE_SOA_DYNAMIC_COLUMN(ImpactParameterProductJpsi, impactParameterProductJpsi, // J/Psi impact parameter for Bs -> J/Psi phi + [](float dcaDauPos, float dcaDauNeg) -> float { return dcaDauPos * dcaDauNeg; }); +DECLARE_SOA_DYNAMIC_COLUMN(ImpactParameterProductK0Star, impactParameterProductK0Star, // K*0 impact parameter for B0 -> J/Psi K*0 + [](float dcaLfTrack0, float dcaLfTrack1) -> float { return dcaLfTrack0 * dcaLfTrack1; }); enum DecayTypeMc : uint8_t { B0ToDplusPiToPiKPiPi = 0, B0ToDsPiToKKPiPi, @@ -2070,6 +2080,61 @@ DECLARE_SOA_DYNAMIC_COLUMN(CtXY, ctXY, //! [](float px0, float py0, float pz0, float px1, float py1, float pz1, float px2, float py2, float pz2, float px3, float py3, float pz3, float xVtxP, float yVtxP, float xVtxS, float yVtxS, const std::array& m) -> float { return RecoDecay::ctXY(std::array{xVtxP, yVtxP}, std::array{xVtxS, yVtxS}, std::array{std::array{px0, py0, pz0}, std::array{px1, py1, pz1}, std::array{px2, py2, pz2}, std::array{px3, py3, pz3}}, m); }); } // namespace hf_cand_4prong +// declare dedicated B0 -> J/Psi K*0 decay candidate table +// convention: prongs 0 and 1 should be J/Psi decay products, 2 and 3 should be K*0 decay products +DECLARE_SOA_TABLE(HfCandB0JPBase, "AOD", "HFCANDB0JPBASE", + // general columns + HFCAND_COLUMNS, + /* prong 2 */ hf_cand::ImpactParameterNormalised2, + hf_cand::PtProng2, + hf_cand::Pt2Prong2, + hf_cand::PVectorProng2, + /* prong 3 */ hf_cand::ImpactParameterNormalised3, + hf_cand::PtProng3, + hf_cand::Pt2Prong3, + hf_cand::PVectorProng3, + // 4-prong specific columns + o2::soa::Index<>, + hf_cand::PxProng0, hf_cand::PyProng0, hf_cand::PzProng0, + hf_cand::PxProng1, hf_cand::PyProng1, hf_cand::PzProng1, + hf_cand::PxProng2, hf_cand::PyProng2, hf_cand::PzProng2, + hf_cand::PxProng3, hf_cand::PyProng3, hf_cand::PzProng3, + hf_cand::ImpactParameter0, hf_cand::ImpactParameter1, hf_cand::ImpactParameter2, hf_cand::ImpactParameter3, + hf_cand::ErrorImpactParameter0, hf_cand::ErrorImpactParameter1, hf_cand::ErrorImpactParameter2, hf_cand::ErrorImpactParameter3, + /* dynamic columns */ + hf_cand_4prong::M, + hf_cand_4prong::M2, + hf_cand_4prong::ImpactParameterProngSqSum, + hf_cand_b0::ImpactParameterProduct, + hf_cand_b0::ImpactParameterProductJpsi, + hf_cand_b0::ImpactParameterProductK0Star, + /* dynamic columns that use candidate momentum components */ + hf_cand::Pt, + hf_cand::Pt2, + hf_cand::P, + hf_cand::P2, + hf_cand::PVector, + hf_cand::Cpa, + hf_cand::CpaXY, + hf_cand::Ct, + hf_cand::ImpactParameterXY, + hf_cand_4prong::MaxNormalisedDeltaIP, + hf_cand::Eta, + hf_cand::Phi, + hf_cand::Y, + hf_cand::E, + hf_cand::E2, + hf_cand_4prong::CtXY); + +// extended table with expression columns that can be used as arguments of dynamic columns +DECLARE_SOA_EXTENDED_TABLE_USER(HfCandB0JPExt, HfCandB0JPBase, "HFCANDB0JPEXT", + hf_cand_4prong::Px, hf_cand_4prong::Py, hf_cand_4prong::Pz); + +DECLARE_SOA_TABLE(HfCandB0JPDaus, "AOD", "HFCANDB0JPDAUS", + hf_cand_b0::Prong0Id, hf_track_index::Prong1Id, hf_track_index::Prong2Id, hf_track_index::Prong3Id); + +using HfCandB0ToJpsi = soa::Join; + // declare dedicated Bs -> J/Psi phi decay candidate table // convention: prongs 0 and 1 should be J/Psi decay products, 2 and 3 should be phi decay products DECLARE_SOA_TABLE(HfCandBsJPBase, "AOD", "HFCANDBSJPBASE",