From 5b8d1b62978d027eec6335a282897ab52a78df4a Mon Sep 17 00:00:00 2001 From: Fabrizio Grosa Date: Thu, 19 Oct 2023 10:28:44 +0200 Subject: [PATCH 1/2] Add the possibility to store D-meson ML scores in reduced B data model --- PWGHF/D2H/DataModel/ReducedDataModel.h | 20 +- .../TableProducer/dataCreatorD0PiReduced.cxx | 479 ++++++++++------- .../dataCreatorDplusPiReduced.cxx | 499 ++++++++++-------- 3 files changed, 580 insertions(+), 418 deletions(-) diff --git a/PWGHF/D2H/DataModel/ReducedDataModel.h b/PWGHF/D2H/DataModel/ReducedDataModel.h index 7a0279e2417..35a58606074 100644 --- a/PWGHF/D2H/DataModel/ReducedDataModel.h +++ b/PWGHF/D2H/DataModel/ReducedDataModel.h @@ -147,9 +147,12 @@ using HfRedTracks = HfRedTracksExt; namespace hf_charm_cand_reduced { -DECLARE_SOA_COLUMN(InvMass, invMass, float); //! Invariant mass of 2prong candidate in GeV/c2 -DECLARE_SOA_COLUMN(InvMassD0, invMassD0, float); //! Invariant mass of 2prong candidate in GeV/c2 -DECLARE_SOA_COLUMN(InvMassD0Bar, invMassD0Bar, float); //! Invariant mass of 2prong candidate in GeV/c2 +DECLARE_SOA_COLUMN(InvMass, invMass, float); //! Invariant mass of 2prong candidate in GeV/c2 +DECLARE_SOA_COLUMN(InvMassD0, invMassD0, float); //! Invariant mass of 2prong candidate in GeV/c2 +DECLARE_SOA_COLUMN(InvMassD0Bar, invMassD0Bar, float); //! Invariant mass of 2prong candidate in GeV/c2 +DECLARE_SOA_COLUMN(MlScoreBkg, mlScoreBkg, float); //! ML score for background class +DECLARE_SOA_COLUMN(MlScorePrompt, mlScorePrompt, float); //! ML score for prompt class +DECLARE_SOA_COLUMN(MlScoreNonprompt, mlScoreNonprompt, float); //! ML score for non-prompt class } // namespace hf_charm_cand_reduced // CAREFUL: need to follow convention [Name = Description + 's'] in DECLARE_SOA_TABLE(Name, "AOD", Description) @@ -170,6 +173,11 @@ DECLARE_SOA_TABLE(HfRed2ProngsCov, "AOD", "HFRED2PRONGSCOV", //! Table with 2pro HFTRACKPARCOV_COLUMNS, o2::soa::Marker<1>); +DECLARE_SOA_TABLE(HfRed2ProngsMl, "AOD", "HFRED2PRONGML", //! Table with 2prong candidate ML scores + hf_charm_cand_reduced::MlScoreBkg, + hf_charm_cand_reduced::MlScorePrompt, + hf_charm_cand_reduced::MlScoreNonprompt); + // CAREFUL: need to follow convention [Name = Description + 's'] in DECLARE_SOA_TABLE(Name, "AOD", Description) // to call DECLARE_SOA_INDEX_COLUMN_FULL later on DECLARE_SOA_TABLE(HfRed3Prongs, "AOD", "HFRED3PRONG", //! Table with 3prong candidate information for reduced workflow @@ -188,6 +196,12 @@ DECLARE_SOA_TABLE(HfRed3ProngsCov, "AOD", "HFRED3PRONGSCOV", //! Table with 3pro HFTRACKPARCOV_COLUMNS, o2::soa::Marker<2>); +DECLARE_SOA_TABLE(HfRed3ProngsMl, "AOD", "HFRED3PRONGML", //! Table with 3prong candidate ML scores + hf_charm_cand_reduced::MlScoreBkg, + hf_charm_cand_reduced::MlScorePrompt, + hf_charm_cand_reduced::MlScoreNonprompt, + o2::soa::Marker<1>); + // Beauty candidates prongs namespace hf_cand_b0_reduced { diff --git a/PWGHF/D2H/TableProducer/dataCreatorD0PiReduced.cxx b/PWGHF/D2H/TableProducer/dataCreatorD0PiReduced.cxx index 17b480699af..cbf86f30fb6 100644 --- a/PWGHF/D2H/TableProducer/dataCreatorD0PiReduced.cxx +++ b/PWGHF/D2H/TableProducer/dataCreatorD0PiReduced.cxx @@ -13,6 +13,7 @@ /// \brief Creation of D0-Pi pairs /// /// \author Antonio Palasciano , Università degli Studi di Bari +/// \author Fabrizio Grosa , CERN #include @@ -56,6 +57,7 @@ struct HfDataCreatorD0PiReduced { Produces hfTrackPidPion; Produces hfCand2Prong; Produces hfCand2ProngCov; + Produces hfCand2ProngMl; Produces rowCandidateConfig; Produces rowHfD0PiMcRecReduced; Produces rowHfBpMcGenReduced; @@ -113,10 +115,12 @@ struct HfDataCreatorD0PiReduced { using TracksPIDWithSel = soa::Join; using TracksPIDWithSelAndMc = soa::Join; using CandsDFiltered = soa::Filtered>; + using CandsDFilteredWithMl = soa::Filtered>; Filter filterSelectCandidates = (aod::hf_sel_candidate_d0::isSelD0 >= selectionFlagD0 || aod::hf_sel_candidate_d0::isSelD0bar >= selectionFlagD0bar); Preslice candsDPerCollision = aod::track_association::collisionId; + Preslice candsDPerCollisionWithMl = aod::track_association::collisionId; Preslice trackIndicesPerCollision = aod::track_association::collisionId; HistogramRegistry registry{"registry"}; @@ -214,234 +218,224 @@ struct HfDataCreatorD0PiReduced { return true; } - template - void runDataCreation(aod::Collisions const& collisions, - CandsDFiltered const& candsD0, + template + void runDataCreation(aod::Collision const& collision, + C const& candsD0, aod::TrackAssoc const& trackIndices, T const& tracks, P const& particlesMc, aod::BCsWithTimestamps const& bcs) { - // store configurables needed for B+ workflow - if (!isHfCandBplusConfigFilled) { - rowCandidateConfig(selectionFlagD0.value, selectionFlagD0bar.value, invMassWindowD0Pi.value); - isHfCandBplusConfigFilled = true; + // helpers for ReducedTables filling + int indexHfReducedCollision = hfReducedCollision.lastIndex() + 1; + // std::map where the key is the track.globalIndex() and + // the value is the track index in the table of the selected pions + std::map selectedTracksPion; + bool fillHfReducedCollision = false; + + auto primaryVertex = getPrimaryVertex(collision); + + // Set the magnetic field from ccdb. + // The static instance of the propagator was already modified in the HFTrackIndexSkimCreator, + // but this is not true when running on Run2 data/MC already converted into AO2Ds. + auto bc = collision.bc_as(); + if (runNumber != bc.runNumber()) { + LOG(info) << ">>>>>>>>>>>> Current run number: " << runNumber; + initCCDB(bc, runNumber, ccdb, isRun2 ? ccdbPathGrp : ccdbPathGrpMag, lut, isRun2); + bz = o2::base::Propagator::Instance()->getNominalBz(); + LOG(info) << ">>>>>>>>>>>> Magnetic field: " << bz; } + df2.setBz(bz); - // handle normalization by the right number of collisions - hfCollisionCounter(collisions.tableSize()); - - for (const auto& collision : collisions) { + auto thisCollId = collision.globalIndex(); + for (const auto& candD0 : candsD0) { + int indexHfCand2Prong = hfCand2Prong.lastIndex() + 1; + bool fillHfCand2Prong = false; + float invMassD0{-1.f}, invMassD0bar{-1.f}; - // helpers for ReducedTables filling - int indexHfReducedCollision = hfReducedCollision.lastIndex() + 1; - // std::map where the key is the track.globalIndex() and - // the value is the track index in the table of the selected pions - std::map selectedTracksPion; - bool fillHfReducedCollision = false; - - auto primaryVertex = getPrimaryVertex(collision); - - // Set the magnetic field from ccdb. - // The static instance of the propagator was already modified in the HFTrackIndexSkimCreator, - // but this is not true when running on Run2 data/MC already converted into AO2Ds. - auto bc = collision.bc_as(); - if (runNumber != bc.runNumber()) { - LOG(info) << ">>>>>>>>>>>> Current run number: " << runNumber; - initCCDB(bc, runNumber, ccdb, isRun2 ? ccdbPathGrp : ccdbPathGrpMag, lut, isRun2); - bz = o2::base::Propagator::Instance()->getNominalBz(); - LOG(info) << ">>>>>>>>>>>> Magnetic field: " << bz; + if (candD0.isSelD0() >= selectionFlagD0) { + invMassD0 = hfHelper.invMassD0ToPiK(candD0); + registry.fill(HIST("hMassD0ToKPi"), invMassD0); } - df2.setBz(bz); - - auto thisCollId = collision.globalIndex(); - auto candsDThisColl = candsD0.sliceBy(candsDPerCollision, thisCollId); - for (const auto& candD0 : candsDThisColl) { - int indexHfCand2Prong = hfCand2Prong.lastIndex() + 1; - bool fillHfCand2Prong = false; - float invMassD0{-1.f}, invMassD0bar{-1.f}; - - if (candD0.isSelD0() >= selectionFlagD0) { - invMassD0 = hfHelper.invMassD0ToPiK(candD0); - registry.fill(HIST("hMassD0ToKPi"), invMassD0); - } - if (candD0.isSelD0bar() >= selectionFlagD0bar) { - invMassD0bar = hfHelper.invMassD0barToKPi(candD0); - registry.fill(HIST("hMassD0ToKPi"), invMassD0bar); - } - registry.fill(HIST("hPtD0"), candD0.pt()); - registry.fill(HIST("hCPAD0"), candD0.cpa()); + if (candD0.isSelD0bar() >= selectionFlagD0bar) { + invMassD0bar = hfHelper.invMassD0barToKPi(candD0); + registry.fill(HIST("hMassD0ToKPi"), invMassD0bar); + } + registry.fill(HIST("hPtD0"), candD0.pt()); + registry.fill(HIST("hCPAD0"), candD0.cpa()); - auto track0 = candD0.template prong0_as(); - auto track1 = candD0.template prong1_as(); - auto trackParCov0 = getTrackParCov(track0); - auto trackParCov1 = getTrackParCov(track1); + auto track0 = candD0.template prong0_as(); + auto track1 = candD0.template prong1_as(); + auto trackParCov0 = getTrackParCov(track0); + auto trackParCov1 = getTrackParCov(track1); - std::array pVec0 = {track0.px(), track0.py(), track0.pz()}; - std::array pVec1 = {track1.px(), track1.py(), track1.pz()}; + std::array pVec0 = {track0.px(), track0.py(), track0.pz()}; + std::array pVec1 = {track1.px(), track1.py(), track1.pz()}; - auto dca0 = o2::dataformats::DCA(track0.dcaXY(), track0.dcaZ(), track0.cYY(), track0.cZY(), track0.cZZ()); - auto dca1 = o2::dataformats::DCA(track1.dcaXY(), track1.dcaZ(), track1.cYY(), track1.cZY(), track1.cZZ()); + auto dca0 = o2::dataformats::DCA(track0.dcaXY(), track0.dcaZ(), track0.cYY(), track0.cZY(), track0.cZZ()); + auto dca1 = o2::dataformats::DCA(track1.dcaXY(), track1.dcaZ(), track1.cYY(), track1.cZY(), track1.cZZ()); - // repropagate tracks to this collision if needed - if (track0.collisionId() != thisCollId) { - trackParCov0.propagateToDCA(primaryVertex, bz, &dca0); - } + // repropagate tracks to this collision if needed + if (track0.collisionId() != thisCollId) { + trackParCov0.propagateToDCA(primaryVertex, bz, &dca0); + } - if (track1.collisionId() != thisCollId) { - trackParCov1.propagateToDCA(primaryVertex, bz, &dca1); - } + if (track1.collisionId() != thisCollId) { + trackParCov1.propagateToDCA(primaryVertex, bz, &dca1); + } - // --------------------------------- - // reconstruct 2-prong secondary vertex (D0) - if (df2.process(trackParCov0, trackParCov1) == 0) { - continue; - } + // --------------------------------- + // reconstruct 2-prong secondary vertex (D0) + if (df2.process(trackParCov0, trackParCov1) == 0) { + continue; + } - const auto& secondaryVertexD0 = df2.getPCACandidate(); - // propagate the 2 prongs to the secondary vertex - trackParCov0.propagateTo(secondaryVertexD0[0], bz); - trackParCov1.propagateTo(secondaryVertexD0[0], bz); + const auto& secondaryVertexD0 = df2.getPCACandidate(); + // propagate the 2 prongs to the secondary vertex + trackParCov0.propagateTo(secondaryVertexD0[0], bz); + trackParCov1.propagateTo(secondaryVertexD0[0], bz); - // update pVec of tracks - df2.getTrack(0).getPxPyPzGlo(pVec0); - df2.getTrack(1).getPxPyPzGlo(pVec1); + // update pVec of tracks + df2.getTrack(0).getPxPyPzGlo(pVec0); + df2.getTrack(1).getPxPyPzGlo(pVec1); - // D0(bar) → π∓ K± - std::array pVecD0 = RecoDecay::pVec(pVec0, pVec1); - auto trackParCovD0 = o2::dataformats::V0(df2.getPCACandidatePos(), pVecD0, df2.calcPCACovMatrixFlat(), trackParCov0, trackParCov1); + // D0(bar) → π∓ K± + std::array pVecD0 = RecoDecay::pVec(pVec0, pVec1); + auto trackParCovD0 = o2::dataformats::V0(df2.getPCACandidatePos(), pVecD0, df2.calcPCACovMatrixFlat(), trackParCov0, trackParCov1); - auto trackIdsThisCollision = trackIndices.sliceBy(trackIndicesPerCollision, thisCollId); - for (const auto& trackId : trackIdsThisCollision) { - auto trackPion = trackId.template track_as(); + for (const auto& trackId : trackIndices) { + auto trackPion = trackId.template track_as(); - // apply selections on pion tracks - if (!isPionSelected(trackPion, track0, track1, candD0) || !isSelectedTrackDCA(trackPion)) { - continue; - } - registry.fill(HIST("hPtPion"), trackPion.pt()); - std::array pVecPion = {trackPion.px(), trackPion.py(), trackPion.pz()}; - // compute invariant mass square and apply selection - auto invMass2D0Pi = RecoDecay::m2(std::array{pVecD0, pVecPion}, std::array{massD0, massPi}); - if ((invMass2D0Pi < invMass2D0PiMin) || (invMass2D0Pi > invMass2D0PiMax)) { - continue; - } + // apply selections on pion tracks + if (!isPionSelected(trackPion, track0, track1, candD0) || !isSelectedTrackDCA(trackPion)) { + continue; + } + registry.fill(HIST("hPtPion"), trackPion.pt()); + std::array pVecPion = {trackPion.px(), trackPion.py(), trackPion.pz()}; + // compute invariant mass square and apply selection + auto invMass2D0Pi = RecoDecay::m2(std::array{pVecD0, pVecPion}, std::array{massD0, massPi}); + if ((invMass2D0Pi < invMass2D0PiMin) || (invMass2D0Pi > invMass2D0PiMax)) { + continue; + } - // fill Pion tracks table - // if information on track already stored, go to next track - if (!selectedTracksPion.count(trackPion.globalIndex())) { - hfTrackPion(trackPion.globalIndex(), indexHfReducedCollision, - trackPion.x(), trackPion.alpha(), - trackPion.y(), trackPion.z(), trackPion.snp(), - trackPion.tgl(), trackPion.signed1Pt()); - hfTrackCovPion(trackPion.cYY(), trackPion.cZY(), trackPion.cZZ(), - trackPion.cSnpY(), trackPion.cSnpZ(), - trackPion.cSnpSnp(), trackPion.cTglY(), trackPion.cTglZ(), - trackPion.cTglSnp(), trackPion.cTglTgl(), - trackPion.c1PtY(), trackPion.c1PtZ(), trackPion.c1PtSnp(), - trackPion.c1PtTgl(), trackPion.c1Pt21Pt2()); - hfTrackPidPion(trackPion.hasTPC(), trackPion.hasTOF(), - trackPion.tpcNSigmaPi(), trackPion.tofNSigmaPi()); - // add trackPion.globalIndex() to a list - // to keep memory of the pions filled in the table and avoid refilling them if they are paired to another D candidate - // and keep track of their index in hfTrackPion for McRec purposes - selectedTracksPion[trackPion.globalIndex()] = hfTrackPion.lastIndex(); - } + // fill Pion tracks table + // if information on track already stored, go to next track + if (!selectedTracksPion.count(trackPion.globalIndex())) { + hfTrackPion(trackPion.globalIndex(), indexHfReducedCollision, + trackPion.x(), trackPion.alpha(), + trackPion.y(), trackPion.z(), trackPion.snp(), + trackPion.tgl(), trackPion.signed1Pt()); + hfTrackCovPion(trackPion.cYY(), trackPion.cZY(), trackPion.cZZ(), + trackPion.cSnpY(), trackPion.cSnpZ(), + trackPion.cSnpSnp(), trackPion.cTglY(), trackPion.cTglZ(), + trackPion.cTglSnp(), trackPion.cTglTgl(), + trackPion.c1PtY(), trackPion.c1PtZ(), trackPion.c1PtSnp(), + trackPion.c1PtTgl(), trackPion.c1Pt21Pt2()); + hfTrackPidPion(trackPion.hasTPC(), trackPion.hasTOF(), + trackPion.tpcNSigmaPi(), trackPion.tofNSigmaPi()); + // add trackPion.globalIndex() to a list + // to keep memory of the pions filled in the table and avoid refilling them if they are paired to another D candidate + // and keep track of their index in hfTrackPion for McRec purposes + selectedTracksPion[trackPion.globalIndex()] = hfTrackPion.lastIndex(); + } - if constexpr (doMc) { - // we check the MC matching to be stored - auto arrayDaughtersD0 = std::array{track0, track1}; - auto arrayDaughtersBplus = std::array{track0, track1, trackPion}; - int8_t sign{0}; - int8_t flag{0}; - // B+ → D0(bar) π+ → (K+ π-) π+ - // Printf("Checking B+ → D0bar π+"); - auto indexRec = RecoDecay::getMatchedMCRec(particlesMc, arrayDaughtersBplus, pdg::Code::kBPlus, std::array{+kPiPlus, +kKPlus, -kPiPlus}, true, &sign, 2); + if constexpr (doMc) { + // we check the MC matching to be stored + auto arrayDaughtersD0 = std::array{track0, track1}; + auto arrayDaughtersBplus = std::array{track0, track1, trackPion}; + int8_t sign{0}; + int8_t flag{0}; + // B+ → D0(bar) π+ → (K+ π-) π+ + // Printf("Checking B+ → D0bar π+"); + auto indexRec = RecoDecay::getMatchedMCRec(particlesMc, arrayDaughtersBplus, pdg::Code::kBPlus, std::array{+kPiPlus, +kKPlus, -kPiPlus}, true, &sign, 2); + if (indexRec > -1) { + // D0bar → K+ π- + // Printf("Checking D0bar → K+ π-"); + indexRec = RecoDecay::getMatchedMCRec(particlesMc, arrayDaughtersD0, pdg::Code::kD0, std::array{+kPiPlus, -kKPlus}, true, &sign, 1); if (indexRec > -1) { - // D0bar → K+ π- - // Printf("Checking D0bar → K+ π-"); - indexRec = RecoDecay::getMatchedMCRec(particlesMc, arrayDaughtersD0, pdg::Code::kD0, std::array{+kPiPlus, -kKPlus}, true, &sign, 1); - if (indexRec > -1) { - flag = sign * BIT(hf_cand_bplus::DecayType::BplusToD0Pi); - } else { - LOGF(info, "WARNING: B+ decays in the expected final state but the condition on the intermediate state is not fulfilled"); - } + flag = sign * BIT(hf_cand_bplus::DecayType::BplusToD0Pi); + } else { + LOGF(info, "WARNING: B+ decays in the expected final state but the condition on the intermediate state is not fulfilled"); } - auto indexMother = RecoDecay::getMother(particlesMc, trackPion.template mcParticle_as

(), pdg::Code::kBPlus, true); - auto particleMother = particlesMc.rawIteratorAt(indexMother); - - rowHfD0PiMcRecReduced(indexHfCand2Prong, selectedTracksPion[trackPion.globalIndex()], flag, particleMother.pt()); } - fillHfCand2Prong = true; - } // pion loop - if (fillHfCand2Prong) { // fill candD0 table only once per D0 candidate - hfCand2Prong(track0.globalIndex(), track1.globalIndex(), - indexHfReducedCollision, - trackParCovD0.getX(), trackParCovD0.getAlpha(), - trackParCovD0.getY(), trackParCovD0.getZ(), trackParCovD0.getSnp(), - trackParCovD0.getTgl(), trackParCovD0.getQ2Pt(), - candD0.xSecondaryVertex(), candD0.ySecondaryVertex(), candD0.zSecondaryVertex(), invMassD0, invMassD0bar); - hfCand2ProngCov(trackParCovD0.getSigmaY2(), trackParCovD0.getSigmaZY(), trackParCovD0.getSigmaZ2(), - trackParCovD0.getSigmaSnpY(), trackParCovD0.getSigmaSnpZ(), - trackParCovD0.getSigmaSnp2(), trackParCovD0.getSigmaTglY(), trackParCovD0.getSigmaTglZ(), - trackParCovD0.getSigmaTglSnp(), trackParCovD0.getSigmaTgl2(), - trackParCovD0.getSigma1PtY(), trackParCovD0.getSigma1PtZ(), trackParCovD0.getSigma1PtSnp(), - trackParCovD0.getSigma1PtTgl(), trackParCovD0.getSigma1Pt2()); - fillHfReducedCollision = true; + auto indexMother = RecoDecay::getMother(particlesMc, trackPion.template mcParticle_as

