From f347b6bf54e06ce7271e963ebf0d065f40781962 Mon Sep 17 00:00:00 2001 From: Antonio Palasciano <52152842+apalasciano@users.noreply.github.com> Date: Mon, 13 May 2024 21:39:13 +0200 Subject: [PATCH] PWGHF: B+ reduced workflow with ML (#5964) * updateDraftBplusReduced * clang * hfflag removal * Extending hfflag removal to B0 --- PWGHF/Core/HfMlResponseBplusToD0Pi.h | 221 +++++ PWGHF/D2H/DataModel/ReducedDataModel.h | 29 + PWGHF/D2H/TableProducer/CMakeLists.txt | 2 +- .../candidateCreatorB0Reduced.cxx | 3 +- .../candidateCreatorBplusReduced.cxx | 277 ++++-- .../candidateSelectorB0ToDPiReduced.cxx | 24 +- .../candidateSelectorBplusToD0PiReduced.cxx | 152 +++- .../dataCreatorCharmHadPiReduced.cxx | 37 +- PWGHF/D2H/Tasks/taskB0.cxx | 6 - PWGHF/D2H/Tasks/taskB0Reduced.cxx | 24 - PWGHF/D2H/Tasks/taskBplus.cxx | 8 +- PWGHF/D2H/Tasks/taskBplusReduced.cxx | 830 +++++++++++++----- .../DataModel/CandidateReconstructionTables.h | 7 +- PWGHF/DataModel/CandidateSelectionTables.h | 6 +- PWGHF/TableProducer/candidateCreatorB0.cxx | 5 +- PWGHF/TableProducer/candidateCreatorBplus.cxx | 5 +- .../candidateSelectorB0ToDPi.cxx | 9 - .../candidateSelectorBplusToD0Pi.cxx | 6 - .../TableProducer/treeCreatorBplusToD0Pi.cxx | 465 ++++++---- 19 files changed, 1535 insertions(+), 581 deletions(-) create mode 100644 PWGHF/Core/HfMlResponseBplusToD0Pi.h diff --git a/PWGHF/Core/HfMlResponseBplusToD0Pi.h b/PWGHF/Core/HfMlResponseBplusToD0Pi.h new file mode 100644 index 00000000000..8d47d4ade89 --- /dev/null +++ b/PWGHF/Core/HfMlResponseBplusToD0Pi.h @@ -0,0 +1,221 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file HfMlResponseBplusToD0Pi.h +/// \brief Class to compute the ML response for B± → D0(bar) π± analysis selections +/// \author Antonio Palasciano , INFN Bari + +#ifndef PWGHF_CORE_HFMLRESPONSEBPLUSTOD0PI_H_ +#define PWGHF_CORE_HFMLRESPONSEBPLUSTOD0PI_H_ + +#include +#include +#include + +#include "PWGHF/Core/HfMlResponse.h" + +// Fill the map of available input features +// the key is the feature's name (std::string) +// the value is the corresponding value in EnumInputFeatures +#define FILL_MAP_BPLUS(FEATURE) \ + { \ +#FEATURE, static_cast < uint8_t>(InputFeaturesBplusToD0Pi::FEATURE) \ + } + +// Check if the index of mCachedIndices (index associated to a FEATURE) +// matches the entry in EnumInputFeatures associated to this FEATURE +// if so, the inputFeatures vector is filled with the FEATURE's value +// by calling the corresponding GETTER from OBJECT +#define CHECK_AND_FILL_VEC_BPLUS_FULL(OBJECT, FEATURE, GETTER) \ + case static_cast(InputFeaturesBplusToD0Pi::FEATURE): { \ + inputFeatures.emplace_back(OBJECT.GETTER()); \ + break; \ + } + +// Check if the index of mCachedIndices (index associated to a FEATURE) +// matches the entry in EnumInputFeatures associated to this FEATURE +// if so, the inputFeatures vector is filled with the FEATURE's value +// by calling the GETTER function taking OBJECT in argument +#define CHECK_AND_FILL_VEC_BPLUS_FUNC(OBJECT, FEATURE, GETTER) \ + case static_cast(InputFeaturesBplusToD0Pi::FEATURE): { \ + inputFeatures.emplace_back(GETTER(OBJECT)); \ + break; \ + } + +// Specific case of CHECK_AND_FILL_VEC_BPLUS_FULL(OBJECT, FEATURE, GETTER) +// where OBJECT is named candidate and FEATURE = GETTER +#define CHECK_AND_FILL_VEC_BPLUS(GETTER) \ + case static_cast(InputFeaturesBplusToD0Pi::GETTER): { \ + inputFeatures.emplace_back(candidate.GETTER()); \ + break; \ + } + +namespace o2::pid_tpc_tof_utils +{ +template +float getTpcTofNSigmaPi1(const T1& prong1) +{ + float defaultNSigma = -999.f; // -999.f is the default value set in TPCPIDResponse.h and PIDTOF.h + + bool hasTpc = prong1.hasTPC(); + bool hasTof = prong1.hasTOF(); + + if (hasTpc && hasTof) { + float tpcNSigma = prong1.tpcNSigmaPi(); + float tofNSigma = prong1.tofNSigmaPi(); + return sqrt(.5f * tpcNSigma * tpcNSigma + .5f * tofNSigma * tofNSigma); + } + if (hasTpc) { + return abs(prong1.tpcNSigmaPi()); + } + if (hasTof) { + return abs(prong1.tofNSigmaPi()); + } + return defaultNSigma; +} +} // namespace o2::pid_tpc_tof_utils + +namespace o2::analysis +{ + +enum class InputFeaturesBplusToD0Pi : uint8_t { + ptProng0 = 0, + ptProng1, + impactParameter0, + impactParameter1, + impactParameterProduct, + chi2PCA, + decayLength, + decayLengthXY, + decayLengthNormalised, + decayLengthXYNormalised, + cpa, + cpaXY, + maxNormalisedDeltaIP, + prong0MlScoreBkg, + prong0MlScorePrompt, + prong0MlScoreNonprompt, + tpcNSigmaPi1, + tofNSigmaPi1, + tpcTofNSigmaPi1 +}; + +template +class HfMlResponseBplusToD0Pi : public HfMlResponse +{ + public: + /// Default constructor + HfMlResponseBplusToD0Pi() = default; + /// Default destructor + virtual ~HfMlResponseBplusToD0Pi() = default; + + /// Method to get the input features vector needed for ML inference + /// \param candidate is the B+ candidate + /// \param prong1 is the candidate's prong1 + /// \return inputFeatures vector + template + std::vector getInputFeatures(T1 const& candidate, + T2 const& prong1) + { + std::vector inputFeatures; + + for (const auto& idx : MlResponse::mCachedIndices) { + if constexpr (withDmesMl) { + switch (idx) { + CHECK_AND_FILL_VEC_BPLUS(ptProng0); + CHECK_AND_FILL_VEC_BPLUS(ptProng1); + CHECK_AND_FILL_VEC_BPLUS(impactParameter0); + CHECK_AND_FILL_VEC_BPLUS(impactParameter1); + CHECK_AND_FILL_VEC_BPLUS(impactParameterProduct); + CHECK_AND_FILL_VEC_BPLUS(chi2PCA); + CHECK_AND_FILL_VEC_BPLUS(decayLength); + CHECK_AND_FILL_VEC_BPLUS(decayLengthXY); + CHECK_AND_FILL_VEC_BPLUS(decayLengthNormalised); + CHECK_AND_FILL_VEC_BPLUS(decayLengthXYNormalised); + CHECK_AND_FILL_VEC_BPLUS(cpa); + CHECK_AND_FILL_VEC_BPLUS(cpaXY); + CHECK_AND_FILL_VEC_BPLUS(maxNormalisedDeltaIP); + CHECK_AND_FILL_VEC_BPLUS(prong0MlScoreBkg); + CHECK_AND_FILL_VEC_BPLUS(prong0MlScorePrompt); + CHECK_AND_FILL_VEC_BPLUS(prong0MlScoreNonprompt); + // TPC PID variable + CHECK_AND_FILL_VEC_BPLUS_FULL(prong1, tpcNSigmaPi1, tpcNSigmaPi); + // TOF PID variable + CHECK_AND_FILL_VEC_BPLUS_FULL(prong1, tofNSigmaPi1, tofNSigmaPi); + // Combined PID variables + CHECK_AND_FILL_VEC_BPLUS_FUNC(prong1, tpcTofNSigmaPi1, o2::pid_tpc_tof_utils::getTpcTofNSigmaPi1); + } + } else { + switch (idx) { + CHECK_AND_FILL_VEC_BPLUS(ptProng0); + CHECK_AND_FILL_VEC_BPLUS(ptProng1); + CHECK_AND_FILL_VEC_BPLUS(impactParameter0); + CHECK_AND_FILL_VEC_BPLUS(impactParameter1); + CHECK_AND_FILL_VEC_BPLUS(impactParameterProduct); + CHECK_AND_FILL_VEC_BPLUS(chi2PCA); + CHECK_AND_FILL_VEC_BPLUS(decayLength); + CHECK_AND_FILL_VEC_BPLUS(decayLengthXY); + CHECK_AND_FILL_VEC_BPLUS(decayLengthNormalised); + CHECK_AND_FILL_VEC_BPLUS(decayLengthXYNormalised); + CHECK_AND_FILL_VEC_BPLUS(cpa); + CHECK_AND_FILL_VEC_BPLUS(cpaXY); + CHECK_AND_FILL_VEC_BPLUS(maxNormalisedDeltaIP); + // TPC PID variable + CHECK_AND_FILL_VEC_BPLUS_FULL(prong1, tpcNSigmaPi1, tpcNSigmaPi); + // TOF PID variable + CHECK_AND_FILL_VEC_BPLUS_FULL(prong1, tofNSigmaPi1, tofNSigmaPi); + // Combined PID variables + CHECK_AND_FILL_VEC_BPLUS_FUNC(prong1, tpcTofNSigmaPi1, o2::pid_tpc_tof_utils::getTpcTofNSigmaPi1); + } + } + } + + return inputFeatures; + } + + protected: + /// Method to fill the map of available input features + void setAvailableInputFeatures() + { + MlResponse::mAvailableInputFeatures = { + FILL_MAP_BPLUS(ptProng0), + FILL_MAP_BPLUS(ptProng1), + FILL_MAP_BPLUS(impactParameter0), + FILL_MAP_BPLUS(impactParameter1), + FILL_MAP_BPLUS(impactParameterProduct), + FILL_MAP_BPLUS(chi2PCA), + FILL_MAP_BPLUS(decayLength), + FILL_MAP_BPLUS(decayLengthXY), + FILL_MAP_BPLUS(decayLengthNormalised), + FILL_MAP_BPLUS(decayLengthXYNormalised), + FILL_MAP_BPLUS(cpa), + FILL_MAP_BPLUS(cpaXY), + FILL_MAP_BPLUS(maxNormalisedDeltaIP), + FILL_MAP_BPLUS(prong0MlScoreBkg), + FILL_MAP_BPLUS(prong0MlScorePrompt), + FILL_MAP_BPLUS(prong0MlScoreNonprompt), + // TPC PID variable + FILL_MAP_BPLUS(tpcNSigmaPi1), + // TOF PID variable + FILL_MAP_BPLUS(tofNSigmaPi1), + // Combined PID variable + FILL_MAP_BPLUS(tpcTofNSigmaPi1)}; + } +}; + +} // namespace o2::analysis + +#undef FILL_MAP_BPLUS +#undef CHECK_AND_FILL_VEC_BPLUS_FULL +#undef CHECK_AND_FILL_VEC_BPLUS_FUNC +#undef CHECK_AND_FILL_VEC_BPLUS + +#endif // PWGHF_CORE_HFMLRESPONSEBPLUSTOD0PI_H_ diff --git a/PWGHF/D2H/DataModel/ReducedDataModel.h b/PWGHF/D2H/DataModel/ReducedDataModel.h index f268e2e90ed..e4dc8c30a38 100644 --- a/PWGHF/D2H/DataModel/ReducedDataModel.h +++ b/PWGHF/D2H/DataModel/ReducedDataModel.h @@ -238,11 +238,20 @@ namespace hf_cand_bplus_reduced { DECLARE_SOA_INDEX_COLUMN_FULL(Prong0, prong0, int, HfRed2Prongs, "_0"); //! Prong0 index DECLARE_SOA_INDEX_COLUMN_FULL(Prong1, prong1, int, HfRedTrackBases, "_1"); //! Prong1 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_bplus_reduced DECLARE_SOA_TABLE(HfRedBplusProngs, "AOD", "HFREDBPPRONG", hf_cand_bplus_reduced::Prong0Id, hf_cand_bplus_reduced::Prong1Id); +DECLARE_SOA_TABLE(HfRedBplusD0Mls, "AOD", "HFREDBPLUSD0ML", //! Table with ML scores for the D+ daughter + hf_cand_bplus_reduced::Prong0MlScoreBkg, + hf_cand_bplus_reduced::Prong0MlScorePrompt, + hf_cand_bplus_reduced::Prong0MlScoreNonprompt, + o2::soa::Marker<1>); + using HfRedCandBplus = soa::Join; namespace hf_b0_mc @@ -340,6 +349,11 @@ DECLARE_SOA_COLUMN(EtaProng0, etaProng0, float); //! Pseudorapidity of the track DECLARE_SOA_COLUMN(PtProng1, ptProng1, float); //! Transverse momentum of the track's prong1 in GeV/c DECLARE_SOA_COLUMN(YProng1, yProng1, float); //! Rapidity of the track's prong1 DECLARE_SOA_COLUMN(EtaProng1, etaProng1, float); //! Pseudorapidity of the track's prong1 + +DECLARE_SOA_COLUMN(PdgCodeBeautyMother, pdgCodeBeautyMother, int); //! Pdg code of beauty mother +DECLARE_SOA_COLUMN(PdgCodeProng0, pdgCodeProng0, int); //! Pdg code of prong0 +DECLARE_SOA_COLUMN(PdgCodeProng1, pdgCodeProng1, int); //! Pdg code of prong1 +DECLARE_SOA_COLUMN(PdgCodeProng2, pdgCodeProng2, int); //! Pdg code of prong2 } // namespace hf_bplus_mc // table with results of reconstruction level MC matching @@ -350,12 +364,27 @@ DECLARE_SOA_TABLE(HfMcRecRedD0Pis, "AOD", "HFMCRECREDD0PI", //! Table with recon hf_cand_bplus::DebugMcRec, hf_bplus_mc::PtMother); +// DECLARE_SOA_EXTENDED_TABLE_USER(ExTable, Tracks, "EXTABLE", +DECLARE_SOA_TABLE(HfMcCheckD0Pis, "AOD", "HFMCCHECKD0PI", //! Table with reconstructed MC information on D0Pi(<-B0) pairs for MC checks in reduced workflow + hf_bplus_mc::PdgCodeBeautyMother, + hf_bplus_mc::PdgCodeProng0, + hf_bplus_mc::PdgCodeProng1, + hf_bplus_mc::PdgCodeProng2, + o2::soa::Marker<1>); + // Table with same size as HFCANDBPLUS DECLARE_SOA_TABLE(HfMcRecRedBps, "AOD", "HFMCRECREDBP", //! Reconstruction-level MC information on B+ candidates for reduced workflow hf_cand_bplus::FlagMcMatchRec, hf_cand_bplus::DebugMcRec, hf_bplus_mc::PtMother); +DECLARE_SOA_TABLE(HfMcCheckBps, "AOD", "HFMCCHECKBP", //! Table with reconstructed MC information on B+ candidates for MC checks in reduced workflow + hf_bplus_mc::PdgCodeBeautyMother, + hf_bplus_mc::PdgCodeProng0, + hf_bplus_mc::PdgCodeProng1, + hf_bplus_mc::PdgCodeProng2, + o2::soa::Marker<2>); + DECLARE_SOA_TABLE(HfMcGenRedBps, "AOD", "HFMCGENREDBP", //! Generation-level MC information on B+ candidates for reduced workflow hf_cand_bplus::FlagMcMatchGen, hf_bplus_mc::PtTrack, diff --git a/PWGHF/D2H/TableProducer/CMakeLists.txt b/PWGHF/D2H/TableProducer/CMakeLists.txt index 19cf1d44e5f..6f139aacc19 100644 --- a/PWGHF/D2H/TableProducer/CMakeLists.txt +++ b/PWGHF/D2H/TableProducer/CMakeLists.txt @@ -35,7 +35,7 @@ o2physics_add_dpl_workflow(candidate-selector-b0-to-d-pi-reduced o2physics_add_dpl_workflow(candidate-selector-bplus-to-d0-pi-reduced SOURCES candidateSelectorBplusToD0PiReduced.cxx - PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2Physics::MLCore COMPONENT_NAME Analysis) # Data creators diff --git a/PWGHF/D2H/TableProducer/candidateCreatorB0Reduced.cxx b/PWGHF/D2H/TableProducer/candidateCreatorB0Reduced.cxx index c12bd8da3b7..d81347d4eea 100644 --- a/PWGHF/D2H/TableProducer/candidateCreatorB0Reduced.cxx +++ b/PWGHF/D2H/TableProducer/candidateCreatorB0Reduced.cxx @@ -193,8 +193,7 @@ struct HfCandidateCreatorB0Reduced { pVecD[0], pVecD[1], pVecD[2], pVecPion[0], pVecPion[1], pVecPion[2], dcaD.getY(), dcaPion.getY(), - std::sqrt(dcaD.getSigmaY2()), std::sqrt(dcaPion.getSigmaY2()), - BIT(hf_cand_b0::DecayType::B0ToDPi)); + std::sqrt(dcaD.getSigmaY2()), std::sqrt(dcaPion.getSigmaY2())); rowCandidateProngs(candD.globalIndex(), trackPion.globalIndex()); diff --git a/PWGHF/D2H/TableProducer/candidateCreatorBplusReduced.cxx b/PWGHF/D2H/TableProducer/candidateCreatorBplusReduced.cxx index 8bdb7ed36d5..7f7fc71bd6e 100644 --- a/PWGHF/D2H/TableProducer/candidateCreatorBplusReduced.cxx +++ b/PWGHF/D2H/TableProducer/candidateCreatorBplusReduced.cxx @@ -37,8 +37,9 @@ using namespace o2::hf_trkcandsel; /// Reconstruction of B+ candidates struct HfCandidateCreatorBplusReduced { - Produces rowCandidateBase; // table defined in CandidateReconstructionTables.h - Produces rowCandidateProngs; // table defined in ReducedDataModel.h + Produces rowCandidateBase; // table defined in CandidateReconstructionTables.h + Produces rowCandidateProngs; // table defined in ReducedDataModel.h + Produces rowCandidateDmesMlScores; // table defined in ReducedDataModel.h // vertexing Configurable propagateToPCA{"propagateToPCA", true, "create tracks version propagated to PCA"}; @@ -61,6 +62,7 @@ struct HfCandidateCreatorBplusReduced { using HfRedCollisionsWithExtras = soa::Join; Preslice> candsDPerCollision = hf_track_index_reduced::hfRedCollisionId; + Preslice> candsDWithMlPerCollision = hf_track_index_reduced::hfRedCollisionId; Preslice> tracksPionPerCollision = hf_track_index_reduced::hfRedCollisionId; std::shared_ptr hCandidates; @@ -68,6 +70,11 @@ struct HfCandidateCreatorBplusReduced { void init(InitContext const&) { + std::array doprocess{doprocessData, doprocessDataWithDmesMl}; + if ((std::accumulate(doprocess.begin(), doprocess.end(), 0)) != 1) { + LOGP(fatal, "Only one process function for data should be enabled at a time."); + } + // invariant-mass window cut massPi = o2::constants::physics::MassPiPlus; massD0 = o2::constants::physics::MassD0; @@ -93,11 +100,112 @@ struct HfCandidateCreatorBplusReduced { setLabelHistoCands(hCandidates); } - void process(HfRedCollisionsWithExtras const& collisions, - soa::Join const& candsD, - soa::Join const& tracksPion, - aod::HfOrigColCounts const& collisionsCounter, - aod::HfCandBpConfigs const& configs) + /// Main function to perform B+ candidate creation + /// \param withDmesMl is the flag to use the table with ML scores for the D0 daughter (only possible if present in the derived data) + /// \param collision the collision + /// \param candsDThisColl B+ candidates in this collision + /// \param tracksPionThisCollision pion tracks in this collision + /// \param invMass2D0PiMin minimum B+ invariant-mass + /// \param invMass2D0PiMax maximum B+ invariant-mass + template + void runCandidateCreation(Coll const& collision, + Cands const& candsDThisColl, + Pions const& tracksPionThisCollision, + const float& invMass2D0PiMin, + const float& invMass2D0PiMax) + { + auto primaryVertex = getPrimaryVertex(collision); + auto covMatrixPV = primaryVertex.getCov(); + + // Set the magnetic field from ccdb + bz = collision.bz(); + df2.setBz(bz); + + for (const auto& candD0 : candsDThisColl) { + auto trackParCovD = getTrackParCov(candD0); + std::array pVecD0 = candD0.pVector(); + + for (const auto& trackPion : tracksPionThisCollision) { + auto trackParCovPi = getTrackParCov(trackPion); + std::array pVecPion = trackPion.pVector(); + + // 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; + } + // --------------------------------- + // reconstruct the 2-prong B+ vertex + hCandidates->Fill(SVFitting::BeforeFit); + try { + if (df2.process(trackParCovD, trackParCovPi) == 0) { + continue; + } + } catch (const std::runtime_error& error) { + LOG(info) << "Run time error found: " << error.what() << ". DCFitterN cannot work, skipping the candidate."; + hCandidates->Fill(SVFitting::Fail); + continue; + } + hCandidates->Fill(SVFitting::FitOk); + // D0Pi passed B+ reconstruction + + // calculate relevant properties + const auto& secondaryVertexBplus = df2.getPCACandidate(); + auto chi2PCA = df2.getChi2AtPCACandidate(); + auto covMatrixPCA = df2.calcPCACovMatrixFlat(); + registry.fill(HIST("hCovSVXX"), covMatrixPCA[0]); + registry.fill(HIST("hCovPVXX"), covMatrixPV[0]); + + // propagate D0 and Pi to the B+ vertex + df2.propagateTracksToVertex(); + // track.getPxPyPzGlo(pVec) modifies pVec of track + df2.getTrack(0).getPxPyPzGlo(pVecD0); // momentum of D0 at the B+ vertex + df2.getTrack(1).getPxPyPzGlo(pVecPion); // momentum of Pi at the B+ vertex + + registry.fill(HIST("hMassBplusToD0Pi"), std::sqrt(invMass2D0Pi)); + + // compute impact parameters of D0 and Pi + o2::dataformats::DCA dcaD0; + o2::dataformats::DCA dcaPion; + trackParCovD.propagateToDCA(primaryVertex, bz, &dcaD0); + trackParCovPi.propagateToDCA(primaryVertex, bz, &dcaPion); + + // get uncertainty of the decay length + double phi, theta; + // getPointDirection modifies phi and theta + getPointDirection(std::array{collision.posX(), collision.posY(), collision.posZ()}, secondaryVertexBplus, phi, theta); + auto errorDecayLength = std::sqrt(getRotatedCovMatrixXX(covMatrixPV, phi, theta) + getRotatedCovMatrixXX(covMatrixPCA, phi, theta)); + auto errorDecayLengthXY = std::sqrt(getRotatedCovMatrixXX(covMatrixPV, phi, 0.) + getRotatedCovMatrixXX(covMatrixPCA, phi, 0.)); + + // fill the candidate table for the B+ here: + rowCandidateBase(collision.globalIndex(), + collision.posX(), collision.posY(), collision.posZ(), + secondaryVertexBplus[0], secondaryVertexBplus[1], secondaryVertexBplus[2], + errorDecayLength, errorDecayLengthXY, + chi2PCA, + pVecD0[0], pVecD0[1], pVecD0[2], + pVecPion[0], pVecPion[1], pVecPion[2], + dcaD0.getY(), dcaPion.getY(), + std::sqrt(dcaD0.getSigmaY2()), std::sqrt(dcaPion.getSigmaY2())); + + rowCandidateProngs(candD0.globalIndex(), trackPion.globalIndex()); + + if constexpr (withDmesMl) { + if (trackPion.signed1Pt() < 0) { + rowCandidateDmesMlScores(candD0.mlScoreBkgMassHypo0(), candD0.mlScorePromptMassHypo0(), candD0.mlScoreNonpromptMassHypo0()); + } else { + rowCandidateDmesMlScores(candD0.mlScoreBkgMassHypo1(), candD0.mlScorePromptMassHypo1(), candD0.mlScoreNonpromptMassHypo1()); + } + } + } // pi loop + } // D0 loop + } // end runCandidateCreation + + void processData(HfRedCollisionsWithExtras const& collisions, + soa::Join const& candsD, + soa::Join const& tracksPion, + aod::HfOrigColCounts const& collisionsCounter, + aod::HfCandBpConfigs const& configs) { // D0Pi invariant-mass window cut for (const auto& config : configs) { @@ -116,105 +224,66 @@ struct HfCandidateCreatorBplusReduced { for (const auto& collision : collisions) { auto thisCollId = collision.globalIndex(); - auto primaryVertex = getPrimaryVertex(collision); - auto covMatrixPV = primaryVertex.getCov(); - + auto candsDThisColl = candsD.sliceBy(candsDPerCollision, thisCollId); + auto tracksPionThisCollision = tracksPion.sliceBy(tracksPionPerCollision, thisCollId); + runCandidateCreation(collision, candsDThisColl, tracksPionThisCollision, invMass2D0PiMin, invMass2D0PiMax); if (ncol % 10000 == 0) { LOG(debug) << ncol << " collisions parsed"; } ncol++; + } + } // processData - // Set the magnetic field from ccdb - bz = collision.bz(); - df2.setBz(bz); + PROCESS_SWITCH(HfCandidateCreatorBplusReduced, processData, "Process data without any ML score", true); + void processDataWithDmesMl(HfRedCollisionsWithExtras const& collisions, + soa::Join const& candsD, + soa::Join const& tracksPion, + aod::HfOrigColCounts const& collisionsCounter, + aod::HfCandBpConfigs const& configs) + { + // D0Pi invariant-mass window cut + for (const auto& config : configs) { + myInvMassWindowD0Pi = config.myInvMassWindowD0Pi(); + } + // invMassWindowD0PiTolerance is used to apply a slightly tighter cut than in D0Pi pair preselection + // to avoid accepting D0Pi pairs that were not formed in D0Pi pair creator + float invMass2D0PiMin = (massBplus - myInvMassWindowD0Pi + invMassWindowD0PiTolerance) * (massBplus - myInvMassWindowD0Pi + invMassWindowD0PiTolerance); + float invMass2D0PiMax = (massBplus + myInvMassWindowD0Pi - invMassWindowD0PiTolerance) * (massBplus + myInvMassWindowD0Pi - invMassWindowD0PiTolerance); + + for (const auto& collisionCounter : collisionsCounter) { + registry.fill(HIST("hEvents"), 1, collisionCounter.originalCollisionCount()); + } + + static int ncol = 0; + for (const auto& collision : collisions) { + auto thisCollId = collision.globalIndex(); auto candsDThisColl = candsD.sliceBy(candsDPerCollision, thisCollId); - for (const auto& candD0 : candsDThisColl) { - auto trackParCovD = getTrackParCov(candD0); - std::array pVecD0 = candD0.pVector(); - - auto tracksPionThisCollision = tracksPion.sliceBy(tracksPionPerCollision, thisCollId); - for (const auto& trackPion : tracksPionThisCollision) { - auto trackParCovPi = getTrackParCov(trackPion); - std::array pVecPion = trackPion.pVector(); - - // 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; - } - // --------------------------------- - // reconstruct the 2-prong B+ vertex - hCandidates->Fill(SVFitting::BeforeFit); - try { - if (df2.process(trackParCovD, trackParCovPi) == 0) { - continue; - } - } catch (const std::runtime_error& error) { - LOG(info) << "Run time error found: " << error.what() << ". DCFitterN cannot work, skipping the candidate."; - hCandidates->Fill(SVFitting::Fail); - continue; - } - hCandidates->Fill(SVFitting::FitOk); - - // D0Pi passed B+ reconstruction - - // calculate relevant properties - const auto& secondaryVertexBplus = df2.getPCACandidate(); - auto chi2PCA = df2.getChi2AtPCACandidate(); - auto covMatrixPCA = df2.calcPCACovMatrixFlat(); - registry.fill(HIST("hCovSVXX"), covMatrixPCA[0]); - registry.fill(HIST("hCovPVXX"), covMatrixPV[0]); - - // propagate D0 and Pi to the B+ vertex - df2.propagateTracksToVertex(); - // track.getPxPyPzGlo(pVec) modifies pVec of track - df2.getTrack(0).getPxPyPzGlo(pVecD0); // momentum of D0 at the B+ vertex - df2.getTrack(1).getPxPyPzGlo(pVecPion); // momentum of Pi at the B+ vertex - - registry.fill(HIST("hMassBplusToD0Pi"), std::sqrt(invMass2D0Pi)); - - // compute impact parameters of D0 and Pi - o2::dataformats::DCA dcaD0; - o2::dataformats::DCA dcaPion; - trackParCovD.propagateToDCA(primaryVertex, bz, &dcaD0); - trackParCovPi.propagateToDCA(primaryVertex, bz, &dcaPion); - - // get uncertainty of the decay length - double phi, theta; - // getPointDirection modifies phi and theta - getPointDirection(std::array{collision.posX(), collision.posY(), collision.posZ()}, secondaryVertexBplus, phi, theta); - auto errorDecayLength = std::sqrt(getRotatedCovMatrixXX(covMatrixPV, phi, theta) + getRotatedCovMatrixXX(covMatrixPCA, phi, theta)); - auto errorDecayLengthXY = std::sqrt(getRotatedCovMatrixXX(covMatrixPV, phi, 0.) + getRotatedCovMatrixXX(covMatrixPCA, phi, 0.)); - - int hfFlag = BIT(hf_cand_bplus::DecayType::BplusToD0Pi); - - // fill the candidate table for the B+ here: - rowCandidateBase(thisCollId, - collision.posX(), collision.posY(), collision.posZ(), - secondaryVertexBplus[0], secondaryVertexBplus[1], secondaryVertexBplus[2], - errorDecayLength, errorDecayLengthXY, - chi2PCA, - pVecD0[0], pVecD0[1], pVecD0[2], - pVecPion[0], pVecPion[1], pVecPion[2], - dcaD0.getY(), dcaPion.getY(), - std::sqrt(dcaD0.getSigmaY2()), std::sqrt(dcaPion.getSigmaY2()), - hfFlag); - - rowCandidateProngs(candD0.globalIndex(), trackPion.globalIndex()); - } // pi loop - } // D0 loop - } // collision loop - } // process -}; // struct + auto tracksPionThisCollision = tracksPion.sliceBy(tracksPionPerCollision, thisCollId); + runCandidateCreation(collision, candsDThisColl, tracksPionThisCollision, invMass2D0PiMin, invMass2D0PiMax); + if (ncol % 10000 == 0) { + LOGP(debug, "collisions parsed {}", ncol); + } + ncol++; + } + } // processDataWithDmesMl + + PROCESS_SWITCH(HfCandidateCreatorBplusReduced, processDataWithDmesMl, "Process data with ML scores of D mesons", false); +}; // struct /// Extends the table base with expression columns and performs MC matching. struct HfCandidateCreatorBplusReducedExpressions { Spawns rowCandidateBPlus; Spawns rowTracksExt; Produces rowBplusMcRec; - - void processMc(HfMcRecRedD0Pis const& rowsD0PiMcRec, HfRedBplusProngs const& candsBplus) + Produces rowBplusMcCheck; + + /// Fill candidate information at MC reconstruction level + /// \param checkDecayTypeMc + /// \param rowsD0PiMcRec MC reco information on D0Pi pairs + /// \param candsBplus prong global indices of B+ candidates + template + void fillBplusMcRec(McRec const& rowsD0PiMcRec, HfRedBplusProngs const& candsBplus) { for (const auto& candBplus : candsBplus) { bool filledMcInfo{false}; @@ -224,14 +293,34 @@ struct HfCandidateCreatorBplusReducedExpressions { } rowBplusMcRec(rowD0PiMcRec.flagMcMatchRec(), rowD0PiMcRec.debugMcRec(), rowD0PiMcRec.ptMother()); filledMcInfo = true; + if constexpr (checkDecayTypeMc) { + rowBplusMcCheck(rowD0PiMcRec.pdgCodeBeautyMother(), + rowD0PiMcRec.pdgCodeProng0(), + rowD0PiMcRec.pdgCodeProng1(), + rowD0PiMcRec.pdgCodeProng2()); + } break; } - if (!filledMcInfo) { // protection to get same size tables in case something went wrong: we created a candidate that was not preselected in the D-Pi creator + if (!filledMcInfo) { // protection to get same size tables in case something went wrong: we created a candidate that was not preselected in the D0-Pi creator rowBplusMcRec(0, -1, -1.f); + if constexpr (checkDecayTypeMc) { + rowBplusMcCheck(-1, -1, -1, -1); + } } } } + + void processMc(HfMcRecRedD0Pis const& rowsD0PiMcRec, HfRedBplusProngs const& candsBplus) + { + fillBplusMcRec(rowsD0PiMcRec, candsBplus); + } PROCESS_SWITCH(HfCandidateCreatorBplusReducedExpressions, processMc, "Process MC", false); + + void processMcWithDecayTypeCheck(soa::Join const& rowsD0PiMcRec, HfRedBplusProngs const& candsBplus) + { + fillBplusMcRec(rowsD0PiMcRec, candsBplus); + } + PROCESS_SWITCH(HfCandidateCreatorBplusReducedExpressions, processMcWithDecayTypeCheck, "Process MC with decay type checks", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) diff --git a/PWGHF/D2H/TableProducer/candidateSelectorB0ToDPiReduced.cxx b/PWGHF/D2H/TableProducer/candidateSelectorB0ToDPiReduced.cxx index 3653471dbbf..5bf352517e4 100644 --- a/PWGHF/D2H/TableProducer/candidateSelectorB0ToDPiReduced.cxx +++ b/PWGHF/D2H/TableProducer/candidateSelectorB0ToDPiReduced.cxx @@ -39,7 +39,7 @@ struct HfCandidateSelectorB0ToDPiReduced { Configurable ptCandMin{"ptCandMin", 0., "Lower bound of candidate pT"}; Configurable ptCandMax{"ptCandMax", 50., "Upper bound of candidate pT"}; // Enable PID - Configurable usePionPid{"usePionPid", 1, "Switch for PID selection for the bachelor pion (0: none, 1: TPC or TOF, 2: TPC and TOF)"}; + Configurable pionPidMethod{"pionPidMethod", 1, "PID selection method for the bachelor pion (0: none, 1: TPC or TOF, 2: TPC and TOF)"}; Configurable acceptPIDNotApplicable{"acceptPIDNotApplicable", true, "Switch to accept Status::NotApplicable [(NotApplicable for one detector) and (NotApplicable or Conditional for the other)] in PID selection"}; // TPC PID Configurable ptPidTpcMin{"ptPidTpcMin", 0.15, "Lower bound of track pT for TPC PID"}; @@ -94,11 +94,11 @@ struct HfCandidateSelectorB0ToDPiReduced { LOGP(fatal, "Only one process function for data should be enabled at a time."); } - if (usePionPid < 0 || usePionPid > 2) { + if (pionPidMethod < 0 || pionPidMethod > 2) { LOGP(fatal, "Invalid PID option in configurable, please set 0 (no PID), 1 (TPC or TOF), or 2 (TPC and TOF)"); } - if (usePionPid) { + if (pionPidMethod) { selectorPion.setRangePtTpc(ptPidTpcMin, ptPidTpcMax); selectorPion.setRangeNSigmaTpc(-nSigmaTpcMax, nSigmaTpcMax); selectorPion.setRangeNSigmaTpcCondTof(-nSigmaTpcCombinedMax, nSigmaTpcCombinedMax); @@ -153,18 +153,6 @@ struct HfCandidateSelectorB0ToDPiReduced { int statusB0ToDPi = 0; auto ptCandB0 = hfCandB0.pt(); - // check if flagged as B0 → D π - if (!TESTBIT(hfCandB0.hfflag(), hf_cand_b0::DecayType::B0ToDPi)) { - hfSelB0ToDPiCandidate(statusB0ToDPi); - if (applyB0Ml) { - hfMlB0ToDPiCandidate(outputMlNotPreselected); - } - if (activateQA) { - registry.fill(HIST("hSelections"), 1, ptCandB0); - } - // LOGF(info, "B0 candidate selection failed at hfflag check"); - continue; - } SETBIT(statusB0ToDPi, SelectionStep::RecoSkims); // RecoSkims = 0 --> statusB0ToDPi = 1 if (activateQA) { registry.fill(HIST("hSelections"), 2 + SelectionStep::RecoSkims, ptCandB0); @@ -198,9 +186,9 @@ struct HfCandidateSelectorB0ToDPiReduced { // track-level PID selection auto trackPi = hfCandB0.template prong1_as(); - if (usePionPid) { + if (pionPidMethod) { int pidTrackPi{TrackSelectorPID::Status::NotApplicable}; - if (usePionPid == 1) { + if (pionPidMethod == 1) { pidTrackPi = selectorPion.statusTpcOrTof(trackPi); } else { pidTrackPi = selectorPion.statusTpcAndTof(trackPi); @@ -229,7 +217,7 @@ struct HfCandidateSelectorB0ToDPiReduced { hfSelB0ToDPiCandidate(statusB0ToDPi); continue; } - SETBIT(statusB0ToDPi, SelectionStep::RecoMl); // RecoML = 3 --> statusB0ToDPi = 15 if usePionPid, 11 otherwise + SETBIT(statusB0ToDPi, SelectionStep::RecoMl); // RecoML = 3 --> statusB0ToDPi = 15 if pionPidMethod, 11 otherwise if (activateQA) { registry.fill(HIST("hSelections"), 2 + SelectionStep::RecoMl, ptCandB0); } diff --git a/PWGHF/D2H/TableProducer/candidateSelectorBplusToD0PiReduced.cxx b/PWGHF/D2H/TableProducer/candidateSelectorBplusToD0PiReduced.cxx index 32ce92edb07..551444b3904 100644 --- a/PWGHF/D2H/TableProducer/candidateSelectorBplusToD0PiReduced.cxx +++ b/PWGHF/D2H/TableProducer/candidateSelectorBplusToD0PiReduced.cxx @@ -20,6 +20,7 @@ #include "Common/Core/TrackSelectorPID.h" #include "PWGHF/Core/HfHelper.h" +#include "PWGHF/Core/HfMlResponseBplusToD0Pi.h" #include "PWGHF/Core/SelectorCuts.h" #include "PWGHF/DataModel/CandidateReconstructionTables.h" #include "PWGHF/DataModel/CandidateSelectionTables.h" @@ -32,11 +33,12 @@ using namespace o2::analysis; struct HfCandidateSelectorBplusToD0PiReduced { Produces hfSelBplusToD0PiCandidate; // table defined in CandidateSelectionTables.h + Produces hfMlBplusToD0PiCandidate; // table defined in CandidateSelectionTables.h Configurable ptCandMin{"ptCandMin", 0., "Lower bound of candidate pT"}; Configurable ptCandMax{"ptCandMax", 50., "Upper bound of candidate pT"}; // Enable PID - Configurable usePid{"usePid", true, "Switch for PID selection at track level"}; + Configurable pionPidMethod{"pionPidMethod", 1, "PID selection method for the bachelor pion (0: none, 1: TPC or TOF, 2: TPC and TOF)"}; Configurable acceptPIDNotApplicable{"acceptPIDNotApplicable", true, "Switch to accept Status::NotApplicable [(NotApplicable for one detector) and (NotApplicable or Conditional for the other)] in PID selection"}; // TPC PID Configurable ptPidTpcMin{"ptPidTpcMin", 0.15, "Lower bound of track pT for TPC PID"}; @@ -51,14 +53,33 @@ struct HfCandidateSelectorBplusToD0PiReduced { // topological cuts Configurable> binsPt{"binsPt", std::vector{hf_cuts_bplus_to_d0_pi::vecBinsPt}, "pT bin limits"}; Configurable> cuts{"cuts", {hf_cuts_bplus_to_d0_pi::cuts[0], hf_cuts_bplus_to_d0_pi::nBinsPt, hf_cuts_bplus_to_d0_pi::nCutVars, hf_cuts_bplus_to_d0_pi::labelsPt, hf_cuts_bplus_to_d0_pi::labelsCutVar}, "B+ candidate selection per pT bin"}; + // D0-meson ML cuts + Configurable> binsPtDmesMl{"binsPtDmesMl", std::vector{hf_cuts_ml::vecBinsPt}, "D0-meson pT bin limits for ML cuts"}; + Configurable> cutsDmesMl{"cutsDmesMl", {hf_cuts_ml::cuts[0], hf_cuts_ml::nBinsPt, hf_cuts_ml::nCutScores, hf_cuts_ml::labelsPt, hf_cuts_ml::labelsDmesCutScore}, "D0-meson ML cuts per pT bin"}; // QA switch Configurable activateQA{"activateQA", false, "Flag to enable QA histogram"}; - - bool selectionFlagDAndUsePidInSync = true; + // B+ ML inference + Configurable applyBplusMl{"applyBplusMl", false, "Flag to apply ML selections"}; + Configurable> binsPtBpMl{"binsPtBpMl", std::vector{hf_cuts_ml::vecBinsPt}, "pT bin limits for ML application"}; + Configurable> cutDirBpMl{"cutDirBpMl", std::vector{hf_cuts_ml::vecCutDir}, "Whether to reject score values greater or smaller than the threshold"}; + Configurable> cutsBpMl{"cutsBpMl", {hf_cuts_ml::cuts[0], hf_cuts_ml::nBinsPt, hf_cuts_ml::nCutScores, hf_cuts_ml::labelsPt, hf_cuts_ml::labelsCutScore}, "ML selections per pT bin"}; + Configurable nClassesBpMl{"nClassesBpMl", (int8_t)hf_cuts_ml::nCutScores, "Number of classes in ML model"}; + Configurable> namesInputFeatures{"namesInputFeatures", std::vector{"feature1", "feature2"}, "Names of ML model input features"}; + // CCDB configuration + Configurable ccdbUrl{"ccdbUrl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; + Configurable> modelPathsCCDB{"modelPathsCCDB", std::vector{"path_ccdb/BDT_BPLUS/"}, "Paths of models on CCDB"}; + Configurable> onnxFileNames{"onnxFileNames", std::vector{"ModelHandler_onnx_BPLUSToD0Pi.onnx"}, "ONNX file names for each pT bin (if not from CCDB full path)"}; + Configurable timestampCCDB{"timestampCCDB", -1, "timestamp of the ONNX file for ML model used to query in CCDB"}; + Configurable loadModelsFromCCDB{"loadModelsFromCCDB", false, "Flag to enable or disable the loading of models from CCDB"}; // variable that will store the value of selectionFlagD (defined in dataCreatorD0PiReduced.cxx) int mySelectionFlagD0 = -1; int mySelectionFlagD0bar = -1; + o2::analysis::HfMlResponseBplusToD0Pi hfMlResponse; + float outputMlNotPreselected = -1.; + std::vector outputMl = {}; + o2::ccdb::CcdbApi ccdbApi; + HfHelper hfHelper; TrackSelectorPi selectorPion; @@ -68,7 +89,16 @@ struct HfCandidateSelectorBplusToD0PiReduced { void init(InitContext const&) { - if (usePid) { + std::array doprocess{doprocessSelection, doprocessSelectionWithDmesMl}; + if ((std::accumulate(doprocess.begin(), doprocess.end(), 0)) != 1) { + LOGP(fatal, "Only one process function for data should be enabled at a time."); + } + + if (pionPidMethod < 0 || pionPidMethod > 2) { + LOGP(fatal, "Invalid PID option in configurable, please set 0 (no PID), 1 (TPC or TOF), or 2 (TPC and TOF)"); + } + + if (pionPidMethod) { selectorPion.setRangePtTpc(ptPidTpcMin, ptPidTpcMax); selectorPion.setRangeNSigmaTpc(-nSigmaTpcMax, nSigmaTpcMax); selectorPion.setRangeNSigmaTpcCondTof(-nSigmaTpcCombinedMax, nSigmaTpcCombinedMax); @@ -90,40 +120,40 @@ struct HfCandidateSelectorBplusToD0PiReduced { registry.get(HIST("hSelections"))->GetXaxis()->SetBinLabel(iBin + 1, labels[iBin].data()); } } + + if (applyBplusMl) { + hfMlResponse.configure(binsPtBpMl, cutsBpMl, cutDirBpMl, nClassesBpMl); + if (loadModelsFromCCDB) { + ccdbApi.init(ccdbUrl); + hfMlResponse.setModelPathsCCDB(onnxFileNames, ccdbApi, modelPathsCCDB, timestampCCDB); + } else { + hfMlResponse.setModelPathsLocal(onnxFileNames); + } + hfMlResponse.cacheInputFeaturesIndices(namesInputFeatures); + hfMlResponse.init(); + } } - void process(HfRedCandBplus const& hfCandBs, - TracksPion const&, - HfCandBpConfigs const& configs) + /// Main function to perform B+ candidate creation + /// \param withDmesMl is the flag to use the table with ML scores for the D- daughter (only possible if present in the derived data) + /// \param hfCandsBp B+ candidates + /// \param pionTracks pion tracks + /// \param configs config inherited from the D0pi data creator + template + void runSelection(Cands const& hfCandsBp, + TracksPion const& pionTracks, + HfCandBpConfigs const& configs) { - // get DplusPi creator configurable + // get D0Pi creator configurable for (const auto& config : configs) { mySelectionFlagD0 = config.mySelectionFlagD0(); mySelectionFlagD0bar = config.mySelectionFlagD0bar(); - - if ((usePid && !mySelectionFlagD0) || (usePid && !mySelectionFlagD0bar)) { - selectionFlagDAndUsePidInSync = false; - LOG(warning) << "PID selections required on B+ daughters (usePid=true) but no PID selections on D candidates were required a priori."; - } - if ((!usePid && mySelectionFlagD0) || (!usePid && mySelectionFlagD0bar)) { - selectionFlagDAndUsePidInSync = false; - LOG(warning) << "No PID selections required on Bp daughters (usePid=false) but PID selections on D candidates were required a priori."; - } } - for (const auto& hfCandBp : hfCandBs) { + for (const auto& hfCandBp : hfCandsBp) { int statusBplus = 0; auto ptCandBplus = hfCandBp.pt(); - // check if flagged as B+ → D π - if (!TESTBIT(hfCandBp.hfflag(), hf_cand_bplus::DecayType::BplusToD0Pi)) { - hfSelBplusToD0PiCandidate(statusBplus); - if (activateQA) { - registry.fill(HIST("hSelections"), 1, ptCandBplus); - } - // LOGF(info, "B+ candidate selection failed at hfflag check"); - continue; - } SETBIT(statusBplus, SelectionStep::RecoSkims); // RecoSkims = 0 --> statusBplus = 1 if (activateQA) { registry.fill(HIST("hSelections"), 2 + SelectionStep::RecoSkims, ptCandBplus); @@ -132,26 +162,44 @@ struct HfCandidateSelectorBplusToD0PiReduced { // topological cuts if (!hfHelper.selectionBplusToD0PiTopol(hfCandBp, cuts, binsPt)) { hfSelBplusToD0PiCandidate(statusBplus); + if (applyBplusMl) { + hfMlBplusToD0PiCandidate(outputMlNotPreselected); + } // LOGF(info, "B+ candidate selection failed at topology selection"); continue; } + + if constexpr (withDmesMl) { // we include it in the topological selections + if (!hfHelper.selectionDmesMlScoresForB(hfCandBp, cutsDmesMl, binsPtDmesMl)) { + hfSelBplusToD0PiCandidate(statusBplus); + if (applyBplusMl) { + hfMlBplusToD0PiCandidate(outputMlNotPreselected); + } + // LOGF(info, "B+ candidate selection failed at D0-meson ML selection"); + continue; + } + } + SETBIT(statusBplus, SelectionStep::RecoTopol); // RecoTopol = 1 --> statusBplus = 3 if (activateQA) { registry.fill(HIST("hSelections"), 2 + SelectionStep::RecoTopol, ptCandBplus); } - // checking if selectionFlagD and usePid are in sync - if (!selectionFlagDAndUsePidInSync) { - hfSelBplusToD0PiCandidate(statusBplus); - continue; - } // track-level PID selection - if (usePid) { - auto trackPi = hfCandBp.prong1_as(); - int pidTrackPi = selectorPion.statusTpcAndTof(trackPi); + auto trackPi = hfCandBp.template prong1_as(); + if (pionPidMethod) { + int pidTrackPi{TrackSelectorPID::Status::NotApplicable}; + if (pionPidMethod == 1) { + pidTrackPi = selectorPion.statusTpcOrTof(trackPi); + } else { + pidTrackPi = selectorPion.statusTpcAndTof(trackPi); + } if (!hfHelper.selectionBplusToD0PiPid(pidTrackPi, acceptPIDNotApplicable.value)) { // LOGF(info, "B+ candidate selection failed at PID selection"); hfSelBplusToD0PiCandidate(statusBplus); + if (applyBplusMl) { + hfMlBplusToD0PiCandidate(outputMlNotPreselected); + } continue; } SETBIT(statusBplus, SelectionStep::RecoPID); // RecoPID = 2 --> statusBplus = 7 @@ -159,10 +207,44 @@ struct HfCandidateSelectorBplusToD0PiReduced { registry.fill(HIST("hSelections"), 2 + SelectionStep::RecoPID, ptCandBplus); } } + if (applyBplusMl) { + // B+ ML selections + std::vector inputFeatures = hfMlResponse.getInputFeatures(hfCandBp, trackPi); + bool isSelectedMl = hfMlResponse.isSelectedMl(inputFeatures, ptCandBplus, outputMl); + hfMlBplusToD0PiCandidate(outputMl[1]); // storing ML score for signal class + + if (!isSelectedMl) { + hfSelBplusToD0PiCandidate(statusBplus); + continue; + } + SETBIT(statusBplus, SelectionStep::RecoMl); // RecoML = 3 --> statusBplus = 15 if pionPidMethod, 11 otherwise + if (activateQA) { + registry.fill(HIST("hSelections"), 2 + SelectionStep::RecoMl, ptCandBplus); + } + } + hfSelBplusToD0PiCandidate(statusBplus); // LOGF(info, "B+ candidate selection passed all selections"); } } + + void processSelection(HfRedCandBplus const& hfCandsBp, + TracksPion const& pionTracks, + HfCandBpConfigs const& configs) + { + runSelection(hfCandsBp, pionTracks, configs); + } // processSelection + + PROCESS_SWITCH(HfCandidateSelectorBplusToD0PiReduced, processSelection, "Process selection without ML scores of D mesons", true); + + void processSelectionWithDmesMl(soa::Join const& hfCandsBp, + TracksPion const& pionTracks, + HfCandBpConfigs const& configs) + { + runSelection(hfCandsBp, pionTracks, configs); + } // processSelectionWithDmesMl + + PROCESS_SWITCH(HfCandidateSelectorBplusToD0PiReduced, processSelectionWithDmesMl, "Process selection with ML scores of D mesons", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) diff --git a/PWGHF/D2H/TableProducer/dataCreatorCharmHadPiReduced.cxx b/PWGHF/D2H/TableProducer/dataCreatorCharmHadPiReduced.cxx index b4b9f12eb23..ffa5c7160bf 100644 --- a/PWGHF/D2H/TableProducer/dataCreatorCharmHadPiReduced.cxx +++ b/PWGHF/D2H/TableProducer/dataCreatorCharmHadPiReduced.cxx @@ -81,6 +81,7 @@ struct HfDataCreatorCharmHadPiReduced { Produces rowCandidateConfigBplus; Produces rowHfD0PiMcRecReduced; + Produces rowHfD0PiMcCheckReduced; Produces rowHfBpMcGenReduced; // vertexing @@ -409,14 +410,13 @@ struct HfDataCreatorCharmHadPiReduced { // B+ → D0(bar) π+ → (K+ π-) π+ auto indexRec = RecoDecay::getMatchedMCRec(particlesMc, std::array{vecDaughtersB[0], vecDaughtersB[1], vecDaughtersB[2]}, Pdg::kBPlus, std::array{+kPiPlus, +kKPlus, -kPiPlus}, true, &sign, 2); if (indexRec > -1) { - // D0bar → K+ π- // Printf("Checking D0bar → K+ π-"); indexRec = RecoDecay::getMatchedMCRec(particlesMc, std::array{vecDaughtersB[0], vecDaughtersB[1]}, Pdg::kD0, std::array{+kPiPlus, -kKPlus}, true, &sign, 1); if (indexRec > -1) { flag = sign * BIT(hf_cand_bplus::DecayType::BplusToD0Pi); } else { debug = 1; - LOGF(debug, "B0 decays in the expected final state but the condition on the intermediate state is not fulfilled"); + LOGF(debug, "B+ decays in the expected final state but the condition on the intermediate state is not fulfilled"); } auto indexMother = RecoDecay::getMother(particlesMc, vecDaughtersB.back().template mcParticle_as(), Pdg::kBPlus, true); @@ -426,6 +426,39 @@ struct HfDataCreatorCharmHadPiReduced { } } rowHfD0PiMcRecReduced(indexHfCandCharm, selectedTracksPion[vecDaughtersB.back().globalIndex()], flag, debug, motherPt); + + // additional checks for correlated backgrounds + if (checkDecayTypeMc) { + // Partly reconstructed decays, i.e. the 3 prongs have a common b-hadron ancestor + // convention: final state particles are prong0,1,2 + if (!flag) { + auto particleProng0 = vecDaughtersB[0].mcParticle(); + auto particleProng1 = vecDaughtersB[1].mcParticle(); + auto particleProng2 = vecDaughtersB[2].mcParticle(); + // b-hadron hypothesis + std::array bHadronMotherHypos = {Pdg::kBPlus, Pdg::kB0, Pdg::kBS}; + + for (const auto& bHadronMotherHypo : bHadronMotherHypos) { + int index0Mother = RecoDecay::getMother(particlesMc, particleProng0, bHadronMotherHypo, true); + int index1Mother = RecoDecay::getMother(particlesMc, particleProng1, bHadronMotherHypo, true); + int index2Mother = RecoDecay::getMother(particlesMc, particleProng2, bHadronMotherHypo, true); + + // look for common b-hadron ancestor + if (index0Mother > -1 && index1Mother > -1 && index2Mother > -1) { + if (index0Mother == index1Mother && index1Mother == index2Mother) { + flag = BIT(hf_cand_bplus::DecayTypeMc::PartlyRecoDecay); + pdgCodeBeautyMother = particlesMc.rawIteratorAt(index0Mother).pdgCode(); + pdgCodeCharmMother = 0; + pdgCodeProng0 = particleProng0.pdgCode(); + pdgCodeProng1 = particleProng1.pdgCode(); + pdgCodeProng2 = particleProng2.pdgCode(); + break; + } + } + } + } + rowHfD0PiMcCheckReduced(pdgCodeBeautyMother, pdgCodeProng0, pdgCodeProng1, pdgCodeProng2); + } } } diff --git a/PWGHF/D2H/Tasks/taskB0.cxx b/PWGHF/D2H/Tasks/taskB0.cxx index 29d59ed912e..43025a47f6b 100644 --- a/PWGHF/D2H/Tasks/taskB0.cxx +++ b/PWGHF/D2H/Tasks/taskB0.cxx @@ -157,9 +157,6 @@ struct HfTaskB0 { TracksWithSel const&) { for (const auto& candidate : candidates) { - if (!TESTBIT(candidate.hfflag(), hf_cand_b0::DecayType::B0ToDPi)) { - continue; - } if (yCandRecoMax >= 0. && std::abs(hfHelper.yB0(candidate)) > yCandRecoMax) { continue; } @@ -197,9 +194,6 @@ struct HfTaskB0 { { // MC rec for (const auto& candidate : candidates) { - if (!TESTBIT(candidate.hfflag(), hf_cand_b0::DecayType::B0ToDPi)) { - continue; - } if (yCandRecoMax >= 0. && std::abs(hfHelper.yB0(candidate)) > yCandRecoMax) { continue; } diff --git a/PWGHF/D2H/Tasks/taskB0Reduced.cxx b/PWGHF/D2H/Tasks/taskB0Reduced.cxx index 9b5a9c2d40a..5e38055f589 100644 --- a/PWGHF/D2H/Tasks/taskB0Reduced.cxx +++ b/PWGHF/D2H/Tasks/taskB0Reduced.cxx @@ -592,9 +592,6 @@ struct HfTaskB0Reduced { TracksPion const&) { for (const auto& candidate : candidates) { - if (!TESTBIT(candidate.hfflag(), hf_cand_b0::DecayType::B0ToDPi)) { - continue; - } if (yCandRecoMax >= 0. && std::abs(hfHelper.yB0(candidate)) > yCandRecoMax) { continue; } @@ -608,9 +605,6 @@ struct HfTaskB0Reduced { TracksPion const&) { for (const auto& candidate : candidates) { - if (!TESTBIT(candidate.hfflag(), hf_cand_b0::DecayType::B0ToDPi)) { - continue; - } if (yCandRecoMax >= 0. && std::abs(hfHelper.yB0(candidate)) > yCandRecoMax) { continue; } @@ -624,9 +618,6 @@ struct HfTaskB0Reduced { TracksPion const&) { for (const auto& candidate : candidates) { - if (!TESTBIT(candidate.hfflag(), hf_cand_b0::DecayType::B0ToDPi)) { - continue; - } if (yCandRecoMax >= 0. && std::abs(hfHelper.yB0(candidate)) > yCandRecoMax) { continue; } @@ -642,9 +633,6 @@ struct HfTaskB0Reduced { { // MC rec for (const auto& candidate : candidates) { - if (!TESTBIT(candidate.hfflag(), hf_cand_b0::DecayType::B0ToDPi)) { - continue; - } if (yCandRecoMax >= 0. && std::abs(hfHelper.yB0(candidate)) > yCandRecoMax) { continue; } @@ -665,9 +653,6 @@ struct HfTaskB0Reduced { { // MC rec for (const auto& candidate : candidates) { - if (!TESTBIT(candidate.hfflag(), hf_cand_b0::DecayType::B0ToDPi)) { - continue; - } if (yCandRecoMax >= 0. && std::abs(hfHelper.yB0(candidate)) > yCandRecoMax) { continue; } @@ -688,9 +673,6 @@ struct HfTaskB0Reduced { { // MC rec for (const auto& candidate : candidates) { - if (!TESTBIT(candidate.hfflag(), hf_cand_b0::DecayType::B0ToDPi)) { - continue; - } if (yCandRecoMax >= 0. && std::abs(hfHelper.yB0(candidate)) > yCandRecoMax) { continue; } @@ -711,9 +693,6 @@ struct HfTaskB0Reduced { { // MC rec for (const auto& candidate : candidates) { - if (!TESTBIT(candidate.hfflag(), hf_cand_b0::DecayType::B0ToDPi)) { - continue; - } if (yCandRecoMax >= 0. && std::abs(hfHelper.yB0(candidate)) > yCandRecoMax) { continue; } @@ -734,9 +713,6 @@ struct HfTaskB0Reduced { { // MC rec for (const auto& candidate : candidates) { - if (!TESTBIT(candidate.hfflag(), hf_cand_b0::DecayType::B0ToDPi)) { - continue; - } if (yCandRecoMax >= 0. && std::abs(hfHelper.yB0(candidate)) > yCandRecoMax) { continue; } diff --git a/PWGHF/D2H/Tasks/taskBplus.cxx b/PWGHF/D2H/Tasks/taskBplus.cxx index 6920054084c..9586ebaa424 100644 --- a/PWGHF/D2H/Tasks/taskBplus.cxx +++ b/PWGHF/D2H/Tasks/taskBplus.cxx @@ -164,9 +164,6 @@ struct HfTaskBplus { { for (const auto& candidate : selectedBPlusCandidates) { - if (!TESTBIT(candidate.hfflag(), hf_cand_bplus::DecayType::BplusToD0Pi)) { - continue; - } if (yCandRecoMax >= 0. && std::abs(hfHelper.yBplus(candidate)) > yCandRecoMax) { continue; } @@ -208,9 +205,6 @@ struct HfTaskBplus { { // MC rec for (const auto& candidate : selectedBPlusCandidatesMC) { - if (!TESTBIT(candidate.hfflag(), hf_cand_bplus::DecayType::BplusToD0Pi)) { - continue; - } if (yCandRecoMax >= 0. && std::abs(hfHelper.yBplus(candidate)) > yCandRecoMax) { continue; } @@ -294,7 +288,7 @@ struct HfTaskBplus { registry.fill(HIST("hPtGenWithRapidityBelowHalf"), ptParticle); } - // reject B0 daughters that are not in geometrical acceptance + // reject B+ daughters that are not in geometrical acceptance if (!isProngInAcceptance(etaProngs[0], ptProngs[0]) || !isProngInAcceptance(etaProngs[1], ptProngs[1])) { continue; } diff --git a/PWGHF/D2H/Tasks/taskBplusReduced.cxx b/PWGHF/D2H/Tasks/taskBplusReduced.cxx index 6f1b775c240..bedaf94c221 100644 --- a/PWGHF/D2H/Tasks/taskBplusReduced.cxx +++ b/PWGHF/D2H/Tasks/taskBplusReduced.cxx @@ -17,6 +17,7 @@ #include "Framework/AnalysisTask.h" #include "Framework/HistogramRegistry.h" #include "Framework/runDataProcessing.h" +#include "Common/Core/RecoDecay.h" #include "PWGHF/Core/HfHelper.h" #include "PWGHF/Core/SelectorCuts.h" @@ -30,8 +31,81 @@ using namespace o2::analysis; using namespace o2::framework; using namespace o2::framework::expressions; +namespace o2::aod +{ +namespace hf_cand_bplus_lite +{ +DECLARE_SOA_COLUMN(PtProng0, ptProng0, float); //! Transverse momentum of prong0 (GeV/c) +DECLARE_SOA_COLUMN(PtProng1, ptProng1, float); //! Transverse momentum of prong1 (GeV/c) +DECLARE_SOA_COLUMN(MProng0, mProng0, float); //! Invariant mass of prong0 (GeV/c) +DECLARE_SOA_COLUMN(M, m, float); //! Invariant mass of candidate (GeV/c2) +DECLARE_SOA_COLUMN(Pt, pt, float); //! Transverse momentum of candidate (GeV/c) +DECLARE_SOA_COLUMN(PtGen, ptGen, float); //! Transverse momentum of candidate (GeV/c) +DECLARE_SOA_COLUMN(P, p, float); //! Momentum of candidate (GeV/c) +DECLARE_SOA_COLUMN(Y, y, float); //! Rapidity of candidate +DECLARE_SOA_COLUMN(Eta, eta, float); //! Pseudorapidity of candidate +DECLARE_SOA_COLUMN(Phi, phi, float); //! Azimuth angle of candidate +DECLARE_SOA_COLUMN(E, e, float); //! Energy of candidate (GeV) +DECLARE_SOA_COLUMN(NSigTpcPi1, nSigTpcPi1, float); //! TPC Nsigma separation for prong1 with pion mass hypothesis +DECLARE_SOA_COLUMN(NSigTofPi1, nSigTofPi1, float); //! TOF Nsigma separation for prong1 with pion mass hypothesis +DECLARE_SOA_COLUMN(DecayLength, decayLength, float); //! Decay length of candidate (cm) +DECLARE_SOA_COLUMN(DecayLengthXY, decayLengthXY, float); //! Transverse decay length of candidate (cm) +DECLARE_SOA_COLUMN(DecayLengthNormalised, decayLengthNormalised, float); //! Normalised decay length of candidate +DECLARE_SOA_COLUMN(DecayLengthXYNormalised, decayLengthXYNormalised, float); //! Normalised transverse decay length of candidate +DECLARE_SOA_COLUMN(ImpactParameterProduct, impactParameterProduct, float); //! Impact parameter product of candidate +DECLARE_SOA_COLUMN(Cpa, cpa, float); //! Cosine pointing angle of candidate +DECLARE_SOA_COLUMN(CpaXY, cpaXY, float); //! Cosine pointing angle of candidate in transverse plane +DECLARE_SOA_COLUMN(MaxNormalisedDeltaIP, maxNormalisedDeltaIP, float); //! Maximum normalized difference between measured and expected impact parameter of candidate prongs +DECLARE_SOA_COLUMN(MlScoreSig, mlScoreSig, float); //! ML score for signal class +} // namespace hf_cand_bplus_lite + +DECLARE_SOA_TABLE(HfRedCandBpLites, "AOD", "HFREDCANDBPLITE", //! Table with some B+ properties + hf_cand::Chi2PCA, + hf_cand_bplus_lite::DecayLength, + hf_cand_bplus_lite::DecayLengthXY, + hf_cand_bplus_lite::DecayLengthNormalised, + hf_cand_bplus_lite::DecayLengthXYNormalised, + hf_cand_bplus_lite::MProng0, + hf_cand_bplus_lite::PtProng0, + hf_cand_bplus_lite::PtProng1, + hf_cand::ImpactParameter0, + hf_cand::ImpactParameter1, + hf_cand_bplus_lite::ImpactParameterProduct, + hf_cand_bplus_lite::NSigTpcPi1, + hf_cand_bplus_lite::NSigTofPi1, + hf_cand_bplus_reduced::Prong0MlScoreBkg, + hf_cand_bplus_reduced::Prong0MlScorePrompt, + hf_cand_bplus_reduced::Prong0MlScoreNonprompt, + hf_cand_bplus_lite::MlScoreSig, + hf_sel_candidate_bplus::IsSelBplusToD0Pi, + hf_cand_bplus_lite::M, + hf_cand_bplus_lite::Pt, + hf_cand_bplus_lite::Cpa, + hf_cand_bplus_lite::CpaXY, + hf_cand_bplus_lite::MaxNormalisedDeltaIP, + hf_cand_bplus_lite::Eta, + hf_cand_bplus_lite::Phi, + hf_cand_bplus_lite::Y, + hf_cand_2prong::FlagMcMatchRec, + hf_cand_2prong::OriginMcRec, + hf_cand_bplus_lite::PtGen); + +DECLARE_SOA_TABLE(HfRedBpMcCheck, "AOD", "HFREDBPMCCHECK", //! Table with MC decay type check + hf_cand_3prong::FlagMcMatchRec, + hf_cand_bplus_lite::MProng0, + hf_cand_bplus_lite::PtProng0, + hf_cand_bplus_lite::M, + hf_cand_bplus_lite::Pt, + hf_cand_bplus_lite::MlScoreSig, + hf_bplus_mc::PdgCodeBeautyMother, + hf_bplus_mc::PdgCodeProng0, + hf_bplus_mc::PdgCodeProng1, + hf_bplus_mc::PdgCodeProng2); +} // namespace o2::aod + // string definitions, used for histogram axis labels const TString stringPt = "#it{p}_{T} (GeV/#it{c})"; +const TString stringPtD = "#it{p}_{T}(D0) (GeV/#it{c});"; const TString bPlusCandTitle = "B+ candidates;"; const TString entries = "entries"; const TString bPlusCandMatch = "B+ candidates (matched);"; @@ -40,33 +114,50 @@ const TString mcParticleMatched = "MC particles (matched);"; /// B+ analysis task struct HfTaskBplusReduced { + Produces hfRedCandBpLite; + Produces hfRedBpMcCheck; + Configurable selectionFlagBplus{"selectionFlagBplus", 1, "Selection Flag for Bplus"}; Configurable yCandGenMax{"yCandGenMax", 0.5, "max. gen particle rapidity"}; Configurable yCandRecoMax{"yCandRecoMax", 0.8, "max. cand. rapidity"}; Configurable etaTrackMax{"etaTrackMax", 0.8, "max. track pseudo-rapidity"}; Configurable ptTrackMin{"ptTrackMin", 0.1, "min. track transverse momentum"}; + Configurable fillHistograms{"fillHistograms", true, "Flag to enable histogram filling"}; + Configurable fillSparses{"fillSparses", false, "Flag to enable sparse filling"}; + Configurable fillTree{"fillTree", false, "Flag to enable tree filling"}; + Configurable fillBackground{"fillBackground", false, "Flag to enable filling of background histograms/sparses/tree (only MC)"}; + Configurable downSampleBkgFactor{"downSampleBkgFactor", 1., "Fraction of background candidates to keep for ML trainings"}; + Configurable ptMaxForDownSample{"ptMaxForDownSample", 10., "Maximum pt for the application of the downsampling factor"}; Configurable> binsPt{"binsPt", std::vector{hf_cuts_bplus_to_d0_pi::vecBinsPt}, "pT bin limits"}; HfHelper hfHelper; Filter filterSelectCandidates = (aod::hf_sel_candidate_bplus::isSelBplusToD0Pi >= selectionFlagBplus); - HistogramRegistry registry{ - "registry", - {{"hPtProng0", bPlusCandTitle + "prong 0 #it{p}_{T} (GeV/#it{c});" + entries, {HistType::kTH1F, {{1000, 0., 50.}}}}, - {"hPtProng1", bPlusCandTitle + "prong 1 #it{p}_{T} (GeV/#it{c});" + entries, {HistType::kTH1F, {{200, 0., 10.}}}}, - {"hPtCand", bPlusCandTitle + "candidate #it{p}_{T} (GeV/#it{c});" + entries, {HistType::kTH1F, {{1000, 0., 50.}}}}, - {"hCentrality", "centrality;centrality percentile;" + entries, {HistType::kTH1F, {{100, 0., 100.}}}}, - {"hPtRecSig", bPlusCandMatch + "candidate #it{p}_{T} (GeV/#it{c});" + entries, {HistType::kTH1F, {{300, 0., 30.}}}}, - {"hPtRecBg", bPlusCandUnmatch + "candidate #it{p}_{T} (GeV/#it{c});" + entries, {HistType::kTH1F, {{300, 0., 30.}}}}, - {"hPtGenSig", bPlusCandMatch + "candidate #it{p}_{T}^{gen.} (GeV/#it{c});" + entries, {HistType::kTH1F, {{300, 0., 30.}}}}, - {"hPtGen", mcParticleMatched + "candidate #it{p}_{T} (GeV/#it{c});" + entries, {HistType::kTH1F, {{300, 0., 30.}}}}}}; + HistogramRegistry registry{"registry"}; + + using TracksPion = soa::Join; void init(InitContext&) { - const AxisSpec axisMass{150, 4.5, 6.0}; - const AxisSpec axisCPA{120, -1.1, 1.1}; - const AxisSpec axisCPAFiner{300, 0.85, 1.0}; + std::array processFuncData{doprocessData, doprocessDataWithDmesMl}; + if ((std::accumulate(processFuncData.begin(), processFuncData.end(), 0)) > 1) { + LOGP(fatal, "Only one process function for data can be enabled at a time."); + } + std::array processFuncMc{doprocessMc, doprocessMcWithDmesMl}; + if ((std::accumulate(processFuncMc.begin(), processFuncMc.end(), 0)) > 1) { + LOGP(fatal, "Only one process function for MC can be enabled at a time."); + } + + if (((doprocessData || doprocessDataWithDmesMl) && fillTree && downSampleBkgFactor >= 1.) || + ((doprocessMc || doprocessMcWithDmesMl) && fillTree && fillBackground && downSampleBkgFactor >= 1.)) { + LOGP(fatal, "Set downSampleBkgFactor below unity when filling tree with background."); + } + + const AxisSpec axisMlScore{100, 0.f, 1.f}; + const AxisSpec axisMassBplus{150, 4.5, 6.0}; + const AxisSpec axisMassD0{300, 1.75f, 2.05f}; + const AxisSpec axisCpa{120, -1.1, 1.1}; const AxisSpec axisPtProng{100, 0., 10.}; const AxisSpec axisD0Prong{200, -0.05, 0.05}; const AxisSpec axisImpParProd{200, -0.001, 0.001}; @@ -74,78 +165,160 @@ struct HfTaskBplusReduced { const AxisSpec axisNormDecLength{40, 0., 20}; const AxisSpec axisEta{100, -2., 2.}; const AxisSpec axisRapidity{100, -2., 2.}; + const AxisSpec axisPtD0{100, 0., 50.}; const AxisSpec axisPtB{(std::vector)binsPt, "#it{p}_{T}^{B^{+}} (GeV/#it{c})"}; - // histograms process - registry.add("hMass", bPlusCandTitle + "inv. mass #bar{D^{0}}#pi^{+} (GeV/#it{c}^{2});" + stringPt, {HistType::kTH2F, {axisMass, axisPtB}}); - registry.add("hDecLength", bPlusCandTitle + "decay length (cm);" + stringPt, {HistType::kTH2F, {axisDecLength, axisPtB}}); - registry.add("hDecLengthXY", bPlusCandTitle + "decay length xy (cm);" + stringPt, {HistType::kTH2F, {axisDecLength, axisPtB}}); - registry.add("hd0Prong0", bPlusCandTitle + "prong 0 DCAxy to prim. vertex (cm);" + stringPt, {HistType::kTH2F, {axisD0Prong, axisPtB}}); - registry.add("hd0Prong1", bPlusCandTitle + "prong 1 DCAxy to prim. vertex (cm);" + stringPt, {HistType::kTH2F, {axisD0Prong, axisPtB}}); - registry.add("hCPA", bPlusCandTitle + "candidate cosine of pointing angle;" + stringPt, {HistType::kTH2F, {axisCPA, axisPtB}}); - registry.add("hCPAxy", bPlusCandTitle + "candidate cosine of pointing angle xy;" + stringPt, {HistType::kTH2F, {axisCPA, axisPtB}}); - registry.add("hEta", bPlusCandTitle + "candidate #it{#eta};" + stringPt, {HistType::kTH2F, {axisEta, axisPtB}}); - registry.add("hRapidity", bPlusCandTitle + "candidate #it{y};" + stringPt, {HistType::kTH2F, {axisRapidity, axisPtB}}); - registry.add("hImpParErr", bPlusCandTitle + "candidate impact parameter error (cm);" + stringPt, {HistType::kTH2F, {{100, -1., 1.}, axisPtB}}); - registry.add("hDecLenErr", bPlusCandTitle + "candidate decay length error (cm);" + stringPt, {HistType::kTH2F, {{100, 0., 1.}, axisPtB}}); - registry.add("hDecLenXYErr", bPlusCandTitle + "candidate decay length xy error (cm);" + stringPt, {HistType::kTH2F, {{100, 0., 1.}, axisPtB}}); - registry.add("hd0d0", bPlusCandTitle + "candidate product of DCAxy to prim. vertex (cm^{2});" + stringPt, {HistType::kTH2F, {axisImpParProd, axisPtB}}); - registry.add("hInvMassD0", bPlusCandTitle + "prong0, D0 inv. mass (GeV/#it{c}^{2});" + stringPt, {HistType::kTH2F, {{500, 1.4, 2.4}, axisPtB}}); - registry.add("hCPAFinerBinning", bPlusCandTitle + "candidate cosine of pointing angle;" + stringPt, {HistType::kTH2F, {axisCPAFiner, axisPtB}}); - registry.add("hCPAxyFinerBinning", bPlusCandTitle + "candidate cosine of pointing angle xy;" + stringPt, {HistType::kTH2F, {axisCPAFiner, axisPtB}}); - // histograms processMC - Gen Level - registry.add("hEtaGen", mcParticleMatched + "candidate #it{#eta}^{gen};" + stringPt, {HistType::kTH2F, {axisEta, axisPtB}}); - registry.add("hYGen", mcParticleMatched + "candidate #it{y}^{gen};" + stringPt, {HistType::kTH2F, {axisRapidity, axisPtB}}); - registry.add("hPtProng0Gen", mcParticleMatched + "prong 0 #it{p}_{T}^{gen} (GeV/#it{c});" + stringPt, {HistType::kTH2F, {axisPtProng, axisPtB}}); - registry.add("hPtProng1Gen", mcParticleMatched + "prong 1 #it{p}_{T}^{gen} (GeV/#it{c});" + stringPt, {HistType::kTH2F, {axisPtProng, axisPtB}}); - registry.add("hYProng0Gen", mcParticleMatched + "prong 0 #it{y}^{gen};" + stringPt, {HistType::kTH2F, {axisRapidity, axisPtB}}); - registry.add("hYProng1Gen", mcParticleMatched + "prong 1 #it{y}^{gen};" + stringPt, {HistType::kTH2F, {axisRapidity, axisPtB}}); - registry.add("hEtaProng0Gen", mcParticleMatched + "prong 0 #it{#eta}^{gen};" + stringPt, {HistType::kTH2F, {axisEta, axisPtB}}); - registry.add("hEtaProng1Gen", mcParticleMatched + "prong 1 #it{#eta}^{gen};" + stringPt, {HistType::kTH2F, {axisEta, axisPtB}}); - registry.add("hPtProngsVsPtBGen", mcParticleMatched + "prong 0 #it{p}_{T}^{gen} (GeV/#it{c});prong 1 #it{p}_{T}^{gen} (GeV/#it{c});" + stringPt, {HistType::kTH3F, {axisPtProng, axisPtProng, axisPtB}}); - registry.add("hYProngsVsBplusGen", mcParticleMatched + "prong 0 #it{y}^{gen};prong 1 #it{y}^{gen};" + stringPt, {HistType::kTH3F, {axisRapidity, axisRapidity, axisPtB}}); - registry.add("hEtaProngsVsBplusGen", mcParticleMatched + "prong 0 #it{#eta}^{gen};prong 1 #it{#eta}^{gen};" + stringPt, {HistType::kTH3F, {axisEta, axisEta, axisPtB}}); - registry.add("hPtGenWithRapidityBelowHalf", "MC particles (generated - |#it{y}^{gen}|<0.5);candidate #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{300, 0., 30.}}}); - registry.add("hPtGenWithProngsInAcceptance", "MC particles (generated-daughters in acceptance);candidate #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{300, 0., 30.}}}); - registry.add("hEtaGenWithProngsInAcceptance", "MC particles (generated-daughters in acceptance);B^{0} candidate #it{#eta}^{gen};entries", {HistType::kTH2F, {{100, -2., 2.}, {(std::vector)binsPt, "#it{p}_{T} (GeV/#it{c})"}}}); - registry.add("hYGenWithProngsInAcceptance", "MC particles (generated-daughters in acceptance);B^{0} candidate #it{y}^{gen};entries", {HistType::kTH2F, {{100, -2., 2.}, {(std::vector)binsPt, "#it{p}_{T} (GeV/#it{c})"}}}); - // histograms processMC - Reco Level - registry.add("hCPARecSig", bPlusCandMatch + "candidate cosine of pointing angle;" + stringPt, {HistType::kTH2F, {axisCPA, axisPtB}}); - registry.add("hCPARecBg", bPlusCandUnmatch + "candidate cosine of pointing angle;" + stringPt, {HistType::kTH2F, {axisCPA, axisPtB}}); - registry.add("hCPAxyRecSig", bPlusCandMatch + "candidate CPAxy;" + stringPt, {HistType::kTH2F, {axisCPA, axisPtB}}); - registry.add("hCPAxyRecBg", bPlusCandUnmatch + "candidate CPAxy;" + stringPt, {HistType::kTH2F, {axisCPA, axisPtB}}); - registry.add("hEtaRecSig", bPlusCandMatch + "candidate #it{#eta};" + stringPt, {HistType::kTH2F, {axisEta, axisPtB}}); - registry.add("hEtaRecBg", bPlusCandUnmatch + "candidate #it{#eta};" + stringPt, {HistType::kTH2F, {axisEta, axisPtB}}); - registry.add("hRapidityRecSig", bPlusCandMatch + "candidate #it{y};" + stringPt, {HistType::kTH2F, {axisRapidity, axisPtB}}); - registry.add("hRapidityRecBg", bPlusCandUnmatch + "candidate #it{#y};" + stringPt, {HistType::kTH2F, {axisRapidity, axisPtB}}); - registry.add("hPtProng0RecSig", bPlusCandMatch + "prong 0 #it{p}_{T} (GeV/#it{c});" + stringPt, {HistType::kTH2F, {axisPtProng, axisPtB}}); - registry.add("hPtProng1RecSig", bPlusCandMatch + "prong 1 #it{p}_{T} (GeV/#it{c});" + stringPt, {HistType::kTH2F, {axisPtProng, axisPtB}}); - registry.add("hPtProngsVsBplusRecBg", bPlusCandUnmatch + "prong 0 #it{p}_{T} (GeV/#it{c});prong 1 #it{p}_{T} (GeV/#it{c});" + stringPt, {HistType::kTH3F, {axisPtProng, axisPtProng, axisPtB}}); - registry.add("hPtProng0RecBg", bPlusCandUnmatch + "prong 0 #it{p}_{T} (GeV/#it{c});" + stringPt, {HistType::kTH2F, {axisPtProng, axisPtB}}); - registry.add("hPtProng1RecBg", bPlusCandUnmatch + "prong 1 #it{p}_{T} (GeV/#it{c});" + stringPt, {HistType::kTH2F, {axisPtProng, axisPtB}}); - registry.add("hMassRecSig", bPlusCandMatch + "inv. mass #bar{D^{0}}#pi^{+} (GeV/#it{c}^{2});" + stringPt, {HistType::kTH2F, {axisMass, axisPtB}}); - registry.add("hMassRecBg", bPlusCandUnmatch + "inv. mass #bar{D^{0}}#pi^{+} (GeV/#it{c}^{2});" + stringPt, {HistType::kTH2F, {axisMass, axisPtB}}); - registry.add("hd0Prong0RecSig", bPlusCandMatch + "prong 0 DCAxy to prim. vertex (cm);" + stringPt, {HistType::kTH2F, {axisD0Prong, axisPtB}}); - registry.add("hd0Prong1RecSig", bPlusCandMatch + "prong 1 DCAxy to prim. vertex (cm);" + stringPt, {HistType::kTH2F, {axisD0Prong, axisPtB}}); - registry.add("hd0Prong0RecBg", bPlusCandUnmatch + "prong 0 DCAxy to prim. vertex (cm);" + stringPt, {HistType::kTH2F, {axisD0Prong, axisPtB}}); - registry.add("hd0Prong1RecBg", bPlusCandUnmatch + "prong 1 DCAxy to prim. vertex (cm);" + stringPt, {HistType::kTH2F, {axisD0Prong, axisPtB}}); - registry.add("hDecLengthRecSig", bPlusCandMatch + "candidate decay length (cm);" + stringPt, {HistType::kTH2F, {axisDecLength, axisPtB}}); - registry.add("hDecLengthXYRecSig", bPlusCandMatch + "candidate decay length xy (cm);" + stringPt, {HistType::kTH2F, {axisDecLength, axisPtB}}); - registry.add("hDecLengthRecBg", bPlusCandUnmatch + "candidate decay length (cm);" + stringPt, {HistType::kTH2F, {axisDecLength, axisPtB}}); - registry.add("hDecLengthXYRecBg", bPlusCandUnmatch + "candidate decay length xy(cm);" + stringPt, {HistType::kTH2F, {axisDecLength, axisPtB}}); - registry.add("hDecLengthNormRecSig", bPlusCandMatch + "candidate normalized decay length (cm);" + stringPt, {HistType::kTH2F, {axisNormDecLength, axisPtB}}); - registry.add("hDecLengthNormRecBg", bPlusCandUnmatch + "candidate normalized decay length (cm);" + stringPt, {HistType::kTH2F, {axisNormDecLength, axisPtB}}); - registry.add("hd0d0RecSig", bPlusCandMatch + "product of DCAxy to prim. vertex (cm^{2});" + stringPt, {HistType::kTH2F, {axisImpParProd, axisPtB}}); - registry.add("hd0d0RecBg", bPlusCandUnmatch + "product of DCAxy to prim. vertex (cm^{2});" + stringPt, {HistType::kTH2F, {axisImpParProd, axisPtB}}); - // MC histograms with finer binning - registry.add("hCPAFinerBinningRecSig", bPlusCandMatch + "candidate cosine of pointing angle;" + stringPt, {HistType::kTH2F, {axisCPAFiner, axisPtB}}); - registry.add("hCPAFinerBinningRecBg", bPlusCandMatch + "candidate cosine of pointing angle;" + stringPt, {HistType::kTH2F, {axisCPAFiner, axisPtB}}); - registry.add("hCPAxyFinerBinningRecSig", bPlusCandMatch + "candidate cosine of pointing angle;" + stringPt, {HistType::kTH2F, {axisCPAFiner, axisPtB}}); - registry.add("hCPAxyFinerBinningRecBg", bPlusCandMatch + "candidate cosine of pointing angle;" + stringPt, {HistType::kTH2F, {axisCPAFiner, axisPtB}}); - // histograms prong0(D0) - Reco Level - registry.add("hCPAD0RecSig", bPlusCandMatch + "prong0 (D^{0}) cosine of pointing angle;#it{p}_{T}(D0) (GeV/#it{c})" + entries, {HistType::kTH2F, {{220, 0., 1.1}, {120, 0., 60.}}}); - registry.add("hCPAD0RecBg", bPlusCandUnmatch + "prong0 (D^{0}) cosine of pointing angle;#it{p}_{T}(D0) (GeV/#it{c})" + entries, {HistType::kTH2F, {{220, 0., 1.1}, {120, 0., 60.}}}); - registry.add("hDecLengthD0RecSig", bPlusCandMatch + "prong0 D^{0} decay length (cm);#it{p}_{T}(D0) (GeV/#it{c})" + entries, {HistType::kTH2F, {{100, 0., 0.5}, {120, 0., 60.}}}); - registry.add("hDecLengthD0RecBg", bPlusCandUnmatch + "prong0 D^{0} candidate decay length (cm);#it{p}_{T}(D0) (GeV/#it{c})" + entries, {HistType::kTH2F, {{100, 0., 0.5}, {120, 0., 60.}}}); + const AxisSpec axisPtPi{100, 0.f, 10.f}; + + if (doprocessData || doprocessDataWithDmesMl || doprocessDataWithBplusMl) { + if (fillHistograms) { + registry.add("hMass", bPlusCandTitle + "inv. mass #bar{D^{0}}#pi^{+} (GeV/#it{c}^{2});" + stringPt, {HistType::kTH2F, {axisMassBplus, axisPtB}}); + registry.add("hPtCand", bPlusCandTitle + "candidate #it{p}_{T} (GeV/#it{c});" + entries, {HistType::kTH1F, {{1000, 0., 50.}}}); + registry.add("hPtProng0", bPlusCandTitle + "prong 0 #it{p}_{T}^{gen} (GeV/#it{c});" + stringPt, {HistType::kTH2F, {axisPtProng, axisPtB}}); + registry.add("hPtProng1", bPlusCandTitle + "prong 1 #it{p}_{T}^{gen} (GeV/#it{c});" + stringPt, {HistType::kTH2F, {axisPtProng, axisPtB}}); + registry.add("hDecLength", bPlusCandTitle + "decay length (cm);" + stringPt, {HistType::kTH2F, {axisDecLength, axisPtB}}); + registry.add("hDecLengthXy", bPlusCandTitle + "decay length xy (cm);" + stringPt, {HistType::kTH2F, {axisDecLength, axisPtB}}); + registry.add("hNormDecLengthXy", bPlusCandTitle + "Norm. decay length xy (cm);" + stringPt, {HistType::kTH2F, {axisNormDecLength, axisPtB}}); + registry.add("hd0Prong0", bPlusCandTitle + "prong 0 DCAxy to prim. vertex (cm);" + stringPt, {HistType::kTH2F, {axisD0Prong, axisPtB}}); + registry.add("hd0Prong1", bPlusCandTitle + "prong 1 DCAxy to prim. vertex (cm);" + stringPt, {HistType::kTH2F, {axisD0Prong, axisPtB}}); + registry.add("hCpa", bPlusCandTitle + "candidate cosine of pointing angle;" + stringPt, {HistType::kTH2F, {axisCpa, axisPtB}}); + registry.add("hCpaXy", bPlusCandTitle + "candidate cosine of pointing angle xy;" + stringPt, {HistType::kTH2F, {axisCpa, axisPtB}}); + registry.add("hEta", bPlusCandTitle + "candidate #it{#eta};" + stringPt, {HistType::kTH2F, {axisEta, axisPtB}}); + registry.add("hRapidity", bPlusCandTitle + "candidate #it{y};" + stringPt, {HistType::kTH2F, {axisRapidity, axisPtB}}); + registry.add("hd0d0", bPlusCandTitle + "candidate product of DCAxy to prim. vertex (cm^{2});" + stringPt, {HistType::kTH2F, {axisImpParProd, axisPtB}}); + registry.add("hInvMassD0", bPlusCandTitle + "prong0, D0 inv. mass (GeV/#it{c}^{2});" + stringPt, {HistType::kTH2F, {axisMassD0, axisPtD0}}); + registry.add("hDecLengthD", bPlusCandTitle + "#it{p}_{T}(D^{0}) (GeV/#it{c});D^{0} candidate decay length (cm);entries", {HistType::kTH2F, {axisPtD0, axisDecLength}}); + registry.add("hDecLengthXyD", bPlusCandTitle + "#it{p}_{T}(D^{0}) (GeV/#it{c});decay length XY (cm);entries", {HistType::kTH2F, {axisPtD0, axisDecLength}}); + registry.add("hCpaD", bPlusCandTitle + "#it{p}_{T}(D^{0}) (GeV/#it{c});D^{0} candidate cos(#vartheta_{P});entries", {HistType::kTH2F, {axisPtD0, axisCpa}}); + registry.add("hCpaXyD", bPlusCandTitle + "#it{p}_{T}(D^{0}) (GeV/#it{c});D^{0} candidate cos(#vartheta_{P}^{XY});entries", {HistType::kTH2F, {axisPtD0, axisCpa}}); + + // ML scores of D0 daughter + if (doprocessDataWithDmesMl) { + registry.add("hMlScoreBkgD", bPlusCandTitle + "#it{p}_{T}(D0) (GeV/#it{c});prong0, D0 ML background score;entries", {HistType::kTH2F, {axisPtD0, axisMlScore}}); + registry.add("hMlScorePromptD", bPlusCandTitle + "#it{p}_{T}(D0) (GeV/#it{c});prong0, D0 ML prompt score;entries", {HistType::kTH2F, {axisPtD0, axisMlScore}}); + registry.add("hMlScoreNonPromptD", bPlusCandTitle + "#it{p}_{T}(D0) (GeV/#it{c});prong0, D0 ML nonprompt score;entries", {HistType::kTH2F, {axisPtD0, axisMlScore}}); + } + + // ML scores of B+ candidate + if (doprocessDataWithBplusMl) { + registry.add("hMlScoreSigBplus", bPlusCandTitle + "#it{p}_{T} (GeV/#it{c});prong0, B^{+} ML signal score;entries", {HistType::kTH2F, {axisPtB, axisMlScore}}); + } + } + if (fillSparses) { + if (!(doprocessDataWithDmesMl || doprocessDataWithBplusMl)) { + registry.add("hMassPtCutVars", bPlusCandTitle + "#it{M} (D0#pi) (GeV/#it{c}^{2});#it{p}_{T}(B^{+}) (GeV/#it{c});B^{+} candidate decay length (cm);B^{+} candidate norm. decay length XY (cm);B^{+} candidate impact parameter product (cm);B^{+} candidate cos(#vartheta_{P});#it{M} (K#pi) (GeV/#it{c}^{2});#it{p}_{T}(D0) (GeV/#it{c});D0 candidate decay length (cm);D0 candidate cos(#vartheta_{P})", {HistType::kTHnSparseF, {axisMassBplus, axisPtB, axisDecLength, axisNormDecLength, axisImpParProd, axisCpa, axisMassD0, axisPtD0, axisDecLength, axisCpa}}); + } else { + registry.add("hMassPtCutVars", bPlusCandTitle + "#it{M} (D0#pi) (GeV/#it{c}^{2});#it{p}_{T}(B^{+}) (GeV/#it{c});B^{+} candidate decay length (cm);B^{+} candidate norm. decay length XY (cm);B^{+} candidate impact parameter product (cm);B^{+} candidate cos(#vartheta_{P});#it{M} (K#pi) (GeV/#it{c}^{2});#it{p}_{T}(D0) (GeV/#it{c});D0 candidate ML score bkg;D0 candidate ML score nonprompt", {HistType::kTHnSparseF, {axisMassBplus, axisPtB, axisDecLength, axisNormDecLength, axisImpParProd, axisCpa, axisMassD0, axisPtD0, axisMlScore, axisMlScore}}); + } + } + } + + // histograms processMC + if (doprocessMc || doprocessMcWithDecayTypeCheck || doprocessMcWithDmesMl || doprocessMcWithBplusMl || doprocessMcWithBplusMlAndDecayTypeCheck) { + if (fillHistograms) { + // Gen Level + registry.add("hEtaGen", mcParticleMatched + "candidate #it{#eta}^{gen};" + stringPt, {HistType::kTH2F, {axisEta, axisPtB}}); + registry.add("hYGen", mcParticleMatched + "candidate #it{y}^{gen};" + stringPt, {HistType::kTH2F, {axisRapidity, axisPtB}}); + registry.add("hPtProng0Gen", mcParticleMatched + "prong 0 #it{p}_{T}^{gen} (GeV/#it{c});" + stringPt, {HistType::kTH2F, {axisPtProng, axisPtB}}); + registry.add("hPtProng1Gen", mcParticleMatched + "prong 1 #it{p}_{T}^{gen} (GeV/#it{c});" + stringPt, {HistType::kTH2F, {axisPtProng, axisPtB}}); + registry.add("hYProng0Gen", mcParticleMatched + "prong 0 #it{y}^{gen};" + stringPt, {HistType::kTH2F, {axisRapidity, axisPtB}}); + registry.add("hYProng1Gen", mcParticleMatched + "prong 1 #it{y}^{gen};" + stringPt, {HistType::kTH2F, {axisRapidity, axisPtB}}); + registry.add("hEtaProng0Gen", mcParticleMatched + "prong 0 #it{#eta}^{gen};" + stringPt, {HistType::kTH2F, {axisEta, axisPtB}}); + registry.add("hEtaProng1Gen", mcParticleMatched + "prong 1 #it{#eta}^{gen};" + stringPt, {HistType::kTH2F, {axisEta, axisPtB}}); + registry.add("hPtProngsVsPtBGen", mcParticleMatched + "prong 0 #it{p}_{T}^{gen} (GeV/#it{c});prong 1 #it{p}_{T}^{gen} (GeV/#it{c});" + stringPt, {HistType::kTH3F, {axisPtProng, axisPtProng, axisPtB}}); + registry.add("hYProngsVsBplusGen", mcParticleMatched + "prong 0 #it{y}^{gen};prong 1 #it{y}^{gen};" + stringPt, {HistType::kTH3F, {axisRapidity, axisRapidity, axisPtB}}); + registry.add("hEtaProngsVsBplusGen", mcParticleMatched + "prong 0 #it{#eta}^{gen};prong 1 #it{#eta}^{gen};" + stringPt, {HistType::kTH3F, {axisEta, axisEta, axisPtB}}); + registry.add("hYGenWithProngsInAcceptance", "MC particles (generated-daughters in acceptance);B^{+} candidate #it{y}^{gen};entries", {HistType::kTH2F, {{100, -2., 2.}, {(std::vector)binsPt, "#it{p}_{T} (GeV/#it{c})"}}}); + // Reco Level - Signal + registry.add("hPtRecSig", bPlusCandMatch + "candidate #it{p}_{T} (GeV/#it{c});" + entries, {HistType::kTH1F, {{300, 0., 30.}}}); + registry.add("hCpaRecSig", bPlusCandMatch + "candidate cosine of pointing angle;" + stringPt, {HistType::kTH2F, {axisCpa, axisPtB}}); + registry.add("hCpaXyRecSig", bPlusCandMatch + "candidate CpaXy;" + stringPt, {HistType::kTH2F, {axisCpa, axisPtB}}); + registry.add("hEtaRecSig", bPlusCandMatch + "candidate #it{#eta};" + stringPt, {HistType::kTH2F, {axisEta, axisPtB}}); + registry.add("hRapidityRecSig", bPlusCandMatch + "candidate #it{y};" + stringPt, {HistType::kTH2F, {axisRapidity, axisPtB}}); + registry.add("hPtProng0RecSig", bPlusCandMatch + "prong 0 #it{p}_{T} (GeV/#it{c});" + stringPt, {HistType::kTH2F, {axisPtProng, axisPtB}}); + registry.add("hPtProng1RecSig", bPlusCandMatch + "prong 1 #it{p}_{T} (GeV/#it{c});" + stringPt, {HistType::kTH2F, {axisPtProng, axisPtB}}); + registry.add("hMassRecSig", bPlusCandMatch + "inv. mass #bar{D^{0}}#pi^{+} (GeV/#it{c}^{2});" + stringPt, {HistType::kTH2F, {axisMassBplus, axisPtB}}); + registry.add("hd0Prong0RecSig", bPlusCandMatch + "prong 0 DCAxy to prim. vertex (cm);" + stringPt, {HistType::kTH2F, {axisD0Prong, axisPtB}}); + registry.add("hd0Prong1RecSig", bPlusCandMatch + "prong 1 DCAxy to prim. vertex (cm);" + stringPt, {HistType::kTH2F, {axisD0Prong, axisPtB}}); + registry.add("hDecLengthRecSig", bPlusCandMatch + "candidate decay length (cm);" + stringPt, {HistType::kTH2F, {axisDecLength, axisPtB}}); + registry.add("hDecLengthXyRecSig", bPlusCandMatch + "candidate decay length xy (cm);" + stringPt, {HistType::kTH2F, {axisDecLength, axisPtB}}); + registry.add("hNormDecLengthXyRecSig", bPlusCandMatch + "candidate normalized decay length (cm);" + stringPt, {HistType::kTH2F, {axisNormDecLength, axisPtB}}); + registry.add("hd0d0RecSig", bPlusCandMatch + "product of DCAxy to prim. vertex (cm^{2});" + stringPt, {HistType::kTH2F, {axisImpParProd, axisPtB}}); + registry.add("hCpaD0RecSig", bPlusCandMatch + "prong0 (D^{0}) cosine of pointing angle;" + stringPtD + entries, {HistType::kTH2F, {{220, 0., 1.1}, {120, 0., 60.}}}); + registry.add("hCpaXyD0RecSig", bPlusCandMatch + "prong0 (D^{0}) cosine of pointing angle;" + stringPtD + entries, {HistType::kTH2F, {{220, 0., 1.1}, {120, 0., 60.}}}); + registry.add("hDecLengthD0RecSig", bPlusCandMatch + "prong0 D^{0} decay length (cm);" + stringPtD + entries, {HistType::kTH2F, {{100, 0., 0.5}, {120, 0., 60.}}}); + registry.add("hDecLengthXyD0RecSig", bPlusCandMatch + "prong0 D^{0} decay length XY (cm);" + stringPtD + entries, {HistType::kTH2F, {{100, 0., 0.5}, {120, 0., 60.}}}); + // background + if (fillBackground) { + registry.add("hPtRecBg", bPlusCandUnmatch + "candidate #it{p}_{T} (GeV/#it{c});" + entries, {HistType::kTH1F, {{300, 0., 30.}}}); + registry.add("hCpaRecBg", bPlusCandUnmatch + "candidate cosine of pointing angle;" + stringPt, {HistType::kTH2F, {axisCpa, axisPtB}}); + registry.add("hCpaXyRecBg", bPlusCandUnmatch + "candidate CpaXy;" + stringPt, {HistType::kTH2F, {axisCpa, axisPtB}}); + registry.add("hEtaRecBg", bPlusCandUnmatch + "candidate #it{#eta};" + stringPt, {HistType::kTH2F, {axisEta, axisPtB}}); + registry.add("hRapidityRecBg", bPlusCandUnmatch + "candidate #it{#y};" + stringPt, {HistType::kTH2F, {axisRapidity, axisPtB}}); + registry.add("hPtProngsVsBplusRecBg", bPlusCandUnmatch + "prong 0 #it{p}_{T} (GeV/#it{c});prong 1 #it{p}_{T} (GeV/#it{c});" + stringPt, {HistType::kTH3F, {axisPtProng, axisPtProng, axisPtB}}); + registry.add("hPtProng0RecBg", bPlusCandUnmatch + "prong 0 #it{p}_{T} (GeV/#it{c});" + stringPt, {HistType::kTH2F, {axisPtProng, axisPtB}}); + registry.add("hPtProng1RecBg", bPlusCandUnmatch + "prong 1 #it{p}_{T} (GeV/#it{c});" + stringPt, {HistType::kTH2F, {axisPtProng, axisPtB}}); + registry.add("hMassRecBg", bPlusCandUnmatch + "inv. mass #bar{D^{0}}#pi^{+} (GeV/#it{c}^{2});" + stringPt, {HistType::kTH2F, {axisMassBplus, axisPtB}}); + registry.add("hd0Prong0RecBg", bPlusCandUnmatch + "prong 0 DCAxy to prim. vertex (cm);" + stringPt, {HistType::kTH2F, {axisD0Prong, axisPtB}}); + registry.add("hd0Prong1RecBg", bPlusCandUnmatch + "prong 1 DCAxy to prim. vertex (cm);" + stringPt, {HistType::kTH2F, {axisD0Prong, axisPtB}}); + registry.add("hDecLengthRecBg", bPlusCandUnmatch + "candidate decay length (cm);" + stringPt, {HistType::kTH2F, {axisDecLength, axisPtB}}); + registry.add("hDecLengthXyRecBg", bPlusCandUnmatch + "candidate decay length xy(cm);" + stringPt, {HistType::kTH2F, {axisDecLength, axisPtB}}); + registry.add("hNormDecLengthXyRecBg", bPlusCandUnmatch + "candidate normalized decay length (cm);" + stringPt, {HistType::kTH2F, {axisNormDecLength, axisPtB}}); + registry.add("hd0d0RecBg", bPlusCandUnmatch + "product of DCAxy to prim. vertex (cm^{2});" + stringPt, {HistType::kTH2F, {axisImpParProd, axisPtB}}); + registry.add("hCpaD0RecBg", bPlusCandUnmatch + "prong0 (D^{0}) cosine of pointing angle;" + stringPtD + entries, {HistType::kTH2F, {{220, 0., 1.1}, {120, 0., 60.}}}); + registry.add("hCpaXyD0RecBg", bPlusCandUnmatch + "prong0 (D^{0}) cosine of pointing angle;" + stringPtD + entries, {HistType::kTH2F, {{220, 0., 1.1}, {120, 0., 60.}}}); + registry.add("hDecLengthD0RecBg", bPlusCandUnmatch + "prong0 D^{0} decay length (cm);" + stringPtD + entries, {HistType::kTH2F, {{100, 0., 0.5}, {120, 0., 60.}}}); + registry.add("hDecLengthXyD0RecBg", bPlusCandUnmatch + "prong0 D^{0} decay length XY (cm);" + stringPtD + entries, {HistType::kTH2F, {{100, 0., 0.5}, {120, 0., 60.}}}); + } + // MC checks + if (doprocessMcWithDecayTypeCheck || doprocessMcWithBplusMlAndDecayTypeCheck) { + constexpr uint8_t kNBinsDecayTypeMc = hf_cand_bplus::DecayTypeMc::NDecayTypeMc; + TString labels[kNBinsDecayTypeMc]; + labels[hf_cand_bplus::DecayTypeMc::BplusToD0PiToKPiPi] = "B^{+} #rightarrow (#overline{D^{0}} #rightarrow K^{#plus} #pi^{#minus}) #pi^{#plus}"; + labels[hf_cand_bplus::DecayTypeMc::PartlyRecoDecay] = "Partly reconstructed decay channel"; + labels[hf_cand_bplus::DecayTypeMc::OtherDecay] = "Other decays"; + static const AxisSpec axisDecayType = {kNBinsDecayTypeMc, 0.5, kNBinsDecayTypeMc + 0.5, ""}; + registry.add("hDecayTypeMc", "DecayType", {HistType::kTH3F, {axisDecayType, axisMassBplus, axisPtB}}); + for (uint8_t iBin = 0; iBin < kNBinsDecayTypeMc; ++iBin) { + registry.get(HIST("hDecayTypeMc"))->GetXaxis()->SetBinLabel(iBin + 1, labels[iBin]); + } + } + // ML scores of D0 daughter + if (doprocessMcWithDmesMl) { + // signal + registry.add("hMlScoreBkgDRecSig", bPlusCandMatch + stringPtD + "prong0, D0 ML background score;entries", {HistType::kTH2F, {axisPtD0, axisMlScore}}); + registry.add("hMlScorePromptDRecSig", bPlusCandMatch + stringPtD + "prong0, D0 ML prompt score;entries", {HistType::kTH2F, {axisPtD0, axisMlScore}}); + registry.add("hMlScoreNonPromptDRecSig", bPlusCandMatch + stringPtD + "prong0, D0 ML nonprompt score;entries", {HistType::kTH2F, {axisPtD0, axisMlScore}}); + // background + registry.add("hMlScoreBkgDRecBg", bPlusCandUnmatch + stringPtD + "prong0, D0 ML background score;entries", {HistType::kTH2F, {axisPtD0, axisMlScore}}); + registry.add("hMlScorePromptDRecBg", bPlusCandUnmatch + stringPtD + "prong0, D0 ML prompt score;entries", {HistType::kTH2F, {axisPtD0, axisMlScore}}); + registry.add("hMlScoreNonPromptDRecBg", bPlusCandUnmatch + stringPtD + "prong0, D0 ML nonprompt score;entries", {HistType::kTH2F, {axisPtD0, axisMlScore}}); + } + // ML scores of B+ candidate + if (doprocessMcWithBplusMl || doprocessMcWithBplusMlAndDecayTypeCheck) { + // signal + registry.add("hMlScoreSigBplusRecSig", bPlusCandMatch + "#it{p}_{T}(B^{+}) (GeV/#it{c});prong0, B^{+} ML signal score;entries", {HistType::kTH2F, {axisPtB, axisMlScore}}); + // background + registry.add("hMlScoreSigBplusRecBg", bPlusCandUnmatch + "#it{p}_{T}(B^{+}) (GeV/#it{c});prong0, B^{+} ML signal score;entries", {HistType::kTH2F, {axisPtB, axisMlScore}}); + } + } + if (fillSparses) { + // gen sparses + registry.add("hPtYGenSig", "B^{+} particles (generated);#it{p}_{T}(B^{+}) (GeV/#it{c});#it{y}(B^{+})", {HistType::kTHnSparseF, {axisPtB, axisEta}}); + registry.add("hPtYWithProngsInAccepanceGenSig", "B^{+} particles (generated-daughters in acceptance);#it{p}_{T}(B^{+}) (GeV/#it{c});#it{y}(B^{+})", {HistType::kTHnSparseF, {axisPtB, axisEta}}); + // reco sparses + if (!doprocessDataWithDmesMl) { + registry.add("hMassPtCutVarsRecSig", bPlusCandMatch + "#it{M} (D^{0}#pi) (GeV/#it{c}^{2});#it{p}_{T}(B^{+}) (GeV/#it{c});B^{+} candidate decay length (cm);B^{+} candidate norm. decay length XY (cm);B^{+} candidate impact parameter product (cm);B^{+} candidate cos(#vartheta_{P});#it{M} (K#pi) (GeV/#it{c}^{2});#it{p}_{T}(D0) (GeV/#it{c});D0 candidate decay length (cm);D0 candidate cos(#vartheta_{P})", {HistType::kTHnSparseF, {axisMassBplus, axisPtB, axisDecLength, axisNormDecLength, axisImpParProd, axisCpa, axisMassD0, axisPtD0, axisDecLength, axisCpa}}); + if (fillBackground) { + registry.add("hMassPtCutVarsRecBg", bPlusCandUnmatch + "#it{M} (D^{0}#pi) (GeV/#it{c}^{2});#it{p}_{T}(B^{+}) (GeV/#it{c});B^{+} candidate decay length (cm);B^{+} candidate norm. decay length XY (cm);B^{+} candidate impact parameter product (cm);B^{+} candidate cos(#vartheta_{P});#it{M} (K#pi) (GeV/#it{c}^{2});#it{p}_{T}(D0) (GeV/#it{c});D0 candidate decay length (cm);D0 candidate cos(#vartheta_{P})", {HistType::kTHnSparseF, {axisMassBplus, axisPtB, axisDecLength, axisNormDecLength, axisImpParProd, axisCpa, axisMassD0, axisPtD0, axisDecLength, axisCpa}}); + } + } else { + registry.add("hMassPtCutVarsRecSig", bPlusCandMatch + "#it{M} (D^{0}#pi) (GeV/#it{c}^{2});#it{p}_{T}(B^{+}) (GeV/#it{c});B^{+} candidate decay length (cm);B^{+} candidate norm. decay length XY (cm);B^{+} candidate impact parameter product (cm);B^{+} candidate cos(#vartheta_{P});#it{M} (K#pi) (GeV/#it{c}^{2});#it{p}_{T}(D0) (GeV/#it{c});D0 candidate ML score bkg;D0 candidate ML score nonprompt", {HistType::kTHnSparseF, {axisMassBplus, axisPtB, axisDecLength, axisNormDecLength, axisImpParProd, axisCpa, axisMassD0, axisPtD0, axisMlScore, axisMlScore}}); + if (fillBackground) { + registry.add("hMassPtCutVarsRecBg", bPlusCandUnmatch + "#it{M} (D^{0}#pi) (GeV/#it{c}^{2});#it{p}_{T}(B^{+}) (GeV/#it{c});B^{+} candidate decay length (cm);B^{+} candidate norm. decay length XY (cm);B^{+} candidate impact parameter product (cm);B^{+} candidate cos(#vartheta_{P});#it{M} (K#pi) (GeV/#it{c}^{2});#it{p}_{T}(D0) (GeV/#it{c});D0 candidate ML score bkg;D0 candidate ML score nonprompt", {HistType::kTHnSparseF, {axisMassBplus, axisPtB, axisDecLength, axisNormDecLength, axisImpParProd, axisCpa, axisMassD0, axisPtD0, axisMlScore, axisMlScore}}); + } + } + } + } } /// Selection of B+ daughter in geometrical acceptance @@ -158,152 +331,411 @@ struct HfTaskBplusReduced { return std::abs(etaProng) <= etaTrackMax && ptProng >= ptTrackMin; } - void process(soa::Filtered> const& candidates, - aod::HfRed2Prongs const&, - aod::HfRedTracks const&) + /// Fill candidate information at reconstruction level + /// \param doMc is the flag to enable the filling with MC information + /// \param withDecayTypeCheck is the flag to enable MC with decay type check + /// \param withDmesMl is the flag to enable the filling with ML scores for the D0 daughter + /// \param withBplusMl is the flag to enable the filling with ML scores for the B+ candidate + /// \param candidate is the B+ candidate + /// \param candidatesD is the table with D0 candidates + template + void fillCand(Cand const& candidate, + aod::HfRed2Prongs const& candidatesD) + { + auto ptCandBplus = candidate.pt(); + auto invMassBplus = hfHelper.invMassBplusToD0Pi(candidate); + auto candD0 = candidate.template prong0_as(); + auto candPi = candidate.template prong1_as(); + auto ptD0 = candidate.ptProng0(); + auto invMassD0 = (candPi.signed1Pt() < 0) ? candD0.invMassD0() : candD0.invMassD0Bar(); + std::array posPv{candidate.posX(), candidate.posY(), candidate.posZ()}; + std::array posSvD{candD0.xSecondaryVertex(), candD0.ySecondaryVertex(), candD0.zSecondaryVertex()}; + std::array momD{candD0.pVector()}; + auto cpaD0 = RecoDecay::cpa(posPv, posSvD, momD); + auto cpaXyD0 = RecoDecay::cpaXY(posPv, posSvD, momD); + auto decLenD0 = RecoDecay::distance(posPv, posSvD); + auto decLenXyD0 = RecoDecay::distanceXY(posPv, posSvD); + + int8_t flagMcMatchRec = 0; + bool isSignal = false; + if constexpr (doMc) { + flagMcMatchRec = candidate.flagMcMatchRec(); + isSignal = TESTBIT(std::abs(flagMcMatchRec), hf_cand_bplus::DecayTypeMc::BplusToD0PiToKPiPi); + } + + if (fillHistograms) { + if constexpr (doMc) { + if (isSignal) { + registry.fill(HIST("hMassRecSig"), invMassBplus, ptCandBplus); + registry.fill(HIST("hPtRecSig"), ptCandBplus); + registry.fill(HIST("hPtProng0RecSig"), ptD0, ptCandBplus); + registry.fill(HIST("hPtProng1RecSig"), candidate.ptProng1(), ptCandBplus); + registry.fill(HIST("hCpaRecSig"), candidate.cpa(), ptCandBplus); + registry.fill(HIST("hCpaXyRecSig"), candidate.cpaXY(), ptCandBplus); + registry.fill(HIST("hEtaRecSig"), candidate.eta(), ptCandBplus); + registry.fill(HIST("hRapidityRecSig"), hfHelper.yBplus(candidate), ptCandBplus); + registry.fill(HIST("hDecLengthRecSig"), candidate.decayLength(), ptCandBplus); + registry.fill(HIST("hDecLengthXyRecSig"), candidate.decayLengthXY(), ptCandBplus); + registry.fill(HIST("hNormDecLengthXyRecSig"), candidate.decayLengthXYNormalised(), ptCandBplus); + registry.fill(HIST("hd0Prong0RecSig"), candidate.impactParameter0(), ptCandBplus); + registry.fill(HIST("hd0Prong1RecSig"), candidate.impactParameter1(), ptCandBplus); + registry.fill(HIST("hd0d0RecSig"), candidate.impactParameterProduct(), ptCandBplus); + registry.fill(HIST("hPtProng0RecSig"), candidate.ptProng0(), ptCandBplus); + registry.fill(HIST("hPtProng1RecSig"), candidate.ptProng1(), ptCandBplus); + registry.fill(HIST("hDecLengthD0RecSig"), decLenD0, candidate.ptProng0()); + registry.fill(HIST("hDecLengthXyD0RecSig"), decLenXyD0, ptD0); + registry.fill(HIST("hCpaD0RecSig"), cpaD0, ptD0); + registry.fill(HIST("hCpaXyD0RecSig"), cpaXyD0, ptD0); + if constexpr (withDecayTypeCheck) { + registry.fill(HIST("hDecayTypeMc"), 1 + hf_cand_bplus::DecayTypeMc::BplusToD0PiToKPiPi, invMassBplus, ptCandBplus); + } + if constexpr (withDmesMl) { + registry.fill(HIST("hMlScoreBkgDRecSig"), ptD0, candidate.prong0MlScoreBkg()); + registry.fill(HIST("hMlScorePromptDRecSig"), ptD0, candidate.prong0MlScorePrompt()); + registry.fill(HIST("hMlScoreNonPromptDRecSig"), ptD0, candidate.prong0MlScoreNonprompt()); + } + if constexpr (withBplusMl) { + registry.fill(HIST("hMlScoreSigBplusRecSig"), ptCandBplus, candidate.mlProbBplusToD0Pi()); + } + } else if (fillBackground) { + registry.fill(HIST("hMassRecBg"), invMassBplus, ptCandBplus); + registry.fill(HIST("hPtRecBg"), ptCandBplus); + registry.fill(HIST("hPtProng0RecBg"), candidate.ptProng0(), ptCandBplus); + registry.fill(HIST("hPtProng1RecBg"), candidate.ptProng1(), ptCandBplus); + registry.fill(HIST("hCpaRecBg"), candidate.cpa(), ptCandBplus); + registry.fill(HIST("hCpaXyRecBg"), candidate.cpaXY(), ptCandBplus); + registry.fill(HIST("hEtaRecBg"), candidate.eta(), ptCandBplus); + registry.fill(HIST("hRapidityRecBg"), hfHelper.yBplus(candidate), ptCandBplus); + registry.fill(HIST("hDecLengthRecBg"), candidate.decayLength(), ptCandBplus); + registry.fill(HIST("hDecLengthXyRecBg"), candidate.decayLengthXY(), ptCandBplus); + registry.fill(HIST("hNormDecLengthXyRecBg"), candidate.decayLengthXYNormalised(), ptCandBplus); + registry.fill(HIST("hd0Prong0RecBg"), candidate.impactParameter0(), ptCandBplus); + registry.fill(HIST("hd0Prong1RecBg"), candidate.impactParameter1(), ptCandBplus); + registry.fill(HIST("hd0d0RecBg"), candidate.impactParameterProduct(), ptCandBplus); + registry.fill(HIST("hDecLengthD0RecBg"), decLenD0, candidate.ptProng0()); + registry.fill(HIST("hInvMassDRecBg"), invMassD0, ptD0); + registry.fill(HIST("hDecLengthDRecBg"), decLenD0, ptD0); + registry.fill(HIST("hDecLengthXyDRecBg"), decLenXyD0, ptD0); + registry.fill(HIST("hCpaDRecBg"), cpaD0, ptD0); + registry.fill(HIST("hCpaXyDRecBg"), cpaXyD0, ptD0); + if constexpr (withDmesMl) { + registry.fill(HIST("hMlScoreBkgDRecBg"), ptD0, candidate.prong0MlScoreBkg()); + registry.fill(HIST("hMlScorePromptDRecBg"), ptD0, candidate.prong0MlScorePrompt()); + registry.fill(HIST("hMlScoreNonPromptDRecBg"), ptD0, candidate.prong0MlScoreNonprompt()); + } + if constexpr (withBplusMl) { + registry.fill(HIST("hMlScoreSigBplusRecBg"), ptCandBplus, candidate.mlProbBplusToD0Pi()); + } + } else if constexpr (withDecayTypeCheck) { + if (TESTBIT(flagMcMatchRec, hf_cand_bplus::DecayTypeMc::PartlyRecoDecay)) { // Partly reconstructed decay channel + registry.fill(HIST("hDecayTypeMc"), 1 + hf_cand_bplus::DecayTypeMc::PartlyRecoDecay, invMassBplus, ptCandBplus); + } else { + registry.fill(HIST("hDecayTypeMc"), 1 + hf_cand_bplus::DecayTypeMc::OtherDecay, invMassBplus, ptCandBplus); + } + } + } else { + registry.fill(HIST("hMass"), invMassBplus, ptCandBplus); + registry.fill(HIST("hPtCand"), ptCandBplus); + registry.fill(HIST("hPtProng0"), ptD0, ptCandBplus); + registry.fill(HIST("hPtProng1"), candidate.ptProng1(), ptCandBplus); + registry.fill(HIST("hd0d0"), candidate.impactParameterProduct(), ptCandBplus); + registry.fill(HIST("hDecLength"), candidate.decayLength(), ptCandBplus); + registry.fill(HIST("hDecLengthXy"), candidate.decayLengthXY(), ptCandBplus); + registry.fill(HIST("hNormDecLengthXy"), candidate.decayLengthXY() / candidate.errorDecayLengthXY(), ptCandBplus); + registry.fill(HIST("hd0Prong0"), candidate.impactParameter0(), ptCandBplus); + registry.fill(HIST("hd0Prong1"), candidate.impactParameter1(), ptCandBplus); + registry.fill(HIST("hCpa"), candidate.cpa(), ptCandBplus); + registry.fill(HIST("hCpaXy"), candidate.cpaXY(), ptCandBplus); + registry.fill(HIST("hEta"), candidate.eta(), ptCandBplus); + registry.fill(HIST("hRapidity"), hfHelper.yBplus(candidate), ptCandBplus); + registry.fill(HIST("hInvMassD0"), invMassD0, ptCandBplus); + registry.fill(HIST("hDecLengthD0"), decLenD0, ptD0); + registry.fill(HIST("hDecLengthXyD0"), decLenXyD0, ptD0); + registry.fill(HIST("hCpaD0"), cpaD0, ptD0); + registry.fill(HIST("hCpaXyD0"), cpaXyD0, ptD0); + if constexpr (withDmesMl) { + registry.fill(HIST("hMlScoreBkgD"), ptD0, candidate.prong0MlScoreBkg()); + registry.fill(HIST("hMlScorePromptD"), ptD0, candidate.prong0MlScorePrompt()); + registry.fill(HIST("hMlScoreNonPromptD"), ptD0, candidate.prong0MlScoreNonprompt()); + } + if constexpr (withBplusMl) { + registry.fill(HIST("hMlScoreSigBplus"), ptCandBplus, candidate.mlProbBplusToD0Pi()); + } + } + } + if (fillSparses) { + if constexpr (doMc) { + if (isSignal) { + if constexpr (withDmesMl) { + registry.fill(HIST("hMassPtCutVarsRecSig"), invMassBplus, ptCandBplus, candidate.decayLength(), candidate.decayLengthXY() / candidate.errorDecayLengthXY(), candidate.impactParameterProduct(), candidate.cpa(), invMassD0, ptD0, candidate.prong0MlScoreBkg(), candidate.prong0MlScoreNonprompt()); + } else { + registry.fill(HIST("hMassPtCutVarsRecSig"), invMassBplus, ptCandBplus, candidate.decayLength(), candidate.decayLengthXY() / candidate.errorDecayLengthXY(), candidate.impactParameterProduct(), candidate.cpa(), invMassD0, ptD0, decLenD0, cpaD0); + } + } else if (fillBackground) { + if constexpr (withDmesMl) { + registry.fill(HIST("hMassPtCutVarsRecBg"), invMassBplus, ptCandBplus, candidate.decayLength(), candidate.decayLengthXY() / candidate.errorDecayLengthXY(), candidate.impactParameterProduct(), candidate.cpa(), invMassD0, ptD0, candidate.prong0MlScoreBkg(), candidate.prong0MlScoreNonprompt()); + } else { + registry.fill(HIST("hMassPtCutVarsRecBg"), invMassBplus, ptCandBplus, candidate.decayLength(), candidate.decayLengthXY() / candidate.errorDecayLengthXY(), candidate.impactParameterProduct(), candidate.cpa(), invMassD0, ptD0, decLenD0, cpaD0); + } + } + } else { + if constexpr (withDmesMl) { + registry.fill(HIST("hMassPtCutVars"), invMassBplus, ptCandBplus, candidate.decayLength(), candidate.decayLengthXY() / candidate.errorDecayLengthXY(), candidate.impactParameterProduct(), candidate.cpa(), invMassD0, ptD0, candidate.prong0MlScoreBkg(), candidate.prong0MlScoreNonprompt()); + } else { + registry.fill(HIST("hMassPtCutVars"), invMassBplus, ptCandBplus, candidate.decayLength(), candidate.decayLengthXY() / candidate.errorDecayLengthXY(), candidate.impactParameterProduct(), candidate.cpa(), invMassD0, ptD0, decLenD0, cpaD0); + } + } + } + if (fillTree) { + float pseudoRndm = ptD0 * 1000. - (int64_t)(ptD0 * 1000); + if (ptCandBplus >= ptMaxForDownSample || pseudoRndm < downSampleBkgFactor) { + float prong0MlScoreBkg = -1.; + float prong0MlScorePrompt = -1.; + float prong0MlScoreNonprompt = -1.; + float candidateMlScoreSig = -1; + if constexpr (withDmesMl) { + prong0MlScoreBkg = candidate.prong0MlScoreBkg(); + prong0MlScorePrompt = candidate.prong0MlScorePrompt(); + prong0MlScoreNonprompt = candidate.prong0MlScoreNonprompt(); + } + if constexpr (withBplusMl) { + candidateMlScoreSig = candidate.mlProbBplusToD0Pi(); + } + auto prong1 = candidate.template prong1_as(); + + float ptMother = -1.; + if constexpr (doMc) { + ptMother = candidate.ptMother(); + } + + hfRedCandBpLite( + candidate.chi2PCA(), + candidate.decayLength(), + candidate.decayLengthXY(), + candidate.decayLengthNormalised(), + candidate.decayLengthXYNormalised(), + invMassD0, + ptD0, + candidate.ptProng1(), + candidate.impactParameter0(), + candidate.impactParameter1(), + candidate.impactParameterProduct(), + prong1.tpcNSigmaPi(), + prong1.tofNSigmaPi(), + prong0MlScoreBkg, + prong0MlScorePrompt, + prong0MlScoreNonprompt, + candidateMlScoreSig, + candidate.isSelBplusToD0Pi(), + invMassBplus, + ptCandBplus, + candidate.cpa(), + candidate.cpaXY(), + candidate.maxNormalisedDeltaIP(), + candidate.eta(), + candidate.phi(), + hfHelper.yBplus(candidate), + flagMcMatchRec, + isSignal, + ptMother); + } + if constexpr (withDecayTypeCheck) { + float candidateMlScoreSig = -1; + if constexpr (withBplusMl) { + candidateMlScoreSig = candidate.mlProbBplusToD0Pi(); + } + hfRedBpMcCheck( + flagMcMatchRec, + invMassD0, + ptD0, + invMassBplus, + ptCandBplus, + candidateMlScoreSig, + candidate.pdgCodeBeautyMother(), + candidate.pdgCodeProng0(), + candidate.pdgCodeProng1(), + candidate.pdgCodeProng2()); + } + } + } + + /// Fill particle histograms (gen MC truth) + void fillCandMcGen(aod::HfMcGenRedBps::iterator const& particle) + { + auto ptParticle = particle.ptTrack(); + auto yParticle = particle.yTrack(); + auto etaParticle = particle.etaTrack(); + if (yCandGenMax >= 0. && std::abs(yParticle) > yCandGenMax) { + return; + } + std::array ptProngs = {particle.ptProng0(), particle.ptProng1()}; + std::array yProngs = {particle.yProng0(), particle.yProng1()}; + std::array etaProngs = {particle.etaProng0(), particle.etaProng1()}; + bool prongsInAcc = isProngInAcceptance(etaProngs[0], ptProngs[0]) && isProngInAcceptance(etaProngs[1], ptProngs[1]); + + if (fillHistograms) { + registry.fill(HIST("hPtProng0Gen"), ptParticle, ptProngs[0]); + registry.fill(HIST("hPtProng1Gen"), ptParticle, ptProngs[1]); + registry.fill(HIST("hYProng0Gen"), ptParticle, yProngs[0]); + registry.fill(HIST("hYProng1Gen"), ptParticle, yProngs[1]); + registry.fill(HIST("hEtaProng0Gen"), ptParticle, etaProngs[0]); + registry.fill(HIST("hEtaProng1Gen"), ptParticle, etaProngs[1]); + + registry.fill(HIST("hYGen"), ptParticle, yParticle); + registry.fill(HIST("hEtaGen"), ptParticle, etaParticle); + + // generated B+ with daughters in geometrical acceptance + if (prongsInAcc) { + registry.fill(HIST("hYGenWithProngsInAcceptance"), ptParticle, yParticle); + } + } + if (fillSparses) { + registry.fill(HIST("hPtYGenSig"), ptParticle, yParticle); + if (prongsInAcc) { + registry.fill(HIST("hPtYWithProngsInAccepanceGenSig"), ptParticle, yParticle); + } + } + } + + // Process functions + void processData(soa::Filtered> const& candidates, + aod::HfRed2Prongs const& candidatesD, + TracksPion const&) { for (const auto& candidate : candidates) { - if (!TESTBIT(candidate.hfflag(), hf_cand_bplus::DecayType::BplusToD0Pi)) { + if (yCandRecoMax >= 0. && std::abs(hfHelper.yBplus(candidate)) > yCandRecoMax) { continue; } + fillCand(candidate, candidatesD); + } // candidate loop + } // processData + PROCESS_SWITCH(HfTaskBplusReduced, processData, "Process data without ML scores for D0 daughter", true); + + void processDataWithDmesMl(soa::Filtered> const& candidates, + aod::HfRed2Prongs const& candidatesD, + TracksPion const&) + { + for (const auto& candidate : candidates) { if (yCandRecoMax >= 0. && std::abs(hfHelper.yBplus(candidate)) > yCandRecoMax) { continue; } + fillCand(candidate, candidatesD); + } // candidate loop + } // processDataWithDmesMl + PROCESS_SWITCH(HfTaskBplusReduced, processDataWithDmesMl, "Process data with ML scores for D0 daughter", false); - auto candD0 = candidate.prong0_as(); - auto candPi = candidate.prong1_as(); - auto ptCandBplus = candidate.pt(); - auto invMassD0 = (candPi.signed1Pt() < 0) ? candD0.invMassD0() : candD0.invMassD0Bar(); - - registry.fill(HIST("hMass"), hfHelper.invMassBplusToD0Pi(candidate), ptCandBplus); - registry.fill(HIST("hPtCand"), ptCandBplus); - registry.fill(HIST("hPtProng0"), candidate.ptProng0()); - registry.fill(HIST("hPtProng1"), candidate.ptProng1()); - registry.fill(HIST("hd0d0"), candidate.impactParameterProduct(), ptCandBplus); - registry.fill(HIST("hDecLength"), candidate.decayLength(), ptCandBplus); - registry.fill(HIST("hDecLengthXY"), candidate.decayLengthXY(), ptCandBplus); - registry.fill(HIST("hd0Prong0"), candidate.impactParameter0(), ptCandBplus); - registry.fill(HIST("hd0Prong1"), candidate.impactParameter1(), ptCandBplus); - registry.fill(HIST("hCPA"), candidate.cpa(), ptCandBplus); - registry.fill(HIST("hCPAxy"), candidate.cpaXY(), ptCandBplus); - registry.fill(HIST("hEta"), candidate.eta(), ptCandBplus); - registry.fill(HIST("hRapidity"), hfHelper.yBplus(candidate), ptCandBplus); - registry.fill(HIST("hImpParErr"), candidate.errorImpactParameter0(), ptCandBplus); - registry.fill(HIST("hImpParErr"), candidate.errorImpactParameter1(), ptCandBplus); - registry.fill(HIST("hDecLenErr"), candidate.errorDecayLength(), ptCandBplus); - registry.fill(HIST("hDecLenXYErr"), candidate.errorDecayLengthXY(), ptCandBplus); - registry.fill(HIST("hInvMassD0"), invMassD0, ptCandBplus); - registry.fill(HIST("hCPAFinerBinning"), candidate.cpa(), ptCandBplus); - registry.fill(HIST("hCPAxyFinerBinning"), candidate.cpaXY(), ptCandBplus); + void processDataWithBplusMl(soa::Filtered> const& candidates, + aod::HfRed2Prongs const& candidatesD, + TracksPion const&) + { + for (const auto& candidate : candidates) { + if (yCandRecoMax >= 0. && std::abs(hfHelper.yBplus(candidate)) > yCandRecoMax) { + continue; + } + fillCand(candidate, candidatesD); } // candidate loop - } // process + } // processDataWithBplusMl + PROCESS_SWITCH(HfTaskBplusReduced, processDataWithBplusMl, "Process data with(out) ML scores for B+ (D0 daughter)", false); - /// B+ MC analysis and fill histograms - void processMc(soa::Join const& candidates, + void processMc(soa::Join const& candidates, aod::HfMcGenRedBps const& mcParticles, - aod::HfRed2Prongs const&) + aod::HfRed2Prongs const& candidatesD, + TracksPion const&) { // MC rec for (const auto& candidate : candidates) { - if (!TESTBIT(candidate.hfflag(), hf_cand_bplus::DecayType::BplusToD0Pi)) { + if (yCandRecoMax >= 0. && std::abs(hfHelper.yBplus(candidate)) > yCandRecoMax) { continue; } + fillCand(candidate, candidatesD); + } // rec + + // MC gen. level + for (const auto& particle : mcParticles) { + fillCandMcGen(particle); + } // gen + } // processMc + PROCESS_SWITCH(HfTaskBplusReduced, processMc, "Process MC without ML scores for B+ and D0 daughter", false); + + void processMcWithDecayTypeCheck(soa::Filtered> const& candidates, + aod::HfMcGenRedBps const& mcParticles, + aod::HfRed2Prongs const& candidatesD, + TracksPion const&) + { + // MC rec + for (const auto& candidate : candidates) { if (yCandRecoMax >= 0. && std::abs(hfHelper.yBplus(candidate)) > yCandRecoMax) { continue; } + fillCand(candidate, candidatesD); + } // rec - auto ptCandBplus = candidate.pt(); - auto candD0 = candidate.prong0_as(); - std::array posPv{candidate.posX(), candidate.posY(), candidate.posZ()}; - std::array posSvD{candD0.xSecondaryVertex(), candD0.ySecondaryVertex(), candD0.zSecondaryVertex()}; - std::array momD{candD0.pVector()}; - auto cospD0 = RecoDecay::cpa(posPv, posSvD, momD); - auto decLenD0 = RecoDecay::distance(posPv, posSvD); - - if (TESTBIT(std::abs(candidate.flagMcMatchRec()), hf_cand_bplus::DecayType::BplusToD0Pi)) { - - registry.fill(HIST("hPtGenSig"), candidate.ptMother()); - registry.fill(HIST("hPtRecSig"), ptCandBplus); - registry.fill(HIST("hCPARecSig"), candidate.cpa(), ptCandBplus); - registry.fill(HIST("hCPAxyRecSig"), candidate.cpaXY(), ptCandBplus); - registry.fill(HIST("hCPAFinerBinningRecSig"), candidate.cpa(), ptCandBplus); - registry.fill(HIST("hCPAxyFinerBinningRecSig"), candidate.cpaXY(), ptCandBplus); - registry.fill(HIST("hEtaRecSig"), candidate.eta(), ptCandBplus); - registry.fill(HIST("hRapidityRecSig"), hfHelper.yBplus(candidate), ptCandBplus); - registry.fill(HIST("hDecLengthRecSig"), candidate.decayLength(), ptCandBplus); - registry.fill(HIST("hDecLengthXYRecSig"), candidate.decayLengthXY(), ptCandBplus); - registry.fill(HIST("hMassRecSig"), hfHelper.invMassBplusToD0Pi(candidate), ptCandBplus); - registry.fill(HIST("hd0Prong0RecSig"), candidate.impactParameter0(), ptCandBplus); - registry.fill(HIST("hd0Prong1RecSig"), candidate.impactParameter1(), ptCandBplus); - registry.fill(HIST("hPtProng0RecSig"), candidate.ptProng0(), ptCandBplus); - registry.fill(HIST("hPtProng1RecSig"), candidate.ptProng1(), ptCandBplus); - registry.fill(HIST("hd0d0RecSig"), candidate.impactParameterProduct(), ptCandBplus); - registry.fill(HIST("hDecLengthNormRecSig"), candidate.decayLengthXYNormalised(), ptCandBplus); - registry.fill(HIST("hCPAD0RecSig"), cospD0, candidate.ptProng0()); - registry.fill(HIST("hDecLengthD0RecSig"), decLenD0, candidate.ptProng0()); - } else { - registry.fill(HIST("hPtRecBg"), ptCandBplus); - registry.fill(HIST("hCPARecBg"), candidate.cpa(), ptCandBplus); - registry.fill(HIST("hCPAxyRecBg"), candidate.cpaXY(), ptCandBplus); - registry.fill(HIST("hCPAFinerBinningRecBg"), candidate.cpa(), ptCandBplus); - registry.fill(HIST("hCPAxyFinerBinningRecBg"), candidate.cpaXY(), ptCandBplus); - registry.fill(HIST("hEtaRecBg"), candidate.eta(), ptCandBplus); - registry.fill(HIST("hRapidityRecBg"), hfHelper.yBplus(candidate), ptCandBplus); - registry.fill(HIST("hDecLengthRecBg"), candidate.decayLength(), ptCandBplus); - registry.fill(HIST("hDecLengthXYRecBg"), candidate.decayLengthXY(), ptCandBplus); - registry.fill(HIST("hMassRecBg"), hfHelper.invMassBplusToD0Pi(candidate), ptCandBplus); - registry.fill(HIST("hd0Prong0RecBg"), candidate.impactParameter0(), ptCandBplus); - registry.fill(HIST("hd0Prong1RecBg"), candidate.impactParameter1(), ptCandBplus); - registry.fill(HIST("hPtProng0RecBg"), candidate.ptProng0(), ptCandBplus); - registry.fill(HIST("hPtProng1RecBg"), candidate.ptProng1(), ptCandBplus); - registry.fill(HIST("hd0d0RecBg"), candidate.impactParameterProduct(), ptCandBplus); - registry.fill(HIST("hDecLengthNormRecBg"), candidate.decayLengthXYNormalised(), ptCandBplus); - registry.fill(HIST("hCPAD0RecBg"), cospD0, candidate.ptProng0()); - registry.fill(HIST("hDecLengthD0RecBg"), decLenD0, candidate.ptProng0()); + // MC gen. level + for (const auto& particle : mcParticles) { + fillCandMcGen(particle); + } // gen + } // processMc + PROCESS_SWITCH(HfTaskBplusReduced, processMcWithDecayTypeCheck, "Process MC with decay type check and without ML scores for B+ and D daughter", false); + + void processMcWithDmesMl(soa::Join const& candidates, + aod::HfMcGenRedBps const& mcParticles, + aod::HfRed2Prongs const& candidatesD, + TracksPion const&) + { + // MC rec + for (const auto& candidate : candidates) { + if (yCandRecoMax >= 0. && std::abs(hfHelper.yBplus(candidate)) > yCandRecoMax) { + continue; } + fillCand(candidate, candidatesD); } // rec // MC gen. level for (const auto& particle : mcParticles) { - auto ptParticle = particle.ptTrack(); - auto yParticle = particle.yTrack(); - auto etaParticle = particle.etaTrack(); - if (yCandGenMax >= 0. && std::abs(yParticle) > yCandGenMax) { + fillCandMcGen(particle); + } // gen + } // processMcWithDmesMl + PROCESS_SWITCH(HfTaskBplusReduced, processMcWithDmesMl, "Process MC with(out) ML scores for D0 daughter (B+)", false); + + void processMcWithBplusMl(soa::Filtered> const& candidates, + aod::HfMcGenRedBps const& mcParticles, + aod::HfRed2Prongs const& candidatesD, + TracksPion const&) + { + // MC rec + for (const auto& candidate : candidates) { + if (yCandRecoMax >= 0. && std::abs(hfHelper.yBplus(candidate)) > yCandRecoMax) { continue; } + fillCand(candidate, candidatesD); + } // rec - std::array ptProngs = {particle.ptProng0(), particle.ptProng1()}; - std::array yProngs = {particle.yProng0(), particle.yProng1()}; - std::array etaProngs = {particle.etaProng0(), particle.etaProng1()}; - - registry.fill(HIST("hPtProng0Gen"), ptProngs[0], ptParticle); - registry.fill(HIST("hPtProng1Gen"), ptProngs[1], ptParticle); - registry.fill(HIST("hPtProngsVsPtBGen"), ptProngs[0], ptProngs[1], ptParticle); - registry.fill(HIST("hYProng0Gen"), yProngs[0], ptParticle); - registry.fill(HIST("hYProng1Gen"), yProngs[1], ptParticle); - registry.fill(HIST("hYProngsVsBplusGen"), yProngs[0], yProngs[1], ptParticle); - registry.fill(HIST("hEtaProng0Gen"), etaProngs[0], ptParticle); - registry.fill(HIST("hEtaProng1Gen"), etaProngs[1], ptParticle); - registry.fill(HIST("hEtaProngsVsBplusGen"), etaProngs[0], etaProngs[1], ptParticle); - - registry.fill(HIST("hPtGen"), ptParticle); - registry.fill(HIST("hYGen"), yParticle, ptParticle); - registry.fill(HIST("hEtaGen"), etaParticle, ptParticle); - - // generated B+ with |y|<0.5 - if (std::abs(yParticle) < 0.5) { - registry.fill(HIST("hPtGenWithRapidityBelowHalf"), ptParticle); - } + // MC gen. level + for (const auto& particle : mcParticles) { + fillCandMcGen(particle); + } // gen + } // processMcWithBplusMl + PROCESS_SWITCH(HfTaskBplusReduced, processMcWithBplusMl, "Process MC with(out) ML scores for B+ (D0 daughter)", false); - // generated B+ with daughters in geometrical acceptance - if (isProngInAcceptance(etaProngs[0], ptProngs[0]) && isProngInAcceptance(etaProngs[1], ptProngs[1])) { - registry.fill(HIST("hPtGenWithProngsInAcceptance"), ptParticle); - registry.fill(HIST("hYGenWithProngsInAcceptance"), yParticle, ptParticle); - registry.fill(HIST("hEtaGenWithProngsInAcceptance"), etaParticle, ptParticle); + void processMcWithBplusMlAndDecayTypeCheck(soa::Filtered> const& candidates, + aod::HfMcGenRedBps const& mcParticles, + aod::HfRed2Prongs const& candidatesD, + TracksPion const&) + { + // MC rec + for (const auto& candidate : candidates) { + if (yCandRecoMax >= 0. && std::abs(hfHelper.yBplus(candidate)) > yCandRecoMax) { + continue; } + fillCand(candidate, candidatesD); + } // rec + + // MC gen. level + for (const auto& particle : mcParticles) { + fillCandMcGen(particle); } // gen - } // process - PROCESS_SWITCH(HfTaskBplusReduced, processMc, "Process MC", false); + } // processMc + PROCESS_SWITCH(HfTaskBplusReduced, processMcWithBplusMlAndDecayTypeCheck, "Process MC with decay type check and with(out) ML scores for B+ (D0 daughter)", false); }; // struct WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) diff --git a/PWGHF/DataModel/CandidateReconstructionTables.h b/PWGHF/DataModel/CandidateReconstructionTables.h index a232d83ef7c..5b21aa5c1b4 100644 --- a/PWGHF/DataModel/CandidateReconstructionTables.h +++ b/PWGHF/DataModel/CandidateReconstructionTables.h @@ -767,6 +767,11 @@ DECLARE_SOA_COLUMN(OriginMcGen, originMcGen, int8_t); // particle origin, DECLARE_SOA_COLUMN(DebugMcRec, debugMcRec, int8_t); // debug flag for mis-association reconstruction level enum DecayType { BplusToD0Pi = 0 }; + +enum DecayTypeMc : uint8_t { BplusToD0PiToKPiPi = 0, + PartlyRecoDecay, + OtherDecay, + NDecayTypeMc }; } // namespace hf_cand_bplus // declare dedicated BPlus decay candidate table @@ -779,7 +784,6 @@ DECLARE_SOA_TABLE(HfCandBplusBase, "AOD", "HFCANDBPLUSBASE", hf_cand::PxProng1, hf_cand::PyProng1, hf_cand::PzProng1, hf_cand::ImpactParameter0, hf_cand::ImpactParameter1, hf_cand::ErrorImpactParameter0, hf_cand::ErrorImpactParameter1, - hf_track_index::HFflag, /* dynamic columns */ hf_cand_2prong::M, hf_cand_2prong::M2, @@ -1473,7 +1477,6 @@ DECLARE_SOA_TABLE(HfCandB0Base, "AOD", "HFCANDB0BASE", hf_cand::PxProng1, hf_cand::PyProng1, hf_cand::PzProng1, hf_cand::ImpactParameter0, hf_cand::ImpactParameter1, hf_cand::ErrorImpactParameter0, hf_cand::ErrorImpactParameter1, - hf_track_index::HFflag, /* dynamic columns */ hf_cand_2prong::M, hf_cand_2prong::M2, diff --git a/PWGHF/DataModel/CandidateSelectionTables.h b/PWGHF/DataModel/CandidateSelectionTables.h index 4776c1067d5..36676d8ecf0 100644 --- a/PWGHF/DataModel/CandidateSelectionTables.h +++ b/PWGHF/DataModel/CandidateSelectionTables.h @@ -263,12 +263,16 @@ DECLARE_SOA_TABLE(HfMlBsToDsPi, "AOD", "HFMLBS", //! namespace hf_sel_candidate_bplus { -DECLARE_SOA_COLUMN(IsSelBplusToD0Pi, isSelBplusToD0Pi, int); //! +DECLARE_SOA_COLUMN(IsSelBplusToD0Pi, isSelBplusToD0Pi, int); //! selection flag on B+ candidate +DECLARE_SOA_COLUMN(MlProbBplusToD0Pi, mlProbBplusToD0Pi, float); //! ML score of B+ candidate for signal class } // namespace hf_sel_candidate_bplus DECLARE_SOA_TABLE(HfSelBplusToD0Pi, "AOD", "HFSELBPLUS", //! hf_sel_candidate_bplus::IsSelBplusToD0Pi); +DECLARE_SOA_TABLE(HfMlBplusToD0Pi, "AOD", "HFMLBPLUS", //! + hf_sel_candidate_bplus::MlProbBplusToD0Pi); + namespace hf_sel_candidate_lb { DECLARE_SOA_COLUMN(IsSelLbToLcPi, isSelLbToLcPi, int); //! diff --git a/PWGHF/TableProducer/candidateCreatorB0.cxx b/PWGHF/TableProducer/candidateCreatorB0.cxx index 7fea6b97b20..ec298631266 100644 --- a/PWGHF/TableProducer/candidateCreatorB0.cxx +++ b/PWGHF/TableProducer/candidateCreatorB0.cxx @@ -342,8 +342,6 @@ struct HfCandidateCreatorB0 { auto errorDecayLength = std::sqrt(getRotatedCovMatrixXX(covMatrixPV, phi, theta) + getRotatedCovMatrixXX(covMatrixPCA, phi, theta)); auto errorDecayLengthXY = std::sqrt(getRotatedCovMatrixXX(covMatrixPV, phi, 0.) + getRotatedCovMatrixXX(covMatrixPCA, phi, 0.)); - int hfFlag = BIT(hf_cand_b0::DecayType::B0ToDPi); - // fill the candidate table for the B0 here: rowCandidateBase(thisCollId, collision.posX(), collision.posY(), collision.posZ(), @@ -353,8 +351,7 @@ struct HfCandidateCreatorB0 { pVecD[0], pVecD[1], pVecD[2], pVecPion[0], pVecPion[1], pVecPion[2], dcaD.getY(), dcaPion.getY(), - std::sqrt(dcaD.getSigmaY2()), std::sqrt(dcaPion.getSigmaY2()), - hfFlag); + std::sqrt(dcaD.getSigmaY2()), std::sqrt(dcaPion.getSigmaY2())); rowCandidateProngs(candD.globalIndex(), trackPion.globalIndex()); } // pi loop diff --git a/PWGHF/TableProducer/candidateCreatorBplus.cxx b/PWGHF/TableProducer/candidateCreatorBplus.cxx index c780cf972e2..695f42d4710 100644 --- a/PWGHF/TableProducer/candidateCreatorBplus.cxx +++ b/PWGHF/TableProducer/candidateCreatorBplus.cxx @@ -335,8 +335,6 @@ struct HfCandidateCreatorBplus { auto errorDecayLength = std::sqrt(getRotatedCovMatrixXX(covMatrixPV, phi, theta) + getRotatedCovMatrixXX(covMatrixPCA, phi, theta)); auto errorDecayLengthXY = std::sqrt(getRotatedCovMatrixXX(covMatrixPV, phi, 0.) + getRotatedCovMatrixXX(covMatrixPCA, phi, 0.)); - int hfFlag = BIT(hf_cand_bplus::DecayType::BplusToD0Pi); - // compute invariant mass square and apply selection auto invMass2D0Pi = RecoDecay::m2(std::array{pVecD0, pVecBach}, std::array{massD0, massPi}); if ((invMass2D0Pi < invMass2D0PiMin) || (invMass2D0Pi > invMass2D0PiMax)) { @@ -353,8 +351,7 @@ struct HfCandidateCreatorBplus { pVecD0[0], pVecD0[1], pVecD0[2], pVecBach[0], pVecBach[1], pVecBach[2], impactParameter0.getY(), impactParameter1.getY(), - std::sqrt(impactParameter0.getSigmaY2()), std::sqrt(impactParameter1.getSigmaY2()), - hfFlag); + std::sqrt(impactParameter0.getSigmaY2()), std::sqrt(impactParameter1.getSigmaY2())); rowCandidateProngs(candD0.globalIndex(), trackPion.globalIndex()); // index D0 and bachelor } // track loop diff --git a/PWGHF/TableProducer/candidateSelectorB0ToDPi.cxx b/PWGHF/TableProducer/candidateSelectorB0ToDPi.cxx index c256265b341..985be6f5aa6 100644 --- a/PWGHF/TableProducer/candidateSelectorB0ToDPi.cxx +++ b/PWGHF/TableProducer/candidateSelectorB0ToDPi.cxx @@ -118,15 +118,6 @@ struct HfCandidateSelectorB0ToDPi { int statusB0ToDPi = 0; auto ptCandB0 = hfCandB0.pt(); - // check if flagged as B0 → D π - if (!TESTBIT(hfCandB0.hfflag(), hf_cand_b0::DecayType::B0ToDPi)) { - hfSelB0ToDPiCandidate(statusB0ToDPi); - if (activateQA) { - registry.fill(HIST("hSelections"), 1, ptCandB0); - } - // LOGF(info, "B0 candidate selection failed at hfflag check"); - continue; - } SETBIT(statusB0ToDPi, SelectionStep::RecoSkims); // RecoSkims = 0 --> statusB0ToDPi = 1 if (activateQA) { registry.fill(HIST("hSelections"), 2 + SelectionStep::RecoSkims, ptCandB0); diff --git a/PWGHF/TableProducer/candidateSelectorBplusToD0Pi.cxx b/PWGHF/TableProducer/candidateSelectorBplusToD0Pi.cxx index a5ca2557962..e8589833f9e 100644 --- a/PWGHF/TableProducer/candidateSelectorBplusToD0Pi.cxx +++ b/PWGHF/TableProducer/candidateSelectorBplusToD0Pi.cxx @@ -142,12 +142,6 @@ struct HfCandidateSelectorBplusToD0Pi { int statusBplus = 0; auto ptCandB = hfCandB.pt(); - // check if flagged as B+ --> D0bar Pi - if (!(hfCandB.hfflag() & 1 << hf_cand_bplus::DecayType::BplusToD0Pi)) { - hfSelBplusToD0PiCandidate(statusBplus); - // LOGF(debug, "B+ candidate selection failed at hfflag check"); - continue; - } SETBIT(statusBplus, SelectionStep::RecoSkims); // RecoSkims = 0 --> statusBplus = 1 if (activateQA) { registry.fill(HIST("hSelections"), 2 + SelectionStep::RecoSkims, ptCandB); diff --git a/PWGHF/TableProducer/treeCreatorBplusToD0Pi.cxx b/PWGHF/TableProducer/treeCreatorBplusToD0Pi.cxx index 64c7a8990c8..fdc18d3b3e0 100644 --- a/PWGHF/TableProducer/treeCreatorBplusToD0Pi.cxx +++ b/PWGHF/TableProducer/treeCreatorBplusToD0Pi.cxx @@ -26,42 +26,38 @@ #include "PWGHF/DataModel/CandidateSelectionTables.h" using namespace o2; -using namespace o2::aod; using namespace o2::framework; +using namespace o2::framework::expressions; namespace o2::aod { namespace full { -DECLARE_SOA_COLUMN(RSecondaryVertex, rSecondaryVertex, float); -DECLARE_SOA_COLUMN(PtProng0, ptProng0, float); -DECLARE_SOA_COLUMN(PProng0, pProng0, float); -DECLARE_SOA_COLUMN(PtProng1, ptProng1, float); -DECLARE_SOA_COLUMN(PProng1, pProng1, float); -// DECLARE_SOA_COLUMN(CandidateSelFlag, candidateSelFlag, int8_t); -DECLARE_SOA_COLUMN(M, m, float); -DECLARE_SOA_COLUMN(Pt, pt, float); -DECLARE_SOA_COLUMN(P, p, float); -DECLARE_SOA_COLUMN(Eta, eta, float); -DECLARE_SOA_COLUMN(Phi, phi, float); -DECLARE_SOA_COLUMN(Y, y, float); -DECLARE_SOA_COLUMN(DecayLength, decayLength, float); -DECLARE_SOA_COLUMN(DecayLengthXY, decayLengthXY, float); -DECLARE_SOA_COLUMN(DecayLengthNormalised, decayLengthNormalised, float); -DECLARE_SOA_COLUMN(DecayLengthXYNormalised, decayLengthXYNormalised, float); -DECLARE_SOA_COLUMN(CPA, cpa, float); -DECLARE_SOA_COLUMN(CPAXY, cpaXY, float); -DECLARE_SOA_COLUMN(ImpactParameterProduct, impactParameterProduct, float); -DECLARE_SOA_COLUMN(ImpactParameter0, impactParameter0, float); -DECLARE_SOA_COLUMN(ImpactParameter1, impactParameter1, float); -DECLARE_SOA_COLUMN(ImpactParameterNormalised0, impactParameterNormalised0, float); -DECLARE_SOA_COLUMN(ImpactParameterNormalised1, impactParameterNormalised1, float); +DECLARE_SOA_COLUMN(RSecondaryVertex, rSecondaryVertex, float); //! Radius of secondary vertex (cm) +DECLARE_SOA_COLUMN(PtProng0, ptProng0, float); //! Transverse momentum of prong0 (GeV/c) +DECLARE_SOA_COLUMN(PProng0, pProng0, float); //! Momentum of prong0 (GeV/c) +DECLARE_SOA_COLUMN(ImpactParameterNormalised0, impactParameterNormalised0, float); //! Normalised impact parameter of prong0 +DECLARE_SOA_COLUMN(PtProng1, ptProng1, float); //! Transverse momentum of prong1 (GeV/c) +DECLARE_SOA_COLUMN(PProng1, pProng1, float); //! Momentum of prong1 (in GeV/c) +DECLARE_SOA_COLUMN(ImpactParameterNormalised1, impactParameterNormalised1, float); //! Normalised impact parameter of prong1 +DECLARE_SOA_COLUMN(CandidateSelFlag, candidateSelFlag, int); //! Selection flag of candidate (output of candidateSelector) +DECLARE_SOA_COLUMN(M, m, float); //! Invariant mass of candidate (GeV/c2) +DECLARE_SOA_COLUMN(Pt, pt, float); //! Transverse momentum of candidate (GeV/c) +DECLARE_SOA_COLUMN(P, p, float); //! Momentum of candidate (GeV/c) +DECLARE_SOA_COLUMN(Y, y, float); //! Rapidity of candidate +DECLARE_SOA_COLUMN(Eta, eta, float); //! Pseudorapidity of candidate +DECLARE_SOA_COLUMN(Phi, phi, float); //! Azimuth angle of candidate +DECLARE_SOA_COLUMN(E, e, float); //! Energy of candidate (GeV) +DECLARE_SOA_COLUMN(NSigTpcPi1, nSigTpcPi1, float); //! TPC Nsigma separation for prong1 with pion mass hypothesis +DECLARE_SOA_COLUMN(NSigTofPi1, nSigTofPi1, float); //! TOF Nsigma separation for prong1 with pion mass hypothesis +DECLARE_SOA_COLUMN(DecayLength, decayLength, float); //! Decay length of candidate (cm) +DECLARE_SOA_COLUMN(DecayLengthXY, decayLengthXY, float); //! Transverse decay length of candidate (cm) +DECLARE_SOA_COLUMN(DecayLengthNormalised, decayLengthNormalised, float); //! Normalised decay length of candidate +DECLARE_SOA_COLUMN(DecayLengthXYNormalised, decayLengthXYNormalised, float); //! Normalised transverse decay length of candidate +DECLARE_SOA_COLUMN(Cpa, cpa, float); //! Cosine pointing angle of candidate +DECLARE_SOA_COLUMN(CpaXY, cpaXY, float); //! Cosine pointing angle of candidate in transverse plane +DECLARE_SOA_COLUMN(MaxNormalisedDeltaIP, maxNormalisedDeltaIP, float); //! Maximum normalized difference between measured and expected impact parameter of candidate prongs DECLARE_SOA_COLUMN(Ct, ct, float); -DECLARE_SOA_COLUMN(Chi2PCA, chi2PCA, float); -DECLARE_SOA_COLUMN(NSigmaTOFBachPi, nSigmaTOFBachPi, float); -DECLARE_SOA_COLUMN(NSigmaTOFBachKa, nSigmaTOFBachKa, float); -DECLARE_SOA_COLUMN(NSigmaTPCBachPi, nSigmaTPCBachPi, float); -DECLARE_SOA_COLUMN(NSigmaTPCBachKa, nSigmaTPCBachKa, float); DECLARE_SOA_COLUMN(MCflag, mcflag, int8_t); // D0 (Prong0) selection variable DECLARE_SOA_COLUMN(D0M, d0M, float); @@ -93,49 +89,57 @@ DECLARE_SOA_COLUMN(NSigmaTPCTrk1Pi, nSigmaTPCTrk1Pi, float); // Events DECLARE_SOA_COLUMN(IsEventReject, isEventReject, int); DECLARE_SOA_COLUMN(RunNumber, runNumber, int); -DECLARE_SOA_INDEX_COLUMN_FULL(Candidate, candidate, int, HfCandBplus, "_0"); -DECLARE_SOA_INDEX_COLUMN(McParticle, mcParticle); -DECLARE_SOA_INDEX_COLUMN(McCollision, mcCollision); } // namespace full // put the arguments into the table DECLARE_SOA_TABLE(HfCandBpFulls, "AOD", "HFCANDBPFULL", + collision::BCId, + collision::NumContrib, + collision::PosX, + collision::PosY, + collision::PosZ, + hf_cand::XSecondaryVertex, + hf_cand::YSecondaryVertex, + hf_cand::ZSecondaryVertex, + hf_cand::ErrorDecayLength, + hf_cand::ErrorDecayLengthXY, + hf_cand::Chi2PCA, full::RSecondaryVertex, - full::PtProng0, - full::PProng0, - full::PtProng1, - full::PProng1, - // full::CandidateSelFlag, - full::M, - full::Pt, - full::P, - full::Ct, - full::Eta, - full::Phi, - full::Y, full::DecayLength, full::DecayLengthXY, full::DecayLengthNormalised, full::DecayLengthXYNormalised, - full::CPA, - full::CPAXY, - full::ImpactParameterProduct, - hf_cand::ImpactParameter0, - hf_cand::ImpactParameter1, full::ImpactParameterNormalised0, + full::PtProng0, + full::PProng0, full::ImpactParameterNormalised1, + full::PtProng1, + full::PProng1, hf_cand::PxProng0, hf_cand::PyProng0, hf_cand::PzProng0, hf_cand::PxProng1, hf_cand::PyProng1, hf_cand::PzProng1, - hf_cand::Chi2PCA, - full::NSigmaTOFBachPi, - full::NSigmaTOFBachKa, - full::NSigmaTPCBachPi, - full::NSigmaTPCBachKa, - full::MCflag, + hf_cand::ImpactParameter0, + hf_cand::ImpactParameter1, + hf_cand::ErrorImpactParameter0, + hf_cand::ErrorImpactParameter1, + full::NSigTpcPi1, + full::NSigTofPi1, + full::CandidateSelFlag, + full::M, + full::Pt, + full::P, + full::Cpa, + full::CpaXY, + full::MaxNormalisedDeltaIP, + full::Ct, + full::Eta, + full::Phi, + full::Y, + full::E, + hf_cand_2prong::FlagMcMatchRec, full::D0M, full::D0PtProng0, full::D0PtProng1, @@ -160,8 +164,7 @@ DECLARE_SOA_TABLE(HfCandBpFulls, "AOD", "HFCANDBPFULL", full::NSigmaTOFTrk1Pi, full::NSigmaTOFTrk1Ka, full::NSigmaTPCTrk1Pi, - full::NSigmaTPCTrk1Ka, - full::CandidateId); + full::NSigmaTPCTrk1Ka); DECLARE_SOA_TABLE(HfCandBpFullEvs, "AOD", "HFCANDBPFULLEV", collision::BCId, @@ -174,13 +177,36 @@ DECLARE_SOA_TABLE(HfCandBpFullEvs, "AOD", "HFCANDBPFULLEV", DECLARE_SOA_TABLE(HfCandBpFullPs, "AOD", "HFCANDBPFULLP", collision::BCId, - full::McCollisionId, full::Pt, full::Eta, full::Phi, full::Y, - full::MCflag, - full::McParticleId); + hf_cand_2prong::FlagMcMatchRec, + hf_cand_2prong::OriginMcGen); + +DECLARE_SOA_TABLE(HfCandBpLites, "AOD", "HFCANDBPLITE", + hf_cand::Chi2PCA, + full::DecayLength, + full::DecayLengthXY, + full::DecayLengthNormalised, + full::DecayLengthXYNormalised, + full::PtProng0, + full::PtProng1, + hf_cand::ImpactParameter0, + hf_cand::ImpactParameter1, + full::NSigTpcPi1, + full::NSigTofPi1, + full::CandidateSelFlag, + full::M, + full::Pt, + full::Cpa, + full::CpaXY, + full::MaxNormalisedDeltaIP, + full::Eta, + full::Phi, + full::Y, + hf_cand_2prong::FlagMcMatchRec, + hf_cand_2prong::OriginMcRec); } // namespace o2::aod @@ -189,153 +215,258 @@ struct HfTreeCreatorBplusToD0Pi { Produces rowCandidateFull; Produces rowCandidateFullEvents; Produces rowCandidateFullParticles; + Produces rowCandidateLite; - Configurable isSignal{"isSignal", 1, "save only MC matched candidates"}; + Configurable selectionFlagBplus{"selectionBplus", 1, "Selection Flag for B+"}; + Configurable fillCandidateLiteTable{"fillCandidateLiteTable", false, "Switch to fill lite table with candidate properties"}; + // parameters for production of training samples + Configurable fillOnlySignal{"fillOnlySignal", false, "Flag to fill derived tables with signal for ML trainings"}; + Configurable fillOnlyBackground{"fillOnlyBackground", false, "Flag to fill derived tables with background for ML trainings"}; + Configurable downSampleBkgFactor{"downSampleBkgFactor", 1., "Fraction of background candidates to keep for ML trainings"}; + Configurable ptMaxForDownSample{"ptMaxForDownSample", 10., "Maximum pt for the application of the downsampling factor"}; HfHelper hfHelper; + using SelectedCandidatesMc = soa::Filtered>; using TracksWPid = soa::Join; + Filter filterSelectCandidates = aod::hf_sel_candidate_bplus::isSelBplusToD0Pi >= selectionFlagBplus; + + Partition recSig = nabs(aod::hf_cand_bplus::flagMcMatchRec) == (int8_t)BIT(aod::hf_cand_bplus::DecayTypeMc::BplusToD0PiToKPiPi); + Partition recBg = nabs(aod::hf_cand_bplus::flagMcMatchRec) != (int8_t)BIT(aod::hf_cand_bplus::DecayTypeMc::BplusToD0PiToKPiPi); + void init(InitContext const&) { } - void process(aod::Collisions const& collisions, - aod::McCollisions const&, - soa::Join const& candidates, - soa::Join const& particles, - TracksWPid const&, - aod::HfCand2Prong const&) + template + void fillEvent(const T& collision, int isEventReject, int runNumber) { + rowCandidateFullEvents( + collision.bcId(), + collision.numContrib(), + collision.posX(), + collision.posY(), + collision.posZ(), + isEventReject, + runNumber); + } + + template + void fillCandidateTable(const T& candidate, const U& prong1) + { + int8_t flagMc = 0; + int8_t originMc = 0; + if constexpr (doMc) { + flagMc = candidate.flagMcMatchRec(); + originMc = candidate.originMcRec(); + } + auto d0Cand = candidate.prong0(); + // adding D0 daughters to the table + auto d0Daughter0 = d0Cand.template prong0_as(); + auto d0Daughter1 = d0Cand.template prong1_as(); + auto invMassD0 = 0.; + if (prong1.signed1Pt() > 0) { + invMassD0 = hfHelper.invMassD0barToKPi(d0Cand); + } else if (prong1.signed1Pt() < 0) { + invMassD0 = hfHelper.invMassD0ToPiK(d0Cand); + } + if (fillCandidateLiteTable) { + rowCandidateLite( + candidate.chi2PCA(), + candidate.decayLength(), + candidate.decayLengthXY(), + candidate.decayLengthNormalised(), + candidate.decayLengthXYNormalised(), + candidate.ptProng0(), + candidate.ptProng1(), + candidate.impactParameter0(), + candidate.impactParameter1(), + prong1.tpcNSigmaPi(), + prong1.tofNSigmaPi(), + candidate.isSelBplusToD0Pi(), + hfHelper.invMassBplusToD0Pi(candidate), + candidate.pt(), + candidate.cpa(), + candidate.cpaXY(), + candidate.maxNormalisedDeltaIP(), + candidate.eta(), + candidate.phi(), + hfHelper.yBplus(candidate), + flagMc, + originMc); + } else { + rowCandidateFull( + prong1.collision().bcId(), + prong1.collision().numContrib(), + candidate.posX(), + candidate.posY(), + candidate.posZ(), + candidate.xSecondaryVertex(), + candidate.ySecondaryVertex(), + candidate.zSecondaryVertex(), + candidate.errorDecayLength(), + candidate.errorDecayLengthXY(), + candidate.chi2PCA(), + candidate.rSecondaryVertex(), + candidate.decayLength(), + candidate.decayLengthXY(), + candidate.decayLengthNormalised(), + candidate.decayLengthXYNormalised(), + candidate.impactParameterNormalised0(), + candidate.ptProng0(), + RecoDecay::p(candidate.pxProng0(), candidate.pyProng0(), candidate.pzProng0()), + candidate.impactParameterNormalised1(), + candidate.ptProng1(), + RecoDecay::p(candidate.pxProng1(), candidate.pyProng1(), candidate.pzProng1()), + candidate.pxProng0(), + candidate.pyProng0(), + candidate.pzProng0(), + candidate.pxProng1(), + candidate.pyProng1(), + candidate.pzProng1(), + candidate.impactParameter0(), + candidate.impactParameter1(), + candidate.errorImpactParameter0(), + candidate.errorImpactParameter1(), + prong1.tpcNSigmaPi(), + prong1.tofNSigmaPi(), + candidate.isSelBplusToD0Pi(), + hfHelper.invMassBplusToD0Pi(candidate), + candidate.pt(), + candidate.p(), + candidate.cpa(), + candidate.cpaXY(), + candidate.maxNormalisedDeltaIP(), + hfHelper.ctBplus(candidate), + candidate.eta(), + candidate.phi(), + hfHelper.yBplus(candidate), + hfHelper.eBplus(candidate), + flagMc, + invMassD0, + d0Cand.ptProng0(), + d0Cand.ptProng1(), + hfHelper.yD0(d0Cand), + d0Cand.eta(), + d0Cand.cpa(), + d0Cand.cpaXY(), + d0Cand.chi2PCA(), + d0Cand.decayLength(), + d0Cand.decayLengthXY(), + d0Cand.decayLengthNormalised(), + d0Cand.decayLengthXYNormalised(), + d0Cand.impactParameterProduct(), + d0Cand.impactParameter0(), + d0Cand.impactParameter1(), + d0Cand.impactParameterNormalised0(), + d0Cand.impactParameterNormalised1(), + d0Daughter0.tofNSigmaPi(), + d0Daughter0.tofNSigmaKa(), + d0Daughter0.tpcNSigmaPi(), + d0Daughter0.tpcNSigmaKa(), + d0Daughter1.tofNSigmaPi(), + d0Daughter1.tofNSigmaKa(), + d0Daughter1.tpcNSigmaPi(), + d0Daughter1.tpcNSigmaKa()); + } + } + void processData(aod::Collisions const& collisions, + soa::Filtered> const& candidates, + TracksWPid const&) + { // Filling event properties rowCandidateFullEvents.reserve(collisions.size()); for (const auto& collision : collisions) { - rowCandidateFullEvents( - collision.bcId(), - collision.numContrib(), - collision.posX(), - collision.posY(), - collision.posZ(), - 0, - 1); + fillEvent(collision, 0, 1); } // Filling candidate properties rowCandidateFull.reserve(candidates.size()); + if (fillCandidateLiteTable) { + rowCandidateLite.reserve(candidates.size()); + } for (const auto& candidate : candidates) { - auto fillTable = [&](int /*CandFlag*/, - // int FunctionSelection, - float FunctionInvMass, - float FunctionCt, - float FunctionY) { - auto d0Cand = candidate.prong0(); - auto piCand = candidate.prong1_as(); - // adding D0 daughters to the table - auto d0Daughter0 = d0Cand.prong0_as(); - auto d0Daughter1 = d0Cand.prong1_as(); - - auto invMassD0 = 0.; - if (piCand.sign() > 0) { - invMassD0 = hfHelper.invMassD0barToKPi(d0Cand); - } else if (piCand.sign() < 0) { - invMassD0 = hfHelper.invMassD0ToPiK(d0Cand); + if (fillOnlyBackground) { + float pseudoRndm = candidate.ptProng1() * 1000. - (int64_t)(candidate.ptProng1() * 1000); + if (candidate.pt() < ptMaxForDownSample && pseudoRndm >= downSampleBkgFactor) { + continue; } + } + auto prong1 = candidate.prong1_as(); + fillCandidateTable(candidate, prong1); + } + } - // if (FunctionSelection >= 1) { - if (std::abs(candidate.flagMcMatchRec()) >= isSignal) { + PROCESS_SWITCH(HfTreeCreatorBplusToD0Pi, processData, "Process data", true); - rowCandidateFull( - candidate.rSecondaryVertex(), - candidate.ptProng0(), - RecoDecay::p(candidate.pxProng0(), candidate.pyProng0(), candidate.pzProng0()), - candidate.ptProng1(), - RecoDecay::p(candidate.pxProng1(), candidate.pyProng1(), candidate.pzProng1()), - // 1 << CandFlag, - FunctionInvMass, - candidate.pt(), - candidate.p(), - FunctionCt, - candidate.eta(), - candidate.phi(), - FunctionY, - candidate.decayLength(), - candidate.decayLengthXY(), - candidate.decayLengthNormalised(), - candidate.decayLengthXYNormalised(), - candidate.cpa(), - candidate.cpaXY(), - candidate.impactParameterProduct(), - candidate.impactParameter0(), - candidate.impactParameter1(), - candidate.impactParameterNormalised0(), - candidate.impactParameterNormalised1(), - candidate.pxProng0(), - candidate.pyProng0(), - candidate.pzProng0(), - candidate.pxProng1(), - candidate.pyProng1(), - candidate.pzProng1(), - candidate.chi2PCA(), - piCand.tofNSigmaPi(), - piCand.tofNSigmaKa(), - piCand.tpcNSigmaPi(), - piCand.tpcNSigmaKa(), - candidate.flagMcMatchRec(), - invMassD0, - d0Cand.ptProng0(), - d0Cand.ptProng1(), - hfHelper.yD0(d0Cand), - d0Cand.eta(), - d0Cand.cpa(), - d0Cand.cpaXY(), - d0Cand.chi2PCA(), - d0Cand.decayLength(), - d0Cand.decayLengthXY(), - d0Cand.decayLengthNormalised(), - d0Cand.decayLengthXYNormalised(), - d0Cand.impactParameterProduct(), - d0Cand.impactParameter0(), - d0Cand.impactParameter1(), - d0Cand.impactParameterNormalised0(), - d0Cand.impactParameterNormalised1(), - d0Daughter0.tofNSigmaPi(), - d0Daughter0.tofNSigmaKa(), - d0Daughter0.tpcNSigmaPi(), - d0Daughter0.tpcNSigmaKa(), - d0Daughter1.tofNSigmaPi(), - d0Daughter1.tofNSigmaKa(), - d0Daughter1.tpcNSigmaPi(), - d0Daughter1.tpcNSigmaKa(), - candidate.globalIndex()); - } - }; + void processMc(aod::Collisions const& collisions, + aod::McCollisions const&, + SelectedCandidatesMc const& candidates, + soa::Join const& particles, + TracksWPid const&) + { + // Filling event properties + rowCandidateFullEvents.reserve(collisions.size()); + for (const auto& collision : collisions) { + fillEvent(collision, 0, 1); + } - // fillTable(0, candidate.isSelBplusToD0Pi(), hfHelper.invMassBplusToD0Pi(candidate), hfHelper.ctBplus(candidate), hfHelper.yBplus(candidate)); - fillTable(0, hfHelper.invMassBplusToD0Pi(candidate), hfHelper.ctBplus(candidate), hfHelper.yBplus(candidate)); + // Filling candidate properties + if (fillOnlySignal) { + rowCandidateFull.reserve(recSig.size()); + if (fillCandidateLiteTable) { + rowCandidateLite.reserve(recSig.size()); + } + for (const auto& candidate : recSig) { + auto prong1 = candidate.prong1_as(); + fillCandidateTable(candidate, prong1); + } + } else if (fillOnlyBackground) { + rowCandidateFull.reserve(recBg.size()); + if (fillCandidateLiteTable) { + rowCandidateLite.reserve(recBg.size()); + } + for (const auto& candidate : recBg) { + float pseudoRndm = candidate.ptProng1() * 1000. - (int64_t)(candidate.ptProng1() * 1000); + if (candidate.pt() < ptMaxForDownSample && pseudoRndm >= downSampleBkgFactor) { + continue; + } + auto prong1 = candidate.prong1_as(); + fillCandidateTable(candidate, prong1); + } + } else { + rowCandidateFull.reserve(candidates.size()); + if (fillCandidateLiteTable) { + rowCandidateLite.reserve(candidates.size()); + } + for (const auto& candidate : candidates) { + auto prong1 = candidate.prong1_as(); + fillCandidateTable(candidate, prong1); + } } // Filling particle properties rowCandidateFullParticles.reserve(particles.size()); for (const auto& particle : particles) { - if (std::abs(particle.flagMcMatchGen()) == 1 << aod::hf_cand_bplus::DecayType::BplusToD0Pi) { + if (std::abs(particle.flagMcMatchGen()) == 1 << aod::hf_cand_bplus::DecayTypeMc::BplusToD0PiToKPiPi) { rowCandidateFullParticles( particle.mcCollision().bcId(), particle.mcCollisionId(), particle.pt(), particle.eta(), particle.phi(), - RecoDecay::y(particle.pVector(), o2::constants::physics::MassBPlus), - particle.flagMcMatchGen(), - particle.globalIndex()); + RecoDecay::y(std::array{particle.px(), particle.py(), particle.pz()}, o2::constants::physics::MassBPlus), + particle.flagMcMatchGen()); } } } + + PROCESS_SWITCH(HfTreeCreatorBplusToD0Pi, processMc, "Process MC", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { - WorkflowSpec workflow; - workflow.push_back(adaptAnalysisTask(cfgc)); - return workflow; + return WorkflowSpec{adaptAnalysisTask(cfgc)}; }