(), pdg::Code::kBPlus, true); + auto particleMother = particlesMc.rawIteratorAt(indexMother); + + rowHfD0PiMcRecReduced(indexHfCand2Prong, selectedTracksPion[trackPion.globalIndex()], flag, particleMother.pt()); } - } // candsD loop - registry.fill(HIST("hEvents"), 1 + Event::Processed); - if (!fillHfReducedCollision) { - registry.fill(HIST("hEvents"), 1 + Event::NoD0PiSelected); - continue; - } - registry.fill(HIST("hEvents"), 1 + Event::D0PiSelected); - // fill collision table if it contains a D0Pi pair a minima - hfReducedCollision(collision.posX(), collision.posY(), collision.posZ(), - collision.covXX(), collision.covXY(), collision.covYY(), - collision.covXZ(), collision.covYZ(), collision.covZZ(), - bz); - } // collision - - if constexpr (doMc) { - // Match generated particles. - for (const auto& particle : particlesMc) { - int8_t sign{0}; - int8_t flag{0}; - // B+ → D0bar π+ - if (RecoDecay::isMatchedMCGen(particlesMc, particle, pdg::Code::kBPlus, std::array{static_cast(pdg::Code::kD0), +kPiPlus}, true)) { - // Match D0bar -> π- K+ - auto candD0MC = particlesMc.rawIteratorAt(particle.daughtersIds().front()); - // Printf("Checking D0bar -> π- K+"); - if (RecoDecay::isMatchedMCGen(particlesMc, candD0MC, static_cast(pdg::Code::kD0), std::array{+kPiPlus, -kKPlus}, true, &sign)) { - flag = sign * BIT(hf_cand_bplus::DecayType::BplusToD0Pi); - } + fillHfCand2Prong = true; + } // pion loop + if (fillHfCand2Prong) { // fill candD0 table only once per D0 candidate + hfCand2Prong(track0.globalIndex(), track1.globalIndex(), + indexHfReducedCollision, + trackParCovD0.getX(), trackParCovD0.getAlpha(), + trackParCovD0.getY(), trackParCovD0.getZ(), trackParCovD0.getSnp(), + trackParCovD0.getTgl(), trackParCovD0.getQ2Pt(), + candD0.xSecondaryVertex(), candD0.ySecondaryVertex(), candD0.zSecondaryVertex(), invMassD0, invMassD0bar); + hfCand2ProngCov(trackParCovD0.getSigmaY2(), trackParCovD0.getSigmaZY(), trackParCovD0.getSigmaZ2(), + trackParCovD0.getSigmaSnpY(), trackParCovD0.getSigmaSnpZ(), + trackParCovD0.getSigmaSnp2(), trackParCovD0.getSigmaTglY(), trackParCovD0.getSigmaTglZ(), + trackParCovD0.getSigmaTglSnp(), trackParCovD0.getSigmaTgl2(), + trackParCovD0.getSigma1PtY(), trackParCovD0.getSigma1PtZ(), trackParCovD0.getSigma1PtSnp(), + trackParCovD0.getSigma1PtTgl(), trackParCovD0.getSigma1Pt2()); + if constexpr (withMl) { + hfCand2ProngMl(candD0.mlProbD0()[0], candD0.mlProbD0()[1], candD0.mlProbD0()[2]); } + fillHfReducedCollision = true; + } + } // candsD loop + registry.fill(HIST("hEvents"), 1 + Event::Processed); + if (!fillHfReducedCollision) { + registry.fill(HIST("hEvents"), 1 + Event::NoD0PiSelected); + return; + } + registry.fill(HIST("hEvents"), 1 + Event::D0PiSelected); + // fill collision table if it contains a D0Pi pair a minima + hfReducedCollision(collision.posX(), collision.posY(), collision.posZ(), + collision.covXX(), collision.covXY(), collision.covYY(), + collision.covXZ(), collision.covYZ(), collision.covZZ(), + bz); + } - // save information for B+ task - if (!TESTBIT(std::abs(flag), hf_cand_bplus::DecayType::BplusToD0Pi)) { - continue; + void runMcGen(aod::McParticles const& particlesMc) + { + // Match generated particles. + for (const auto& particle : particlesMc) { + int8_t sign{0}; + int8_t flag{0}; + // B+ → D0bar π+ + if (RecoDecay::isMatchedMCGen(particlesMc, particle, pdg::Code::kBPlus, std::array{static_cast(pdg::Code::kD0), +kPiPlus}, true)) { + // Match D0bar -> π- K+ + auto candD0MC = particlesMc.rawIteratorAt(particle.daughtersIds().front()); + // Printf("Checking D0bar -> π- K+"); + if (RecoDecay::isMatchedMCGen(particlesMc, candD0MC, static_cast(pdg::Code::kD0), std::array{+kPiPlus, -kKPlus}, true, &sign)) { + flag = sign * BIT(hf_cand_bplus::DecayType::BplusToD0Pi); } + } - auto ptParticle = particle.pt(); - auto yParticle = RecoDecay::y(std::array{particle.px(), particle.py(), particle.pz()}, massBplus); - auto etaParticle = particle.eta(); - - std::array ptProngs; - std::array yProngs; - std::array etaProngs; - int counter = 0; - for (const auto& daught : particle.template daughters_as

()) { - ptProngs[counter] = daught.pt(); - etaProngs[counter] = daught.eta(); - yProngs[counter] = RecoDecay::y(std::array{daught.px(), daught.py(), daught.pz()}, pdg->Mass(daught.pdgCode())); - counter++; - } - rowHfBpMcGenReduced(flag, ptParticle, yParticle, etaParticle, - ptProngs[0], yProngs[0], etaProngs[0], - ptProngs[1], yProngs[1], etaProngs[1]); - } // gen - } + // save information for B+ task + if (!TESTBIT(std::abs(flag), hf_cand_bplus::DecayType::BplusToD0Pi)) { + continue; + } + + auto ptParticle = particle.pt(); + auto yParticle = RecoDecay::y(std::array{particle.px(), particle.py(), particle.pz()}, massBplus); + 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(std::array{daught.px(), daught.py(), daught.pz()}, pdg->Mass(daught.pdgCode())); + counter++; + } + rowHfBpMcGenReduced(flag, ptParticle, yParticle, etaParticle, + ptProngs[0], yProngs[0], etaProngs[0], + ptProngs[1], yProngs[1], etaProngs[1]); + } // gen } void processData(aod::Collisions const& collisions, @@ -450,9 +444,47 @@ struct HfDataCreatorD0PiReduced { TracksPIDWithSel const& tracks, aod::BCsWithTimestamps const& bcs) { - runDataCreation(collisions, candsD0, trackIndices, tracks, tracks, bcs); + // store configurables needed for B+ workflow + if (!isHfCandBplusConfigFilled) { + rowCandidateConfig(selectionFlagD0.value, selectionFlagD0bar.value, invMassWindowD0Pi.value); + isHfCandBplusConfigFilled = true; + } + + // handle normalization by the right number of collisions + hfCollisionCounter(collisions.tableSize()); + + for (const auto& collision : collisions) { + auto thisCollId = collision.globalIndex(); + auto candsDThisColl = candsD0.sliceBy(candsDPerCollision, thisCollId); + auto trackIdsThisCollision = trackIndices.sliceBy(trackIndicesPerCollision, thisCollId); + runDataCreation(collision, candsDThisColl, trackIdsThisCollision, tracks, tracks, bcs); + } } - PROCESS_SWITCH(HfDataCreatorD0PiReduced, processData, "Process without MC info", true); + PROCESS_SWITCH(HfDataCreatorD0PiReduced, processData, "Process without MC info and without ML info", true); + + void processDataWithMl(aod::Collisions const& collisions, + CandsDFilteredWithMl const& candsD0, + aod::TrackAssoc const& trackIndices, + TracksPIDWithSel const& tracks, + aod::BCsWithTimestamps const& bcs) + { + // store configurables needed for B+ workflow + if (!isHfCandBplusConfigFilled) { + rowCandidateConfig(selectionFlagD0.value, selectionFlagD0bar.value, invMassWindowD0Pi.value); + isHfCandBplusConfigFilled = true; + } + + // handle normalization by the right number of collisions + hfCollisionCounter(collisions.tableSize()); + + for (const auto& collision : collisions) { + auto thisCollId = collision.globalIndex(); + auto candsDThisColl = candsD0.sliceBy(candsDPerCollision, thisCollId); + auto trackIdsThisCollision = trackIndices.sliceBy(trackIndicesPerCollision, thisCollId); + runDataCreation(collision, candsDThisColl, trackIdsThisCollision, tracks, tracks, bcs); + } + } + PROCESS_SWITCH(HfDataCreatorD0PiReduced, processDataWithMl, "Process without MC info and with ML info", false); void processMc(aod::Collisions const& collisions, CandsDFiltered const& candsD0, @@ -461,9 +493,52 @@ struct HfDataCreatorD0PiReduced { aod::McParticles const& particlesMc, aod::BCsWithTimestamps const& bcs) { - runDataCreation(collisions, candsD0, trackIndices, tracks, particlesMc, bcs); + // store configurables needed for B+ workflow + if (!isHfCandBplusConfigFilled) { + rowCandidateConfig(selectionFlagD0.value, selectionFlagD0bar.value, invMassWindowD0Pi.value); + isHfCandBplusConfigFilled = true; + } + + // handle normalization by the right number of collisions + hfCollisionCounter(collisions.tableSize()); + + for (const auto& collision : collisions) { + auto thisCollId = collision.globalIndex(); + auto candsDThisColl = candsD0.sliceBy(candsDPerCollision, thisCollId); + auto trackIdsThisCollision = trackIndices.sliceBy(trackIndicesPerCollision, thisCollId); + runDataCreation(collision, candsDThisColl, trackIdsThisCollision, tracks, particlesMc, bcs); + } + + runMcGen(particlesMc); + } + PROCESS_SWITCH(HfDataCreatorD0PiReduced, processMc, "Process with MC info and without ML info", false); + + void processMcWithMl(aod::Collisions const& collisions, + CandsDFilteredWithMl const& candsD0, + aod::TrackAssoc const& trackIndices, + TracksPIDWithSelAndMc const& tracks, + aod::McParticles const& particlesMc, + aod::BCsWithTimestamps const& bcs) + { + // store configurables needed for B+ workflow + if (!isHfCandBplusConfigFilled) { + rowCandidateConfig(selectionFlagD0.value, selectionFlagD0bar.value, invMassWindowD0Pi.value); + isHfCandBplusConfigFilled = true; + } + + // handle normalization by the right number of collisions + hfCollisionCounter(collisions.tableSize()); + + for (const auto& collision : collisions) { + auto thisCollId = collision.globalIndex(); + auto candsDThisColl = candsD0.sliceBy(candsDPerCollision, thisCollId); + auto trackIdsThisCollision = trackIndices.sliceBy(trackIndicesPerCollision, thisCollId); + runDataCreation(collision, candsDThisColl, trackIdsThisCollision, tracks, particlesMc, bcs); + } + + runMcGen(particlesMc); } - PROCESS_SWITCH(HfDataCreatorD0PiReduced, processMc, "Process with MC info", false); + PROCESS_SWITCH(HfDataCreatorD0PiReduced, processMcWithMl, "Process with MC info and with ML info", false); }; // struct diff --git a/PWGHF/D2H/TableProducer/dataCreatorDplusPiReduced.cxx b/PWGHF/D2H/TableProducer/dataCreatorDplusPiReduced.cxx index 359e9a8c86d..db2fec04dfb 100644 --- a/PWGHF/D2H/TableProducer/dataCreatorDplusPiReduced.cxx +++ b/PWGHF/D2H/TableProducer/dataCreatorDplusPiReduced.cxx @@ -13,6 +13,7 @@ /// \brief Creation of Dplus-Pi pairs /// /// \author Alexandre Bigot , IPHC Strasbourg +/// \author Fabrizio Grosa , CERN #include @@ -56,6 +57,7 @@ struct HfDataCreatorDplusPiReduced { Produces hfTrackPidPion; Produces hfCand3Prong; Produces hfCand3ProngCov; + Produces hfCand3ProngMl; Produces rowCandidateConfig; Produces rowHfDPiMcRecReduced; Produces rowHfB0McGenReduced; @@ -112,10 +114,12 @@ struct HfDataCreatorDplusPiReduced { using TracksPIDWithSel = soa::Join; using TracksPIDWithSelAndMc = soa::Join; using CandsDFiltered = soa::Filtered>; + using CandsDFilteredWithMl = soa::Filtered>; Filter filterSelectCandidates = (aod::hf_sel_candidate_dplus::isSelDplusToPiKPi >= selectionFlagD); Preslice candsDPerCollision = aod::track_association::collisionId; + Preslice candsDPerCollisionWithMl = aod::track_association::collisionId; Preslice trackIndicesPerCollision = aod::track_association::collisionId; HistogramRegistry registry{"registry"}; @@ -211,241 +215,231 @@ struct HfDataCreatorDplusPiReduced { return true; } - template - void runDataCreation(aod::Collisions const& collisions, - CandsDFiltered const& candsD, + template + void runDataCreation(aod::Collision const& collision, + C const& candsD, aod::TrackAssoc const& trackIndices, T const& tracks, P const& particlesMc, aod::BCsWithTimestamps const& bcs) { - // store configurables needed for B0 workflow - if (!isHfCandB0ConfigFilled) { - rowCandidateConfig(selectionFlagD.value, invMassWindowDPi.value); - isHfCandB0ConfigFilled = true; + // helpers for ReducedTables filling + int indexHfReducedCollision = hfReducedCollision.lastIndex() + 1; + // std::map where the key is the track.globalIndex() and + // the value is the track index in the table of the selected pions + std::map selectedTracksPion; + bool fillHfReducedCollision = false; + + auto primaryVertex = getPrimaryVertex(collision); + + // Set the magnetic field from ccdb. + // The static instance of the propagator was already modified in the HFTrackIndexSkimCreator, + // but this is not true when running on Run2 data/MC already converted into AO2Ds. + auto bc = collision.bc_as(); + if (runNumber != bc.runNumber()) { + LOG(info) << ">>>>>>>>>>>> Current run number: " << runNumber; + initCCDB(bc, runNumber, ccdb, isRun2 ? ccdbPathGrp : ccdbPathGrpMag, lut, isRun2); + bz = o2::base::Propagator::Instance()->getNominalBz(); + LOG(info) << ">>>>>>>>>>>> Magnetic field: " << bz; } + df3.setBz(bz); + + auto thisCollId = collision.globalIndex(); + for (const auto& candD : candsD) { + int indexHfCand3Prong = hfCand3Prong.lastIndex() + 1; + bool fillHfCand3Prong = false; + float invMassD = hfHelper.invMassDplusToPiKPi(candD); + + registry.fill(HIST("hMassDToPiKPi"), invMassD); + registry.fill(HIST("hPtD"), candD.pt()); + registry.fill(HIST("hCPAD"), candD.cpa()); + + // track0 <-> pi, track1 <-> K, track2 <-> pi + auto track0 = candD.template prong0_as(); + auto track1 = candD.template prong1_as(); + auto track2 = candD.template prong2_as(); + auto trackParCov0 = getTrackParCov(track0); + auto trackParCov1 = getTrackParCov(track1); + auto trackParCov2 = getTrackParCov(track2); + + std::array pVec0 = {track0.px(), track0.py(), track0.pz()}; + std::array pVec1 = {track1.px(), track1.py(), track1.pz()}; + std::array pVec2 = {track2.px(), track2.py(), track2.pz()}; + + auto dca0 = o2::dataformats::DCA(track0.dcaXY(), track0.dcaZ(), track0.cYY(), track0.cZY(), track0.cZZ()); + auto dca1 = o2::dataformats::DCA(track1.dcaXY(), track1.dcaZ(), track1.cYY(), track1.cZY(), track1.cZZ()); + auto dca2 = o2::dataformats::DCA(track2.dcaXY(), track2.dcaZ(), track2.cYY(), track2.cZY(), track2.cZZ()); + + // repropagate tracks to this collision if needed + if (track0.collisionId() != thisCollId) { + trackParCov0.propagateToDCA(primaryVertex, bz, &dca0); + } - // handle normalization by the right number of collisions - hfCollisionCounter(collisions.tableSize()); + if (track1.collisionId() != thisCollId) { + trackParCov1.propagateToDCA(primaryVertex, bz, &dca1); + } - for (const auto& collision : collisions) { + if (track2.collisionId() != thisCollId) { + trackParCov2.propagateToDCA(primaryVertex, bz, &dca2); + } - // helpers for ReducedTables filling - int indexHfReducedCollision = hfReducedCollision.lastIndex() + 1; - // std::map where the key is the track.globalIndex() and - // the value is the track index in the table of the selected pions - std::map selectedTracksPion; - bool fillHfReducedCollision = false; - - auto primaryVertex = getPrimaryVertex(collision); - - // Set the magnetic field from ccdb. - // The static instance of the propagator was already modified in the HFTrackIndexSkimCreator, - // but this is not true when running on Run2 data/MC already converted into AO2Ds. - auto bc = collision.bc_as(); - if (runNumber != bc.runNumber()) { - LOG(info) << ">>>>>>>>>>>> Current run number: " << runNumber; - initCCDB(bc, runNumber, ccdb, isRun2 ? ccdbPathGrp : ccdbPathGrpMag, lut, isRun2); - bz = o2::base::Propagator::Instance()->getNominalBz(); - LOG(info) << ">>>>>>>>>>>> Magnetic field: " << bz; + // --------------------------------- + // reconstruct 3-prong secondary vertex (D±) + if (df3.process(trackParCov0, trackParCov1, trackParCov2) == 0) { + continue; } - df3.setBz(bz); - auto thisCollId = collision.globalIndex(); - auto candsDThisColl = candsD.sliceBy(candsDPerCollision, thisCollId); - for (const auto& candD : candsDThisColl) { - int indexHfCand3Prong = hfCand3Prong.lastIndex() + 1; - bool fillHfCand3Prong = false; - float invMassD = hfHelper.invMassDplusToPiKPi(candD); - - registry.fill(HIST("hMassDToPiKPi"), invMassD); - registry.fill(HIST("hPtD"), candD.pt()); - registry.fill(HIST("hCPAD"), candD.cpa()); - - // track0 <-> pi, track1 <-> K, track2 <-> pi - auto track0 = candD.template prong0_as(); - auto track1 = candD.template prong1_as(); - auto track2 = candD.template prong2_as(); - auto trackParCov0 = getTrackParCov(track0); - auto trackParCov1 = getTrackParCov(track1); - auto trackParCov2 = getTrackParCov(track2); - - std::array pVec0 = {track0.px(), track0.py(), track0.pz()}; - std::array pVec1 = {track1.px(), track1.py(), track1.pz()}; - std::array pVec2 = {track2.px(), track2.py(), track2.pz()}; - - auto dca0 = o2::dataformats::DCA(track0.dcaXY(), track0.dcaZ(), track0.cYY(), track0.cZY(), track0.cZZ()); - auto dca1 = o2::dataformats::DCA(track1.dcaXY(), track1.dcaZ(), track1.cYY(), track1.cZY(), track1.cZZ()); - auto dca2 = o2::dataformats::DCA(track2.dcaXY(), track2.dcaZ(), track2.cYY(), track2.cZY(), track2.cZZ()); - - // repropagate tracks to this collision if needed - if (track0.collisionId() != thisCollId) { - trackParCov0.propagateToDCA(primaryVertex, bz, &dca0); - } + const auto& secondaryVertexD = df3.getPCACandidate(); + // propagate the 3 prongs to the secondary vertex + trackParCov0.propagateTo(secondaryVertexD[0], bz); + trackParCov1.propagateTo(secondaryVertexD[0], bz); + trackParCov2.propagateTo(secondaryVertexD[0], bz); - if (track1.collisionId() != thisCollId) { - trackParCov1.propagateToDCA(primaryVertex, bz, &dca1); - } + // update pVec of tracks + df3.getTrack(0).getPxPyPzGlo(pVec0); + df3.getTrack(1).getPxPyPzGlo(pVec1); + df3.getTrack(2).getPxPyPzGlo(pVec2); - if (track2.collisionId() != thisCollId) { - trackParCov2.propagateToDCA(primaryVertex, bz, &dca2); - } + // D∓ → π∓ K± π∓ + std::array pVecPiK = RecoDecay::pVec(pVec0, pVec1); + std::array pVecD = RecoDecay::pVec(pVec0, pVec1, pVec2); + auto trackParCovPiK = o2::dataformats::V0(df3.getPCACandidatePos(), pVecPiK, df3.calcPCACovMatrixFlat(), trackParCov0, trackParCov1); + auto trackParCovD = o2::dataformats::V0(df3.getPCACandidatePos(), pVecD, df3.calcPCACovMatrixFlat(), trackParCovPiK, trackParCov2); - // --------------------------------- - // reconstruct 3-prong secondary vertex (D±) - if (df3.process(trackParCov0, trackParCov1, trackParCov2) == 0) { + for (const auto& trackId : trackIndices) { + auto trackPion = trackId.template track_as(); + + // apply selections on pion tracks + if (!isPionSelected(trackPion, track0, track1, track2) || !isSelectedTrackDCA(trackPion)) { + continue; + } + registry.fill(HIST("hPtPion"), trackPion.pt()); + std::array pVecPion = {trackPion.px(), trackPion.py(), trackPion.pz()}; + // compute invariant mass square and apply selection + auto invMass2DPi = RecoDecay::m2(std::array{pVecD, pVecPion}, std::array{massD, massPi}); + if ((invMass2DPi < invMass2DPiMin) || (invMass2DPi > invMass2DPiMax)) { continue; } - const auto& secondaryVertexD = df3.getPCACandidate(); - // propagate the 3 prongs to the secondary vertex - trackParCov0.propagateTo(secondaryVertexD[0], bz); - trackParCov1.propagateTo(secondaryVertexD[0], bz); - trackParCov2.propagateTo(secondaryVertexD[0], bz); - - // update pVec of tracks - df3.getTrack(0).getPxPyPzGlo(pVec0); - df3.getTrack(1).getPxPyPzGlo(pVec1); - df3.getTrack(2).getPxPyPzGlo(pVec2); - - // D∓ → π∓ K± π∓ - std::array pVecPiK = RecoDecay::pVec(pVec0, pVec1); - std::array pVecD = RecoDecay::pVec(pVec0, pVec1, pVec2); - auto trackParCovPiK = o2::dataformats::V0(df3.getPCACandidatePos(), pVecPiK, df3.calcPCACovMatrixFlat(), trackParCov0, trackParCov1); - auto trackParCovD = o2::dataformats::V0(df3.getPCACandidatePos(), pVecD, df3.calcPCACovMatrixFlat(), trackParCovPiK, trackParCov2); - - auto trackIdsThisCollision = trackIndices.sliceBy(trackIndicesPerCollision, thisCollId); - for (const auto& trackId : trackIdsThisCollision) { - auto trackPion = trackId.template track_as(); - - // apply selections on pion tracks - if (!isPionSelected(trackPion, track0, track1, track2) || !isSelectedTrackDCA(trackPion)) { - continue; - } - registry.fill(HIST("hPtPion"), trackPion.pt()); - std::array pVecPion = {trackPion.px(), trackPion.py(), trackPion.pz()}; - // compute invariant mass square and apply selection - auto invMass2DPi = RecoDecay::m2(std::array{pVecD, pVecPion}, std::array{massD, massPi}); - if ((invMass2DPi < invMass2DPiMin) || (invMass2DPi > invMass2DPiMax)) { - continue; - } - - // fill Pion tracks table - // if information on track already stored, go to next track - if (!selectedTracksPion.count(trackPion.globalIndex())) { - hfTrackPion(trackPion.globalIndex(), indexHfReducedCollision, - trackPion.x(), trackPion.alpha(), - trackPion.y(), trackPion.z(), trackPion.snp(), - trackPion.tgl(), trackPion.signed1Pt()); - hfTrackCovPion(trackPion.cYY(), trackPion.cZY(), trackPion.cZZ(), - trackPion.cSnpY(), trackPion.cSnpZ(), - trackPion.cSnpSnp(), trackPion.cTglY(), trackPion.cTglZ(), - trackPion.cTglSnp(), trackPion.cTglTgl(), - trackPion.c1PtY(), trackPion.c1PtZ(), trackPion.c1PtSnp(), - trackPion.c1PtTgl(), trackPion.c1Pt21Pt2()); - hfTrackPidPion(trackPion.hasTPC(), trackPion.hasTOF(), - trackPion.tpcNSigmaPi(), trackPion.tofNSigmaPi()); - // add trackPion.globalIndex() to a list - // to keep memory of the pions filled in the table and avoid refilling them if they are paired to another D candidate - // and keep track of their index in hfTrackPion for McRec purposes - selectedTracksPion[trackPion.globalIndex()] = hfTrackPion.lastIndex(); - } + // fill Pion tracks table + // if information on track already stored, go to next track + if (!selectedTracksPion.count(trackPion.globalIndex())) { + hfTrackPion(trackPion.globalIndex(), indexHfReducedCollision, + trackPion.x(), trackPion.alpha(), + trackPion.y(), trackPion.z(), trackPion.snp(), + trackPion.tgl(), trackPion.signed1Pt()); + hfTrackCovPion(trackPion.cYY(), trackPion.cZY(), trackPion.cZZ(), + trackPion.cSnpY(), trackPion.cSnpZ(), + trackPion.cSnpSnp(), trackPion.cTglY(), trackPion.cTglZ(), + trackPion.cTglSnp(), trackPion.cTglTgl(), + trackPion.c1PtY(), trackPion.c1PtZ(), trackPion.c1PtSnp(), + trackPion.c1PtTgl(), trackPion.c1Pt21Pt2()); + hfTrackPidPion(trackPion.hasTPC(), trackPion.hasTOF(), + trackPion.tpcNSigmaPi(), trackPion.tofNSigmaPi()); + // add trackPion.globalIndex() to a list + // to keep memory of the pions filled in the table and avoid refilling them if they are paired to another D candidate + // and keep track of their index in hfTrackPion for McRec purposes + selectedTracksPion[trackPion.globalIndex()] = hfTrackPion.lastIndex(); + } - if constexpr (doMc) { - // we check the MC matching to be stored - auto arrayDaughtersD = std::array{track0, track1, track2}; - auto arrayDaughtersB0 = std::array{track0, track1, track2, trackPion}; - int8_t sign{0}; - int8_t flag{0}; - int8_t debug{0}; - // B0 → D- π+ → (π- K+ π-) π+ - auto indexRec = RecoDecay::getMatchedMCRec(particlesMc, arrayDaughtersB0, pdg::Code::kB0, std::array{-kPiPlus, +kKPlus, -kPiPlus, +kPiPlus}, true, &sign, 3); + if constexpr (doMc) { + // we check the MC matching to be stored + auto arrayDaughtersD = std::array{track0, track1, track2}; + auto arrayDaughtersB0 = std::array{track0, track1, track2, trackPion}; + int8_t sign{0}; + int8_t flag{0}; + int8_t debug{0}; + // B0 → D- π+ → (π- K+ π-) π+ + auto indexRec = RecoDecay::getMatchedMCRec(particlesMc, arrayDaughtersB0, pdg::Code::kB0, std::array{-kPiPlus, +kKPlus, -kPiPlus, +kPiPlus}, true, &sign, 3); + if (indexRec > -1) { + // D- → π- K+ π- + // Printf("Checking D- → π- K+ π-"); + indexRec = RecoDecay::getMatchedMCRec(particlesMc, arrayDaughtersD, pdg::Code::kDMinus, std::array{-kPiPlus, +kKPlus, -kPiPlus}, true, &sign, 2); if (indexRec > -1) { - // D- → π- K+ π- - // Printf("Checking D- → π- K+ π-"); - indexRec = RecoDecay::getMatchedMCRec(particlesMc, arrayDaughtersD, pdg::Code::kDMinus, std::array{-kPiPlus, +kKPlus, -kPiPlus}, true, &sign, 2); - if (indexRec > -1) { - flag = sign * BIT(hf_cand_b0::DecayType::B0ToDPi); - } else { - debug = 1; - LOGF(debug, "B0 decays in the expected final state but the condition on the intermediate state is not fulfilled"); - } + flag = sign * BIT(hf_cand_b0::DecayType::B0ToDPi); + } else { + debug = 1; + LOGF(debug, "B0 decays in the expected final state but the condition on the intermediate state is not fulfilled"); } - auto indexMother = RecoDecay::getMother(particlesMc, trackPion.template mcParticle_as

(), pdg::Code::kB0, true); - auto particleMother = particlesMc.rawIteratorAt(indexMother); - rowHfDPiMcRecReduced(indexHfCand3Prong, selectedTracksPion[trackPion.globalIndex()], flag, debug, particleMother.pt()); } - fillHfCand3Prong = true; - } // pion loop - if (fillHfCand3Prong) { // fill candDplus table only once per D candidate - hfCand3Prong(track0.globalIndex(), track1.globalIndex(), track2.globalIndex(), - indexHfReducedCollision, - trackParCovD.getX(), trackParCovD.getAlpha(), - trackParCovD.getY(), trackParCovD.getZ(), trackParCovD.getSnp(), - trackParCovD.getTgl(), trackParCovD.getQ2Pt(), - candD.xSecondaryVertex(), candD.ySecondaryVertex(), candD.zSecondaryVertex(), invMassD); - hfCand3ProngCov(trackParCovD.getSigmaY2(), trackParCovD.getSigmaZY(), trackParCovD.getSigmaZ2(), - trackParCovD.getSigmaSnpY(), trackParCovD.getSigmaSnpZ(), - trackParCovD.getSigmaSnp2(), trackParCovD.getSigmaTglY(), trackParCovD.getSigmaTglZ(), - trackParCovD.getSigmaTglSnp(), trackParCovD.getSigmaTgl2(), - trackParCovD.getSigma1PtY(), trackParCovD.getSigma1PtZ(), trackParCovD.getSigma1PtSnp(), - trackParCovD.getSigma1PtTgl(), trackParCovD.getSigma1Pt2()); - fillHfReducedCollision = true; + auto indexMother = RecoDecay::getMother(particlesMc, trackPion.template mcParticle_as

(), pdg::Code::kB0, true); + auto particleMother = particlesMc.rawIteratorAt(indexMother); + rowHfDPiMcRecReduced(indexHfCand3Prong, selectedTracksPion[trackPion.globalIndex()], flag, debug, particleMother.pt()); } - } // candsD loop - - registry.fill(HIST("hEvents"), 1 + Event::Processed); - if (!fillHfReducedCollision) { - registry.fill(HIST("hEvents"), 1 + Event::NoDPiSelected); - continue; - } - registry.fill(HIST("hEvents"), 1 + Event::DPiSelected); - // fill collision table if it contains a DPi pair a minima - hfReducedCollision(collision.posX(), collision.posY(), collision.posZ(), - collision.covXX(), collision.covXY(), collision.covYY(), - collision.covXZ(), collision.covYZ(), collision.covZZ(), - bz); - } // collision - - if constexpr (doMc) { - // Match generated particles. - for (const auto& particle : particlesMc) { - int8_t sign{0}; - int8_t flag{0}; - // B0 → D- π+ - if (RecoDecay::isMatchedMCGen(particlesMc, particle, pdg::Code::kB0, std::array{-static_cast(pdg::Code::kDPlus), +kPiPlus}, true)) { - // Match D- -> π- K+ π- - auto candDMC = particlesMc.rawIteratorAt(particle.daughtersIds().front()); - // Printf("Checking D- -> π- K+ π-"); - if (RecoDecay::isMatchedMCGen(particlesMc, candDMC, -static_cast(pdg::Code::kDPlus), std::array{-kPiPlus, +kKPlus, -kPiPlus}, true, &sign)) { - flag = sign * BIT(hf_cand_b0::DecayType::B0ToDPi); - } + fillHfCand3Prong = true; + } // pion loop + if (fillHfCand3Prong) { // fill candDplus table only once per D candidate + hfCand3Prong(track0.globalIndex(), track1.globalIndex(), track2.globalIndex(), + indexHfReducedCollision, + trackParCovD.getX(), trackParCovD.getAlpha(), + trackParCovD.getY(), trackParCovD.getZ(), trackParCovD.getSnp(), + trackParCovD.getTgl(), trackParCovD.getQ2Pt(), + candD.xSecondaryVertex(), candD.ySecondaryVertex(), candD.zSecondaryVertex(), invMassD); + hfCand3ProngCov(trackParCovD.getSigmaY2(), trackParCovD.getSigmaZY(), trackParCovD.getSigmaZ2(), + trackParCovD.getSigmaSnpY(), trackParCovD.getSigmaSnpZ(), + trackParCovD.getSigmaSnp2(), trackParCovD.getSigmaTglY(), trackParCovD.getSigmaTglZ(), + trackParCovD.getSigmaTglSnp(), trackParCovD.getSigmaTgl2(), + trackParCovD.getSigma1PtY(), trackParCovD.getSigma1PtZ(), trackParCovD.getSigma1PtSnp(), + trackParCovD.getSigma1PtTgl(), trackParCovD.getSigma1Pt2()); + if constexpr (withMl) { + hfCand3ProngMl(candD.mlProbDplusToPiKPi()[0], candD.mlProbDplusToPiKPi()[1], candD.mlProbDplusToPiKPi()[2]); } + fillHfReducedCollision = true; + } + } // candsD loop - // save information for B0 task - if (!TESTBIT(std::abs(flag), hf_cand_b0::DecayType::B0ToDPi)) { - continue; - } + registry.fill(HIST("hEvents"), 1 + Event::Processed); + if (!fillHfReducedCollision) { + registry.fill(HIST("hEvents"), 1 + Event::NoDPiSelected); + return; + } + registry.fill(HIST("hEvents"), 1 + Event::DPiSelected); + // fill collision table if it contains a DPi pair a minima + hfReducedCollision(collision.posX(), collision.posY(), collision.posZ(), + collision.covXX(), collision.covXY(), collision.covYY(), + collision.covXZ(), collision.covYZ(), collision.covZZ(), + bz); + } - auto ptParticle = particle.pt(); - auto yParticle = RecoDecay::y(std::array{particle.px(), particle.py(), particle.pz()}, massB0); - auto etaParticle = particle.eta(); - - std::array ptProngs; - std::array yProngs; - std::array etaProngs; - int counter = 0; - for (const auto& daught : particle.template daughters_as

()) { - ptProngs[counter] = daught.pt(); - etaProngs[counter] = daught.eta(); - yProngs[counter] = RecoDecay::y(std::array{daught.px(), daught.py(), daught.pz()}, pdg->Mass(daught.pdgCode())); - counter++; + void runMcGen(aod::McParticles const& particlesMc) + { + // Match generated particles. + for (const auto& particle : particlesMc) { + int8_t sign{0}; + int8_t flag{0}; + // B0 → D- π+ + if (RecoDecay::isMatchedMCGen(particlesMc, particle, pdg::Code::kB0, std::array{-static_cast(pdg::Code::kDPlus), +kPiPlus}, true)) { + // Match D- -> π- K+ π- + auto candDMC = particlesMc.rawIteratorAt(particle.daughtersIds().front()); + // Printf("Checking D- -> π- K+ π-"); + if (RecoDecay::isMatchedMCGen(particlesMc, candDMC, -static_cast(pdg::Code::kDPlus), std::array{-kPiPlus, +kKPlus, -kPiPlus}, true, &sign)) { + flag = sign * BIT(hf_cand_b0::DecayType::B0ToDPi); } - rowHfB0McGenReduced(flag, ptParticle, yParticle, etaParticle, - ptProngs[0], yProngs[0], etaProngs[0], - ptProngs[1], yProngs[1], etaProngs[1]); - } // gen - } + } + + // save information for B0 task + if (!TESTBIT(std::abs(flag), hf_cand_b0::DecayType::B0ToDPi)) { + continue; + } + + auto ptParticle = particle.pt(); + auto yParticle = RecoDecay::y(std::array{particle.px(), particle.py(), particle.pz()}, 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(std::array{daught.px(), daught.py(), daught.pz()}, pdg->Mass(daught.pdgCode())); + counter++; + } + rowHfB0McGenReduced(flag, ptParticle, yParticle, etaParticle, + ptProngs[0], yProngs[0], etaProngs[0], + ptProngs[1], yProngs[1], etaProngs[1]); + } // gen } void processData(aod::Collisions const& collisions, @@ -454,9 +448,47 @@ struct HfDataCreatorDplusPiReduced { TracksPIDWithSel const& tracks, aod::BCsWithTimestamps const& bcs) { - runDataCreation(collisions, candsD, trackIndices, tracks, tracks, bcs); + // store configurables needed for B0 workflow + if (!isHfCandB0ConfigFilled) { + rowCandidateConfig(selectionFlagD.value, invMassWindowDPi.value); + isHfCandB0ConfigFilled = true; + } + + // handle normalization by the right number of collisions + hfCollisionCounter(collisions.tableSize()); + + for (const auto& collision : collisions) { + auto thisCollId = collision.globalIndex(); + auto candsDThisColl = candsD.sliceBy(candsDPerCollision, thisCollId); + auto trackIdsThisCollision = trackIndices.sliceBy(trackIndicesPerCollision, thisCollId); + runDataCreation(collision, candsDThisColl, trackIdsThisCollision, tracks, tracks, bcs); + } } - PROCESS_SWITCH(HfDataCreatorDplusPiReduced, processData, "Process without MC info", true); + PROCESS_SWITCH(HfDataCreatorDplusPiReduced, processData, "Process without MC info and without ML info", true); + + void processDataWithMl(aod::Collisions const& collisions, + CandsDFilteredWithMl const& candsD, + aod::TrackAssoc const& trackIndices, + TracksPIDWithSel const& tracks, + aod::BCsWithTimestamps const& bcs) + { + // store configurables needed for B0 workflow + if (!isHfCandB0ConfigFilled) { + rowCandidateConfig(selectionFlagD.value, invMassWindowDPi.value); + isHfCandB0ConfigFilled = true; + } + + // handle normalization by the right number of collisions + hfCollisionCounter(collisions.tableSize()); + + for (const auto& collision : collisions) { + auto thisCollId = collision.globalIndex(); + auto candsDThisColl = candsD.sliceBy(candsDPerCollisionWithMl, thisCollId); + auto trackIdsThisCollision = trackIndices.sliceBy(trackIndicesPerCollision, thisCollId); + runDataCreation(collision, candsDThisColl, trackIdsThisCollision, tracks, tracks, bcs); + } + } + PROCESS_SWITCH(HfDataCreatorDplusPiReduced, processDataWithMl, "Process without MC info and with ML info", false); void processMc(aod::Collisions const& collisions, CandsDFiltered const& candsD, @@ -465,9 +497,50 @@ struct HfDataCreatorDplusPiReduced { aod::McParticles const& particlesMc, aod::BCsWithTimestamps const& bcs) { - runDataCreation(collisions, candsD, trackIndices, tracks, particlesMc, bcs); + // store configurables needed for B0 workflow + if (!isHfCandB0ConfigFilled) { + rowCandidateConfig(selectionFlagD.value, invMassWindowDPi.value); + isHfCandB0ConfigFilled = true; + } + + // handle normalization by the right number of collisions + hfCollisionCounter(collisions.tableSize()); + + for (const auto& collision : collisions) { + auto thisCollId = collision.globalIndex(); + auto candsDThisColl = candsD.sliceBy(candsDPerCollision, thisCollId); + auto trackIdsThisCollision = trackIndices.sliceBy(trackIndicesPerCollision, thisCollId); + runDataCreation(collision, candsDThisColl, trackIdsThisCollision, tracks, particlesMc, bcs); + } + runMcGen(particlesMc); + } + PROCESS_SWITCH(HfDataCreatorDplusPiReduced, processMc, "Process with MC info and without ML info", false); + + void processMcWithMl(aod::Collisions const& collisions, + CandsDFilteredWithMl const& candsD, + aod::TrackAssoc const& trackIndices, + TracksPIDWithSelAndMc const& tracks, + aod::McParticles const& particlesMc, + aod::BCsWithTimestamps const& bcs) + { + // store configurables needed for B0 workflow + if (!isHfCandB0ConfigFilled) { + rowCandidateConfig(selectionFlagD.value, invMassWindowDPi.value); + isHfCandB0ConfigFilled = true; + } + + // handle normalization by the right number of collisions + hfCollisionCounter(collisions.tableSize()); + + for (const auto& collision : collisions) { + auto thisCollId = collision.globalIndex(); + auto candsDThisColl = candsD.sliceBy(candsDPerCollisionWithMl, thisCollId); + auto trackIdsThisCollision = trackIndices.sliceBy(trackIndicesPerCollision, thisCollId); + runDataCreation(collision, candsD, trackIndices, tracks, particlesMc, bcs); + } + runMcGen(particlesMc); } - PROCESS_SWITCH(HfDataCreatorDplusPiReduced, processMc, "Process with MC info", false); + PROCESS_SWITCH(HfDataCreatorDplusPiReduced, processMcWithMl, "Process with MC info and with ML info", false); }; // struct WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) From 897c7f584875bfcb9ea89259e2c5742a5f093572 Mon Sep 17 00:00:00 2001 From: Fabrizio Grosa Date: Fri, 20 Oct 2023 11:22:43 +0200 Subject: [PATCH 2/2] Remove defaults to template arguments --- PWGHF/D2H/TableProducer/dataCreatorD0PiReduced.cxx | 2 +- PWGHF/D2H/TableProducer/dataCreatorDplusPiReduced.cxx | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/PWGHF/D2H/TableProducer/dataCreatorD0PiReduced.cxx b/PWGHF/D2H/TableProducer/dataCreatorD0PiReduced.cxx index cbf86f30fb6..64d51116fed 100644 --- a/PWGHF/D2H/TableProducer/dataCreatorD0PiReduced.cxx +++ b/PWGHF/D2H/TableProducer/dataCreatorD0PiReduced.cxx @@ -218,7 +218,7 @@ struct HfDataCreatorD0PiReduced { return true; } - template + template void runDataCreation(aod::Collision const& collision, C const& candsD0, aod::TrackAssoc const& trackIndices, diff --git a/PWGHF/D2H/TableProducer/dataCreatorDplusPiReduced.cxx b/PWGHF/D2H/TableProducer/dataCreatorDplusPiReduced.cxx index db2fec04dfb..e92be945f36 100644 --- a/PWGHF/D2H/TableProducer/dataCreatorDplusPiReduced.cxx +++ b/PWGHF/D2H/TableProducer/dataCreatorDplusPiReduced.cxx @@ -215,7 +215,7 @@ struct HfDataCreatorDplusPiReduced { return true; } - template + template void runDataCreation(aod::Collision const& collision, C const& candsD, aod::TrackAssoc const& trackIndices,