From 4d404449adbd3914c32da1fd515ae21d49ad13ee Mon Sep 17 00:00:00 2001 From: Fabrizio Grosa Date: Tue, 27 Feb 2024 16:08:08 +0100 Subject: [PATCH 1/5] PWGHF: unify data creators for B hadrons --- PWGHF/D2H/DataModel/ReducedDataModel.h | 36 +- PWGHF/D2H/TableProducer/CMakeLists.txt | 9 +- .../candidateCreatorB0Reduced.cxx | 12 +- .../candidateCreatorBplusReduced.cxx | 4 +- .../dataCreatorCharmHadPiReduced.cxx | 929 ++++++++++++++++++ .../TableProducer/dataCreatorD0PiReduced.cxx | 568 ----------- .../dataCreatorDplusPiReduced.cxx | 645 ------------ 7 files changed, 963 insertions(+), 1240 deletions(-) create mode 100644 PWGHF/D2H/TableProducer/dataCreatorCharmHadPiReduced.cxx delete mode 100644 PWGHF/D2H/TableProducer/dataCreatorD0PiReduced.cxx delete mode 100644 PWGHF/D2H/TableProducer/dataCreatorDplusPiReduced.cxx diff --git a/PWGHF/D2H/DataModel/ReducedDataModel.h b/PWGHF/D2H/DataModel/ReducedDataModel.h index 9fea8ebf8c6..66ec55a161c 100644 --- a/PWGHF/D2H/DataModel/ReducedDataModel.h +++ b/PWGHF/D2H/DataModel/ReducedDataModel.h @@ -42,10 +42,13 @@ DECLARE_SOA_COLUMN(OriginalCollisionCount, originalCollisionCount, int); //! Siz } // namespace hf_reduced_collision DECLARE_SOA_TABLE(HfRedCollisions, "AOD", "HFREDCOLLISION", //! Table with collision for reduced workflow - soa::Index<>, + o2::soa::Index<>, collision::PosX, collision::PosY, collision::PosZ, + o2::soa::Marker<1>); + +DECLARE_SOA_TABLE(HfRedCollExtras, "AOD", "HFREDCOLLEXTRA", //! Table with collision extras for reduced workflow collision::CovXX, collision::CovXY, collision::CovYY, @@ -147,12 +150,15 @@ using HfRedTracks = HfRedTracksExt; namespace hf_charm_cand_reduced { -DECLARE_SOA_COLUMN(InvMass, invMass, float); //! Invariant mass of 2prong candidate in GeV/c2 -DECLARE_SOA_COLUMN(InvMassD0, invMassD0, float); //! Invariant mass of 2prong candidate in GeV/c2 -DECLARE_SOA_COLUMN(InvMassD0Bar, invMassD0Bar, float); //! Invariant mass of 2prong candidate in GeV/c2 -DECLARE_SOA_COLUMN(MlScoreBkg, mlScoreBkg, float); //! ML score for background class -DECLARE_SOA_COLUMN(MlScorePrompt, mlScorePrompt, float); //! ML score for prompt class -DECLARE_SOA_COLUMN(MlScoreNonprompt, mlScoreNonprompt, float); //! ML score for non-prompt class +DECLARE_SOA_COLUMN(InvMass, invMass, float); //! Invariant mass of 2prong candidate in GeV/c2 +DECLARE_SOA_COLUMN(InvMassD0, invMassD0, float); //! Invariant mass of 2prong candidate in GeV/c2 +DECLARE_SOA_COLUMN(InvMassD0Bar, invMassD0Bar, float); //! Invariant mass of 2prong candidate in GeV/c2 +DECLARE_SOA_COLUMN(MlScoreBkgMassHypo0, mlScoreBkgMassHypo0, float); //! ML score for background class (mass hypothesis 0) +DECLARE_SOA_COLUMN(MlScorePromptMassHypo0, mlScorePromptMassHypo0, float); //! ML score for prompt class (mass hypothesis 0) +DECLARE_SOA_COLUMN(MlScoreNonpromptMassHypo0, mlScoreNonpromptMassHypo0, float); //! ML score for non-prompt class (mass hypothesis 0) +DECLARE_SOA_COLUMN(MlScoreBkgMassHypo1, mlScoreBkgMassHypo1, float); //! ML score for background class (mass hypothesis 1) +DECLARE_SOA_COLUMN(MlScorePromptMassHypo1, mlScorePromptMassHypo1, float); //! ML score for prompt class (mass hypothesis 1) +DECLARE_SOA_COLUMN(MlScoreNonpromptMassHypo1, mlScoreNonpromptMassHypo1, float); //! ML score for non-prompt class (mass hypothesis 1) } // namespace hf_charm_cand_reduced // CAREFUL: need to follow convention [Name = Description + 's'] in DECLARE_SOA_TABLE(Name, "AOD", Description) @@ -174,9 +180,12 @@ DECLARE_SOA_TABLE(HfRed2ProngsCov, "AOD", "HFRED2PRONGSCOV", //! Table with 2pro o2::soa::Marker<1>); DECLARE_SOA_TABLE(HfRed2ProngsMl, "AOD", "HFRED2PRONGML", //! Table with 2prong candidate ML scores - hf_charm_cand_reduced::MlScoreBkg, - hf_charm_cand_reduced::MlScorePrompt, - hf_charm_cand_reduced::MlScoreNonprompt); + hf_charm_cand_reduced::MlScoreBkgMassHypo0, + hf_charm_cand_reduced::MlScorePromptMassHypo0, + hf_charm_cand_reduced::MlScoreNonpromptMassHypo0, + hf_charm_cand_reduced::MlScoreBkgMassHypo1, + hf_charm_cand_reduced::MlScorePromptMassHypo1, + hf_charm_cand_reduced::MlScoreNonpromptMassHypo1); // CAREFUL: need to follow convention [Name = Description + 's'] in DECLARE_SOA_TABLE(Name, "AOD", Description) // to call DECLARE_SOA_INDEX_COLUMN_FULL later on @@ -197,10 +206,9 @@ DECLARE_SOA_TABLE(HfRed3ProngsCov, "AOD", "HFRED3PRONGSCOV", //! Table with 3pro o2::soa::Marker<2>); DECLARE_SOA_TABLE(HfRed3ProngsMl, "AOD", "HFRED3PRONGML", //! Table with 3prong candidate ML scores - hf_charm_cand_reduced::MlScoreBkg, - hf_charm_cand_reduced::MlScorePrompt, - hf_charm_cand_reduced::MlScoreNonprompt, - o2::soa::Marker<1>); + hf_charm_cand_reduced::MlScoreBkgMassHypo0, + hf_charm_cand_reduced::MlScorePromptMassHypo0, + hf_charm_cand_reduced::MlScoreNonpromptMassHypo0); // Beauty candidates prongs namespace hf_cand_b0_reduced diff --git a/PWGHF/D2H/TableProducer/CMakeLists.txt b/PWGHF/D2H/TableProducer/CMakeLists.txt index 515c73b6b03..3d3e2103fd3 100644 --- a/PWGHF/D2H/TableProducer/CMakeLists.txt +++ b/PWGHF/D2H/TableProducer/CMakeLists.txt @@ -35,13 +35,8 @@ o2physics_add_dpl_workflow(candidate-selector-bplus-to-d0-pi-reduced # Data creators -o2physics_add_dpl_workflow(data-creator-d0-pi-reduced - SOURCES dataCreatorD0PiReduced.cxx - PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2::DCAFitter - COMPONENT_NAME Analysis) - -o2physics_add_dpl_workflow(data-creator-dplus-pi-reduced - SOURCES dataCreatorDplusPiReduced.cxx +o2physics_add_dpl_workflow(data-creator-charm-had-pi-reduced + SOURCES dataCreatorCharmHadPiReduced.cxx PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2::DCAFitter COMPONENT_NAME Analysis) diff --git a/PWGHF/D2H/TableProducer/candidateCreatorB0Reduced.cxx b/PWGHF/D2H/TableProducer/candidateCreatorB0Reduced.cxx index fca3ee11082..ce08d3bdf8b 100644 --- a/PWGHF/D2H/TableProducer/candidateCreatorB0Reduced.cxx +++ b/PWGHF/D2H/TableProducer/candidateCreatorB0Reduced.cxx @@ -58,6 +58,8 @@ struct HfCandidateCreatorB0Reduced { o2::vertexing::DCAFitterN<2> df2; // fitter for B vertex (2-prong vertex fitter) + 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; @@ -99,8 +101,8 @@ struct HfCandidateCreatorB0Reduced { /// \param tracksPionThisCollision pion tracks in this collision /// \param invMass2DPiMin minimum B0 invariant-mass /// \param invMass2DPiMax maximum B0 invariant-mass - template - void runCandidateCreation(aod::HfRedCollisions::iterator const& collision, + template + void runCandidateCreation(Coll const& collision, Cands const& candsDThisColl, Pions const& tracksPionThisCollision, const float& invMass2DPiMin, @@ -181,13 +183,13 @@ struct HfCandidateCreatorB0Reduced { rowCandidateProngs(candD.globalIndex(), trackPion.globalIndex()); if constexpr (withDmesMl) { - rowCandidateDmesMlScores(candD.mlScoreBkg(), candD.mlScorePrompt(), candD.mlScoreNonprompt()); + rowCandidateDmesMlScores(candD.mlScoreBkgMassHypo0(), candD.mlScorePromptMassHypo0(), candD.mlScoreNonpromptMassHypo0()); } } // pi loop } // D loop } - void processData(aod::HfRedCollisions const& collisions, + void processData(HfRedCollisionsWithExtras const& collisions, soa::Join const& candsD, soa::Join const& tracksPion, aod::HfOrigColCounts const& collisionsCounter, @@ -221,7 +223,7 @@ struct HfCandidateCreatorB0Reduced { PROCESS_SWITCH(HfCandidateCreatorB0Reduced, processData, "Process data without any ML score", true); - void processDataWithDmesMl(aod::HfRedCollisions const& collisions, + void processDataWithDmesMl(HfRedCollisionsWithExtras const& collisions, soa::Join const& candsD, soa::Join const& tracksPion, aod::HfOrigColCounts const& collisionsCounter, diff --git a/PWGHF/D2H/TableProducer/candidateCreatorBplusReduced.cxx b/PWGHF/D2H/TableProducer/candidateCreatorBplusReduced.cxx index d7ab973c2ec..a959e738b8a 100644 --- a/PWGHF/D2H/TableProducer/candidateCreatorBplusReduced.cxx +++ b/PWGHF/D2H/TableProducer/candidateCreatorBplusReduced.cxx @@ -56,6 +56,8 @@ struct HfCandidateCreatorBplusReduced { double bz{0.}; o2::vertexing::DCAFitterN<2> df2; // fitter for B vertex (2-prong vertex fitter) + using HfRedCollisionsWithExtras = soa::Join; + Preslice> candsDPerCollision = hf_track_index_reduced::hfRedCollisionId; Preslice> tracksPionPerCollision = hf_track_index_reduced::hfRedCollisionId; @@ -84,7 +86,7 @@ struct HfCandidateCreatorBplusReduced { registry.add("hEvents", "Events;;entries", HistType::kTH1F, {{1, 0.5, 1.5}}); } - void process(aod::HfRedCollisions const& collisions, + void process(HfRedCollisionsWithExtras const& collisions, soa::Join const& candsD, soa::Join const& tracksPion, aod::HfOrigColCounts const& collisionsCounter, diff --git a/PWGHF/D2H/TableProducer/dataCreatorCharmHadPiReduced.cxx b/PWGHF/D2H/TableProducer/dataCreatorCharmHadPiReduced.cxx new file mode 100644 index 00000000000..6b8d018e874 --- /dev/null +++ b/PWGHF/D2H/TableProducer/dataCreatorCharmHadPiReduced.cxx @@ -0,0 +1,929 @@ +// 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 dataCreatorCharmHadPiReduced.cxx +/// \brief Creation of CharmHadron-Pi pairs for Beauty hadron analyses +/// +/// \author Alexandre Bigot , IPHC Strasbourg +/// \author Antonio Palasciano , Università degli Studi di Bari +/// \author Fabrizio Grosa , CERN + +#include + +#include "CommonConstants/PhysicsConstants.h" +#include "DCAFitter/DCAFitterN.h" +#include "Framework/AnalysisTask.h" +#include "Framework/O2DatabasePDGPlugin.h" +#include "Framework/runDataProcessing.h" +#include "ReconstructionDataFormats/DCA.h" + +#include "Common/Core/trackUtilities.h" +#include "Common/DataModel/CollisionAssociationTables.h" + +#include "PWGHF/Core/HfHelper.h" +#include "PWGHF/DataModel/CandidateReconstructionTables.h" +#include "PWGHF/DataModel/CandidateSelectionTables.h" +#include "PWGHF/Utils/utilsBfieldCCDB.h" +#include "PWGHF/D2H/DataModel/ReducedDataModel.h" + +using namespace o2; +using namespace o2::analysis; +using namespace o2::aod; +using namespace o2::constants::physics; +using namespace o2::framework; +using namespace o2::framework::expressions; + +enum Event : uint8_t { + Processed = 0, + NoCharmHadPiSelected, + CharmHadPiSelected, + kNEvent +}; + +enum DecayChannel : uint8_t { + B0ToDminusPi = 0, + BplusToD0barPi, +}; + +/// Creation of CharmHad-Pi pairs for Beauty hadrons +struct HfDataCreatorCharmHadPiReduced { + // Produces AOD tables to store track information + // collision related tables + Produces hfReducedCollision; + Produces hfReducedCollExtra; + Produces hfCollisionCounter; + // Pi bachelor related tables + Produces hfTrackPion; + Produces hfTrackCovPion; + Produces hfTrackPidPion; + // charm hadron related tables + Produces hfCand2Prong; + Produces hfCand2ProngCov; + Produces hfCand2ProngMl; + Produces hfCand3Prong; + Produces hfCand3ProngCov; + Produces hfCand3ProngMl; + // B-hadron config and MC related tables + Produces rowCandidateConfigB0; + Produces rowHfDPiMcRecReduced; + Produces rowHfDPiMcCheckReduced; + Produces rowHfB0McGenReduced; + + Produces rowCandidateConfigBplus; + Produces rowHfD0PiMcRecReduced; + Produces rowHfBpMcGenReduced; + + // vertexing + // Configurable bz{"bz", 5., "magnetic field"}; + Configurable propagateToPCA{"propagateToPCA", true, "create tracks version propagated to PCA"}; + Configurable useAbsDCA{"useAbsDCA", false, "Minimise abs. distance rather than chi2"}; + Configurable useWeightedFinalPCA{"useWeightedFinalPCA", false, "Recalculate vertex position using track covariances, effective only if useAbsDCA is true"}; + Configurable maxR{"maxR", 200., "reject PCA's above this radius"}; + Configurable maxDZIni{"maxDZIni", 4., "reject (if>0) PCA candidate if tracks DZ exceeds threshold"}; + Configurable minParamChange{"minParamChange", 1.e-3, "stop iterations if largest change of any B0 is smaller than this"}; + Configurable minRelChi2Change{"minRelChi2Change", 0.9, "stop iterations is chi2/chi2old > this"}; + // selection + Configurable usePionIsGlobalTrackWoDCA{"usePionIsGlobalTrackWoDCA", true, "check isGlobalTrackWoDCA status for pions, for Run3 studies"}; + Configurable ptPionMin{"ptPionMin", 0.5, "minimum pion pT threshold (GeV/c)"}; + Configurable> binsPtPion{"binsPtPion", std::vector{hf_cuts_single_track::vecBinsPtTrack}, "track pT bin limits for pion DCA XY pT-dependent cut"}; + Configurable> cutsTrackPionDCA{"cutsTrackPionDCA", {hf_cuts_single_track::cutsTrack[0], hf_cuts_single_track::nBinsPtTrack, hf_cuts_single_track::nCutVarsTrack, hf_cuts_single_track::labelsPtTrack, hf_cuts_single_track::labelsCutVarTrack}, "Single-track selections per pT bin for pions"}; + Configurable invMassWindowCharmHadPi{"invMassWindowCharmHadPi", 0.3, "invariant-mass window for CharmHad-Pi pair preselections (GeV/c2)"}; + Configurable selectionFlagDplus{"selectionFlagDplus", 7, "Selection Flag for D+"}; + Configurable selectionFlagD0{"selectionFlagD0", 1, "Selection Flag for D0"}; + Configurable selectionFlagD0bar{"selectionFlagD0bar", 1, "Selection Flag for D0bar"}; + + // magnetic field setting from CCDB + Configurable ccdbUrl{"ccdbUrl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; + Configurable ccdbPathGrpMag{"ccdbPathGrpMag", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object (Run 3)"}; + + // MC extra + Configurable checkDecayTypeMc{"checkDecayTypeMc", false, "flag to enable MC checks on decay type"}; + + HfHelper hfHelper; + + // CCDB service + Service ccdb; + o2::base::Propagator::MatCorrType noMatCorr = o2::base::Propagator::MatCorrType::USEMatCorrNONE; + int runNumber; + + // O2DatabasePDG service + Service pdg; + + double massPi{0.}; + double massC{0.}; + double massB{0.}; + double invMass2ChHadPiMin{0.}; + double invMass2ChHadPiMax{0.}; + double bz{0.}; + + bool isHfCandBhadConfigFilled = false; + + // Fitter to redo D-vertex to get extrapolated daughter tracks (2/3-prong vertex filter) + o2::vertexing::DCAFitterN<3> df3; + o2::vertexing::DCAFitterN<2> df2; + + using TracksPidPi = soa::Join; + using TracksPidWithSel = soa::Join; + using TracksPidWithSelAndMc = soa::Join; + + using CandsDplusFiltered = soa::Filtered>; + using CandsDplusFilteredWithMl = soa::Filtered>; + using CandsD0Filtered = soa::Filtered>; + using CandsD0FilteredWithMl = soa::Filtered>; + + Filter filterSelectDplusCandidates = (aod::hf_sel_candidate_dplus::isSelDplusToPiKPi >= selectionFlagDplus); + Filter filterSelectDzeroCandidates = (aod::hf_sel_candidate_d0::isSelD0 >= selectionFlagD0 || aod::hf_sel_candidate_d0::isSelD0bar >= selectionFlagD0bar); + + Preslice candsDplusPerCollision = aod::track_association::collisionId; + Preslice candsDplusPerCollisionWithMl = aod::track_association::collisionId; + Preslice candsD0PerCollision = aod::track_association::collisionId; + Preslice candsD0PerCollisionWithMl = aod::track_association::collisionId; + Preslice trackIndicesPerCollision = aod::track_association::collisionId; + + HistogramRegistry registry{"registry"}; + + void init(InitContext const&) + { + // histograms + constexpr int kNBinsEvents = kNEvent; + std::string labels[kNBinsEvents]; + labels[Event::Processed] = "processed"; + labels[Event::NoCharmHadPiSelected] = "without CharmHad-Pi pairs"; + labels[Event::CharmHadPiSelected] = "with CharmHad-Pi pairs"; + static const AxisSpec axisEvents = {kNBinsEvents, 0.5, kNBinsEvents + 0.5, ""}; + registry.add("hEvents", "Events;;entries", HistType::kTH1F, {axisEvents}); + for (int iBin = 0; iBin < kNBinsEvents; iBin++) { + registry.get(HIST("hEvents"))->GetXaxis()->SetBinLabel(iBin + 1, labels[iBin].data()); + } + + registry.add("hMassCharmHad0", "D^{#minus} candidates;inv. mass (p^{#minus} K^{#plus} #pi^{#minus}) (GeV/#it{c}^{2});entries", {HistType::kTH1F, {{500, 0., 5.}}}); + registry.add("hMassCharmHad1", "D^{#minus} candidates;inv. mass (p^{#minus} K^{#plus} #pi^{#minus}) (GeV/#it{c}^{2});entries", {HistType::kTH1F, {{500, 0., 5.}}}); + registry.add("hPtCharmHad", "D^{#minus} candidates;D^{#minus} candidate #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{100, 0., 10.}}}); + registry.add("hPtPion", "#pi^{#plus} candidates;#pi^{#plus} candidate #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{100, 0., 10.}}}); + registry.add("hCpaCharmHad", "D^{#minus} candidates;D^{#minus} cosine of pointing angle;entries", {HistType::kTH1F, {{110, -1.1, 1.1}}}); + + // Initialize fitter + df2.setPropagateToPCA(propagateToPCA); + df2.setMaxR(maxR); + df2.setMaxDZIni(maxDZIni); + df2.setMinParamChange(minParamChange); + df2.setMinRelChi2Change(minRelChi2Change); + df2.setUseAbsDCA(useAbsDCA); + df2.setWeightedFinalPCA(useWeightedFinalPCA); + df2.setMatCorrType(noMatCorr); + + df3.setPropagateToPCA(propagateToPCA); + df3.setMaxR(maxR); + df3.setMaxDZIni(maxDZIni); + df3.setMinParamChange(minParamChange); + df3.setMinRelChi2Change(minRelChi2Change); + df3.setUseAbsDCA(useAbsDCA); + df3.setWeightedFinalPCA(useWeightedFinalPCA); + df3.setMatCorrType(noMatCorr); + + // Configure CCDB access + ccdb->setURL(ccdbUrl); + ccdb->setCaching(true); + ccdb->setLocalObjectValidityChecking(); + runNumber = 0; + + // invariant-mass window cut + massPi = MassPiPlus; + if (doprocessDplusPiData || doprocessDplusPiDataWithMl || doprocessDplusPiMc || doprocessDplusPiMcWithMl) { + massC = MassDMinus; + massB = MassB0; + } else if (doprocessD0PiData || doprocessD0PiDataWithMl || doprocessD0PiMc || doprocessD0PiMcWithMl) { + massC = MassD0; + massB = MassBPlus; + } + invMass2ChHadPiMin = (massB - invMassWindowCharmHadPi) * (massB - invMassWindowCharmHadPi); + invMass2ChHadPiMax = (massB + invMassWindowCharmHadPi) * (massB + invMassWindowCharmHadPi); + } + + /// Pion selection (D Pi <-- B0) + /// \param trackPion is a track with the pion hypothesis + /// \param trackParCovPion is the track parametrisation of the pion + /// \param dcaPion is the 2-D array with track DCAs of the pion + /// \param charmDautracks charm-hadron daughter tracks + /// \return true if trackPion passes all cuts + template + bool isPionSelected(const T1& trackPion, const T2& trackParCovPion, const T3& dcaPion, const std::vector& charmDautracks) + { + // check isGlobalTrackWoDCA status for pions if wanted + if (usePionIsGlobalTrackWoDCA && !trackPion.isGlobalTrackWoDCA()) { + return false; + } + // minimum pT selection + if (trackParCovPion.getPt() < ptPionMin || !isSelectedTrackDCA(trackParCovPion, dcaPion)) { + return false; + } + // reject pions that are charm-hadron daughters + for (const auto& track : charmDautracks) { + if (trackPion.globalIndex() == track.globalIndex()) { + return false; + } + } + + return true; + } + + /// Single-track cuts for pions on dcaXY + /// \param trackPar is the track parametrisation + /// \param dca is the 2-D array with track DCAs + /// \return true if track passes all cuts + template + bool isSelectedTrackDCA(const T1& trackPar, const T2& dca) + { + auto pTBinTrack = findBin(binsPtPion, trackPar.getPt()); + if (pTBinTrack == -1) { + return false; + } + + if (std::abs(dca[0]) < cutsTrackPionDCA->get(pTBinTrack, "min_dcaxytoprimary")) { + return false; // minimum DCAxy + } + if (std::abs(dca[0]) > cutsTrackPionDCA->get(pTBinTrack, "max_dcaxytoprimary")) { + return false; // maximum DCAxy + } + return true; + } + + /// Function for filling MC reco information in the tables + /// \param particlesMc is the table with MC particles + /// \param vecDaughtersB is the vector with all daughter tracks (bachelor pion in last position) + /// \param indexHfCandCharm is the index of the charm-hadron candidate + /// \param selectedTracksPion is the map with the indices of selected bachelor pion tracks + template + void fillMcRecoInfo(const PParticles& particlesMc, + const std::vector& vecDaughtersB, + int& indexHfCandCharm, + std::map selectedTracksPion) { + + // we check the MC matching to be stored + int8_t sign{0}; + int8_t flag{0}; + int8_t debug{0}; + int pdgCodeBeautyMother{-1}; + int pdgCodeCharmMother{-1}; + int pdgCodeProng0{0}; + int pdgCodeProng1{0}; + int pdgCodeProng2{0}; + int pdgCodeProng3{0}; + float motherPt{-1.f}; + + if constexpr (decChannel == DecayChannel::B0ToDminusPi) { + // B0 → D- π+ → (π- K+ π-) π+ + auto indexRec = RecoDecay::getMatchedMCRec(particlesMc, std::array{vecDaughtersB[0], vecDaughtersB[1], vecDaughtersB[2], vecDaughtersB[3]}, Pdg::kB0, std::array{-kPiPlus, +kKPlus, -kPiPlus, +kPiPlus}, true, &sign, 3); + if (indexRec > -1) { + // D- → π- K+ π- + // Printf("Checking D- → π- K+ π-"); + indexRec = RecoDecay::getMatchedMCRec(particlesMc, std::array{vecDaughtersB[0], vecDaughtersB[1], vecDaughtersB[2]}, Pdg::kDMinus, std::array{-kPiPlus, +kKPlus, -kPiPlus}, true, &sign, 2); + if (indexRec > -1) { + flag = sign * BIT(hf_cand_b0::DecayTypeMc::B0ToDplusPiToPiKPiPi); + } else { + debug = 1; + LOGF(debug, "B0 decays in the expected final state but the condition on the intermediate state is not fulfilled"); + } + + auto indexMother = RecoDecay::getMother(particlesMc, vecDaughtersB.back().template mcParticle_as(), Pdg::kB0, true); + if (indexMother >= 0) { + auto particleMother = particlesMc.rawIteratorAt(indexMother); + motherPt = particleMother.pt(); + } + } + rowHfDPiMcRecReduced(indexHfCandCharm, selectedTracksPion[vecDaughtersB.back().globalIndex()], flag, debug, motherPt); + + // additional checks for correlated backgrounds + if (checkDecayTypeMc) { + // B0 → Ds- π+ → (K- K+ π-) π+ + if (!flag) { + indexRec = RecoDecay::getMatchedMCRec(particlesMc, std::array{vecDaughtersB[0], vecDaughtersB[1], vecDaughtersB[2], vecDaughtersB[3]}, Pdg::kB0, std::array{-kKPlus, +kKPlus, -kPiPlus, +kPiPlus}, true, &sign, 3); + if (indexRec > -1) { + // Ds- → K- K+ π- + indexRec = RecoDecay::getMatchedMCRec(particlesMc, std::array{vecDaughtersB[0], vecDaughtersB[1], vecDaughtersB[2]}, -Pdg::kDS, std::array{-kKPlus, +kKPlus, -kPiPlus}, true, &sign, 2); + if (indexRec > -1) { + flag = sign * BIT(hf_cand_b0::DecayTypeMc::B0ToDsPiToKKPiPi); + } + } + } + // Partly reconstructed decays, i.e. the 4 prongs have a common b-hadron ancestor + // convention: final state particles are prong0,1,2,3 + if (!flag) { + auto particleProng0 = vecDaughtersB[0].mcParticle(); + auto particleProng1 = vecDaughtersB[1].mcParticle(); + auto particleProng2 = vecDaughtersB[2].mcParticle(); + auto particleProng3 = vecDaughtersB[3].mcParticle(); + // b-hadron hypothesis + std::array bHadronMotherHypos = {Pdg::kB0, Pdg::kBS, Pdg::kLambdaB0}; + // c-hadron hypothesis + std::array cHadronMotherHypos = {Pdg::kDPlus, Pdg::kDS, Pdg::kDStar}; + + 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); + int index3Mother = RecoDecay::getMother(particlesMc, particleProng3, bHadronMotherHypo, true); + + // look for common b-hadron ancestor + if (index0Mother > -1 && index1Mother > -1 && index2Mother > -1 && index3Mother > -1) { + if (index0Mother == index1Mother && index1Mother == index2Mother && index2Mother == index3Mother) { + flag = BIT(hf_cand_b0::DecayTypeMc::PartlyRecoDecay); + pdgCodeBeautyMother = particlesMc.rawIteratorAt(index0Mother).pdgCode(); + pdgCodeCharmMother = 0; + pdgCodeProng0 = particleProng0.pdgCode(); + pdgCodeProng1 = particleProng1.pdgCode(); + pdgCodeProng2 = particleProng2.pdgCode(); + pdgCodeProng3 = particleProng3.pdgCode(); + // look for common c-hadron mother among prongs 0, 1 and 2 + for (const auto& cHadronMotherHypo : cHadronMotherHypos) { + int8_t depthMax = 2; + if (cHadronMotherHypo == Pdg::kDStar) { // to include D* -> D π0/γ and D* -> D0 π + depthMax += 1; + } + int index0CharmMother = RecoDecay::getMother(particlesMc, particleProng0, cHadronMotherHypo, true, &sign, depthMax); + int index1CharmMother = RecoDecay::getMother(particlesMc, particleProng1, cHadronMotherHypo, true, &sign, depthMax); + int index2CharmMother = RecoDecay::getMother(particlesMc, particleProng2, cHadronMotherHypo, true, &sign, depthMax); + if (index0CharmMother > -1 && index1CharmMother > -1 && index2CharmMother > -1) { + if (index0CharmMother == index1CharmMother && index1CharmMother == index2CharmMother) { + // pdgCodeCharmMother = + // Pdg::kDPlus (if D+ is the mother and does not come from D*+) + // Pdg::kDPlus + Pdg::kDStar (if D+ is the mother and D*+ -> D+ π0/γ) + // Pdg::kDStar (if D*+ is the mother and D*+ -> D0 π+) + // Pdg::kDS (if Ds is the mother) + pdgCodeCharmMother += std::abs(particlesMc.rawIteratorAt(index0CharmMother).pdgCode()); + } + } + } + break; + } + } + } + } + rowHfDPiMcCheckReduced(pdgCodeBeautyMother, pdgCodeCharmMother, pdgCodeProng0, pdgCodeProng1, pdgCodeProng2, pdgCodeProng3); + } + } + else if constexpr (decChannel == DecayChannel::BplusToD0barPi) { + // 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"); + } + + auto indexMother = RecoDecay::getMother(particlesMc, vecDaughtersB.back().template mcParticle_as(), Pdg::kBPlus, true); + if (indexMother >= 0) { + auto particleMother = particlesMc.rawIteratorAt(indexMother); + motherPt = particleMother.pt(); + } + } + rowHfD0PiMcRecReduced(indexHfCandCharm, selectedTracksPion[vecDaughtersB.back().globalIndex()], flag, debug, motherPt); + } + } + + template + void runDataCreation(aod::Collision const& collision, + CCharmCands const& candsC, + aod::TrackAssoc const& trackIndices, + TTracks const& tracks, + PParticles const& particlesMc, + aod::BCsWithTimestamps const& bcs) + { + // helpers for ReducedTables filling + int indexHfReducedCollision = hfReducedCollision.lastIndex() + 1; + // std::map where the key is the track.globalIndex() and + // the value is the track index in the table of the selected pions + std::map selectedTracksPion; + bool fillHfReducedCollision = false; + + auto primaryVertex = getPrimaryVertex(collision); + + // Set the magnetic field from ccdb. + // The static instance of the propagator was already modified in the HFTrackIndexSkimCreator, + // but this is not true when running on Run2 data/MC already converted into AO2Ds. + auto bc = collision.bc_as(); + if (runNumber != bc.runNumber()) { + LOG(info) << ">>>>>>>>>>>> Current run number: " << runNumber; + o2::parameters::GRPMagField* grpo = ccdb->getForTimeStamp(ccdbPathGrpMag, bc.timestamp()); + if (grpo == nullptr) { + LOGF(fatal, "Run 3 GRP object (type o2::parameters::GRPMagField) is not available in CCDB for run=%d at timestamp=%llu", bc.runNumber(), bc.timestamp()); + } + o2::base::Propagator::initFieldFromGRP(grpo); + bz = o2::base::Propagator::Instance()->getNominalBz(); + LOG(info) << ">>>>>>>>>>>> Magnetic field: " << bz; + runNumber = bc.runNumber(); + } + df3.setBz(bz); + + auto thisCollId = collision.globalIndex(); + for (const auto& candC : candsC) { + int indexHfCandCharm{-1}; + float invMassC0{-1.f}, invMassC1{-1.f}; + if constexpr (decChannel == DecayChannel::B0ToDminusPi) { + indexHfCandCharm = hfCand3Prong.lastIndex() + 1; + invMassC0 = hfHelper.invMassDplusToPiKPi(candC); + registry.fill(HIST("hMassCharmHad"), invMassC0); + } else if constexpr (decChannel == DecayChannel::BplusToD0barPi){ + indexHfCandCharm = hfCand2Prong.lastIndex() + 1; + if (candC.isSelD0() >= selectionFlagD0) { + invMassC0 = hfHelper.invMassD0ToPiK(candC); + registry.fill(HIST("hMassCharmHad0"), invMassC0); + } + if (candC.isSelD0bar() >= selectionFlagD0bar) { + invMassC1 = hfHelper.invMassD0barToKPi(candC); + registry.fill(HIST("hMassCharmHad1"), invMassC1); + } + } + registry.fill(HIST("hPtCharmHad"), candC.pt()); + registry.fill(HIST("hCpaCharmHad"), candC.cpa()); + + bool fillHfCandCharm = false; + + std::vector charmHadDauTracks{candC.template prong0_as(), candC.template prong1_as()}; + o2::track::TrackParCov trackParCov0 = getTrackParCov(charmHadDauTracks[0]); + o2::track::TrackParCov trackParCov1 = getTrackParCov(charmHadDauTracks[1]); + o2::track::TrackParCov trackParCov2{}; + + std::array pVec0 = {charmHadDauTracks[0].px(), charmHadDauTracks[0].py(), charmHadDauTracks[0].pz()}; + std::array pVec1 = {charmHadDauTracks[1].px(), charmHadDauTracks[1].py(), charmHadDauTracks[1].pz()}; + std::array pVec2{}; + + auto dca0 = o2::dataformats::DCA(charmHadDauTracks[0].dcaXY(), charmHadDauTracks[0].dcaZ(), charmHadDauTracks[0].cYY(), charmHadDauTracks[0].cZY(), charmHadDauTracks[0].cZZ()); + auto dca1 = o2::dataformats::DCA(charmHadDauTracks[1].dcaXY(), charmHadDauTracks[1].dcaZ(), charmHadDauTracks[1].cYY(), charmHadDauTracks[1].cZY(), charmHadDauTracks[1].cZZ()); + + // repropagate tracks to this collision if needed + if (charmHadDauTracks[0].collisionId() != thisCollId) { + trackParCov0.propagateToDCA(primaryVertex, bz, &dca0); + } + + if (charmHadDauTracks[1].collisionId() != thisCollId) { + trackParCov1.propagateToDCA(primaryVertex, bz, &dca1); + } + + // third track, if it's a 3-prong + if constexpr (decChannel == DecayChannel::B0ToDminusPi) { + charmHadDauTracks.push_back(candC.template prong2_as()); + trackParCov2 = getTrackParCov(charmHadDauTracks[2]); + pVec2 = std::array{charmHadDauTracks[2].px(), charmHadDauTracks[2].py(), charmHadDauTracks[2].pz()}; + auto dca2 = o2::dataformats::DCA(charmHadDauTracks[2].dcaXY(), charmHadDauTracks[2].dcaZ(), charmHadDauTracks[2].cYY(), charmHadDauTracks[2].cZY(), charmHadDauTracks[2].cZZ()); + + if (charmHadDauTracks[2].collisionId() != thisCollId) { + trackParCov2.propagateToDCA(primaryVertex, bz, &dca2); + } + } + + // --------------------------------- + // reconstruct charm candidate secondary vertex + o2::track::TrackParCov trackParCovCharmHad{}; + std::array pVecCharm{}; + if constexpr (decChannel == DecayChannel::B0ToDminusPi) { // D∓ → π∓ K± π∓ + if (df3.process(trackParCov0, trackParCov1, trackParCov2) == 0) { + continue; + } + auto secondaryVertexCharm = df3.getPCACandidate(); + trackParCov0.propagateTo(secondaryVertexCharm[0], bz); + trackParCov1.propagateTo(secondaryVertexCharm[0], bz); + trackParCov2.propagateTo(secondaryVertexCharm[0], bz); + df3.getTrack(0).getPxPyPzGlo(pVec0); + df3.getTrack(1).getPxPyPzGlo(pVec1); + df3.getTrack(2).getPxPyPzGlo(pVec2); + pVecCharm = RecoDecay::pVec(pVec0, pVec1, pVec2); + trackParCovCharmHad.setAbsCharge(charmHadDauTracks[1].sign()); // to be sure + } else if constexpr (decChannel == DecayChannel::BplusToD0barPi) { // D0(bar) → K± π∓ + if (df2.process(trackParCov0, trackParCov1) == 0) { + continue; + } + auto secondaryVertexCharm = df2.getPCACandidate(); + trackParCov0.propagateTo(secondaryVertexCharm[0], bz); + trackParCov1.propagateTo(secondaryVertexCharm[0], bz); + df2.getTrack(0).getPxPyPzGlo(pVec0); + df2.getTrack(1).getPxPyPzGlo(pVec1); + pVecCharm = RecoDecay::pVec(pVec0, pVec1); + trackParCovCharmHad = df2.createParentTrackParCov(); + trackParCovCharmHad.setAbsCharge(0); // to be sure + } + + for (const auto& trackId : trackIndices) { + auto trackPion = trackId.template track_as(); + + // apply selections on pion tracks + auto trackParCovPion = getTrackParCov(trackPion); + o2::gpu::gpustd::array dcaPion{trackPion.dcaXY(), trackPion.dcaZ()}; + std::array pVecPion = {trackPion.px(), trackPion.py(), trackPion.pz()}; + if (trackPion.collisionId() != thisCollId) { + o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, trackParCovPion, 2.f, noMatCorr, &dcaPion); + getPxPyPz(trackParCovPion, pVecPion); + } + + // reject pi D with same sign as D + if constexpr (decChannel == DecayChannel::B0ToDminusPi) { // D∓ → π∓ K± π∓ + if (trackPion.sign() * charmHadDauTracks[0].sign() > 0) { + continue; + } + } else if constexpr (decChannel == DecayChannel::BplusToD0barPi) { // D0(bar) → K± π∓ + if (!((candC.isSelD0() >= selectionFlagD0 && trackPion.sign() < 0) || (candC.isSelD0bar() >= selectionFlagD0bar && trackPion.sign() > 0))) { + continue; + } + } + + // apply selections on pion tracks + if (!isPionSelected(trackPion, trackParCovPion, dcaPion, charmHadDauTracks)) { + continue; + } + + registry.fill(HIST("hPtPion"), trackParCovPion.getPt()); + // compute invariant mass square and apply selection + auto invMass2DPi = RecoDecay::m2(std::array{pVecCharm, pVecPion}, std::array{massC, massPi}); + if ((invMass2DPi < invMass2ChHadPiMin) || (invMass2DPi > invMass2ChHadPiMax)) { + continue; + } + + // fill Pion tracks table + // if information on track already stored, go to next track + if (!selectedTracksPion.count(trackPion.globalIndex())) { + hfTrackPion(trackPion.globalIndex(), indexHfReducedCollision, + trackParCovPion.getX(), trackParCovPion.getAlpha(), + trackParCovPion.getY(), trackParCovPion.getZ(), trackParCovPion.getSnp(), + trackParCovPion.getTgl(), trackParCovPion.getQ2Pt()); + hfTrackCovPion(trackParCovPion.getSigmaY2(), trackParCovPion.getSigmaZY(), trackParCovPion.getSigmaZ2(), + trackParCovPion.getSigmaSnpY(), trackParCovPion.getSigmaSnpZ(), + trackParCovPion.getSigmaSnp2(), trackParCovPion.getSigmaTglY(), trackParCovPion.getSigmaTglZ(), + trackParCovPion.getSigmaTglSnp(), trackParCovPion.getSigmaTgl2(), + trackParCovPion.getSigma1PtY(), trackParCovPion.getSigma1PtZ(), trackParCovPion.getSigma1PtSnp(), + trackParCovPion.getSigma1PtTgl(), trackParCovPion.getSigma1Pt2()); + hfTrackPidPion(trackPion.hasTPC(), trackPion.hasTOF(), + trackPion.tpcNSigmaPi(), trackPion.tofNSigmaPi()); + // add trackPion.globalIndex() to a list + // to keep memory of the pions filled in the table and avoid refilling them if they are paired to another D candidate + // and keep track of their index in hfTrackPion for McRec purposes + selectedTracksPion[trackPion.globalIndex()] = hfTrackPion.lastIndex(); + } + + if constexpr (doMc) { + std::vector beautyHadDauTracks{}; + std::copy(charmHadDauTracks.begin(), charmHadDauTracks.end(), beautyHadDauTracks.begin()); + beautyHadDauTracks.push_back(trackPion); + fillMcRecoInfo(particlesMc, beautyHadDauTracks, indexHfCandCharm, selectedTracksPion); + } + fillHfCandCharm = true; + } // pion loop + if (fillHfCandCharm) { // fill candCplus table only once per D candidate + if constexpr (decChannel == DecayChannel::B0ToDminusPi) { // D∓ → π∓ K± π∓ + hfCand3Prong(charmHadDauTracks[0].globalIndex(), charmHadDauTracks[1].globalIndex(), charmHadDauTracks[2].globalIndex(), + indexHfReducedCollision, + trackParCovCharmHad.getX(), trackParCovCharmHad.getAlpha(), + trackParCovCharmHad.getY(), trackParCovCharmHad.getZ(), trackParCovCharmHad.getSnp(), + trackParCovCharmHad.getTgl(), trackParCovCharmHad.getQ2Pt(), + candC.xSecondaryVertex(), candC.ySecondaryVertex(), candC.zSecondaryVertex(), invMassC0); + hfCand3ProngCov(trackParCovCharmHad.getSigmaY2(), trackParCovCharmHad.getSigmaZY(), trackParCovCharmHad.getSigmaZ2(), + trackParCovCharmHad.getSigmaSnpY(), trackParCovCharmHad.getSigmaSnpZ(), + trackParCovCharmHad.getSigmaSnp2(), trackParCovCharmHad.getSigmaTglY(), trackParCovCharmHad.getSigmaTglZ(), + trackParCovCharmHad.getSigmaTglSnp(), trackParCovCharmHad.getSigmaTgl2(), + trackParCovCharmHad.getSigma1PtY(), trackParCovCharmHad.getSigma1PtZ(), trackParCovCharmHad.getSigma1PtSnp(), + trackParCovCharmHad.getSigma1PtTgl(), trackParCovCharmHad.getSigma1Pt2()); + if constexpr (withMl) { + hfCand3ProngMl(candC.mlProbDplusToPiKPi()[0], candC.mlProbDplusToPiKPi()[1], candC.mlProbDplusToPiKPi()[2]); + } + } else if constexpr (decChannel == DecayChannel::BplusToD0barPi) { // D0(bar) → K± π∓ + hfCand2Prong(charmHadDauTracks[0].globalIndex(), charmHadDauTracks[1].globalIndex(), + indexHfReducedCollision, + trackParCovCharmHad.getX(), trackParCovCharmHad.getAlpha(), + trackParCovCharmHad.getY(), trackParCovCharmHad.getZ(), trackParCovCharmHad.getSnp(), + trackParCovCharmHad.getTgl(), trackParCovCharmHad.getQ2Pt(), + candC.xSecondaryVertex(), candC.ySecondaryVertex(), candC.zSecondaryVertex(), invMassC0, invMassC1); + hfCand2ProngCov(trackParCovCharmHad.getSigmaY2(), trackParCovCharmHad.getSigmaZY(), trackParCovCharmHad.getSigmaZ2(), + trackParCovCharmHad.getSigmaSnpY(), trackParCovCharmHad.getSigmaSnpZ(), + trackParCovCharmHad.getSigmaSnp2(), trackParCovCharmHad.getSigmaTglY(), trackParCovCharmHad.getSigmaTglZ(), + trackParCovCharmHad.getSigmaTglSnp(), trackParCovCharmHad.getSigmaTgl2(), + trackParCovCharmHad.getSigma1PtY(), trackParCovCharmHad.getSigma1PtZ(), trackParCovCharmHad.getSigma1PtSnp(), + trackParCovCharmHad.getSigma1PtTgl(), trackParCovCharmHad.getSigma1Pt2()); + if constexpr (withMl) { + std::array bdtScores = {-1.f, -1.f, -1.f, -1.f, -1.f, -1.f}; + if (candC.mlProbD0().size() == 3) { + std::copy(candC.mlProbD0().begin(), candC.mlProbD0().end(), bdtScores.begin()); + } + if (candC.mlProbD0bar().size() == 3) { + std::copy(candC.mlProbD0bar().begin(), candC.mlProbD0bar().end(), bdtScores.begin() + 3); + } + + hfCand2ProngMl(bdtScores[0], bdtScores[1], bdtScores[2], bdtScores[3], bdtScores[4], bdtScores[5]); + } + } + fillHfReducedCollision = true; + } + } // candsC loop + + registry.fill(HIST("hEvents"), 1 + Event::Processed); + if (!fillHfReducedCollision) { + registry.fill(HIST("hEvents"), 1 + Event::NoCharmHadPiSelected); + return; + } + registry.fill(HIST("hEvents"), 1 + Event::CharmHadPiSelected); + // fill collision table if it contains a DPi pair a minima + hfReducedCollision(collision.posX(), collision.posY(), collision.posZ()); + hfReducedCollExtra(collision.covXX(), collision.covXY(), collision.covYY(), + collision.covXZ(), collision.covYZ(), collision.covZZ(), + bz); + } + + template + void runMcGen(aod::McParticles const& particlesMc) + { + // Match generated particles. + for (const auto& particle : particlesMc) { + int8_t sign{0}; + int8_t flag{0}; + if constexpr (decayChannel == DecayChannel::B0ToDminusPi) { + // B0 → D- π+ + if (RecoDecay::isMatchedMCGen(particlesMc, particle, Pdg::kB0, std::array{-static_cast(Pdg::kDPlus), +kPiPlus}, true)) { + // Match D- -> π- K+ π- + auto candCMC = particlesMc.rawIteratorAt(particle.daughtersIds().front()); + // Printf("Checking D- -> π- K+ π-"); + if (RecoDecay::isMatchedMCGen(particlesMc, candCMC, -static_cast(Pdg::kDPlus), std::array{-kPiPlus, +kKPlus, -kPiPlus}, true, &sign)) { + flag = sign * BIT(hf_cand_b0::DecayType::B0ToDPi); + } + } + + // save information for B0 task + if (!TESTBIT(std::abs(flag), hf_cand_b0::DecayType::B0ToDPi)) { + continue; + } + + auto ptParticle = particle.pt(); + auto yParticle = RecoDecay::y(std::array{particle.px(), particle.py(), particle.pz()}, massB); + auto etaParticle = particle.eta(); + + std::array ptProngs; + std::array yProngs; + std::array etaProngs; + int counter = 0; + for (const auto& daught : particle.daughters_as()) { + ptProngs[counter] = daught.pt(); + etaProngs[counter] = daught.eta(); + yProngs[counter] = RecoDecay::y(std::array{daught.px(), daught.py(), daught.pz()}, pdg->Mass(daught.pdgCode())); + counter++; + } + rowHfB0McGenReduced(flag, ptParticle, yParticle, etaParticle, + ptProngs[0], yProngs[0], etaProngs[0], + ptProngs[1], yProngs[1], etaProngs[1]); + } else if constexpr (decayChannel == DecayChannel::BplusToD0barPi) { + // B+ → D0bar π+ + if (RecoDecay::isMatchedMCGen(particlesMc, particle, Pdg::kBPlus, std::array{static_cast(Pdg::kD0), +kPiPlus}, true)) { + // Match D0bar -> π- K+ + auto candD0MC = particlesMc.rawIteratorAt(particle.daughtersIds().front()); + // Printf("Checking D0bar -> π- K+"); + if (RecoDecay::isMatchedMCGen(particlesMc, candD0MC, static_cast(Pdg::kD0), std::array{+kPiPlus, -kKPlus}, true, &sign)) { + flag = sign * BIT(hf_cand_bplus::DecayType::BplusToD0Pi); + } + } + + // save information for B+ task + if (!TESTBIT(std::abs(flag), hf_cand_bplus::DecayType::BplusToD0Pi)) { + continue; + } + + auto ptParticle = particle.pt(); + auto yParticle = RecoDecay::y(std::array{particle.px(), particle.py(), particle.pz()}, massB); + auto etaParticle = particle.eta(); + + std::array ptProngs; + std::array yProngs; + std::array etaProngs; + int counter = 0; + for (const auto& daught : particle.daughters_as()) { + ptProngs[counter] = daught.pt(); + etaProngs[counter] = daught.eta(); + yProngs[counter] = RecoDecay::y(std::array{daught.px(), daught.py(), daught.pz()}, pdg->Mass(daught.pdgCode())); + counter++; + } + rowHfBpMcGenReduced(flag, ptParticle, yParticle, etaParticle, + ptProngs[0], yProngs[0], etaProngs[0], + ptProngs[1], yProngs[1], etaProngs[1]); + } + } // gen + } + + //////////////////////////////////////////////////////////////////////////////////////////////////// + // PROCESS FUNCTIONS FOR DATA + + void processDplusPiData(aod::Collisions const& collisions, + CandsDplusFiltered const& candsC, + aod::TrackAssoc const& trackIndices, + TracksPidWithSel const& tracks, + aod::BCsWithTimestamps const& bcs) + { + // store configurables needed for B0 workflow + if (!isHfCandBhadConfigFilled) { + rowCandidateConfigB0(selectionFlagDplus.value, invMassWindowCharmHadPi.value); + isHfCandBhadConfigFilled = true; + } + + // handle normalization by the right number of collisions + hfCollisionCounter(collisions.tableSize()); + + for (const auto& collision : collisions) { + auto thisCollId = collision.globalIndex(); + auto candsCThisColl = candsC.sliceBy(candsDplusPerCollision, thisCollId); + auto trackIdsThisCollision = trackIndices.sliceBy(trackIndicesPerCollision, thisCollId); + runDataCreation(collision, candsCThisColl, trackIdsThisCollision, tracks, tracks, bcs); + } + } + PROCESS_SWITCH(HfDataCreatorCharmHadPiReduced, processDplusPiData, "Process DplusPi without MC info and without ML info", true); + + void processDplusPiDataWithMl(aod::Collisions const& collisions, + CandsDplusFilteredWithMl const& candsC, + aod::TrackAssoc const& trackIndices, + TracksPidWithSel const& tracks, + aod::BCsWithTimestamps const& bcs) + { + // store configurables needed for B0 workflow + if (!isHfCandBhadConfigFilled) { + rowCandidateConfigB0(selectionFlagDplus.value, invMassWindowCharmHadPi.value); + isHfCandBhadConfigFilled = true; + } + + // handle normalization by the right number of collisions + hfCollisionCounter(collisions.tableSize()); + + for (const auto& collision : collisions) { + auto thisCollId = collision.globalIndex(); + auto candsCThisColl = candsC.sliceBy(candsDplusPerCollisionWithMl, thisCollId); + auto trackIdsThisCollision = trackIndices.sliceBy(trackIndicesPerCollision, thisCollId); + runDataCreation(collision, candsCThisColl, trackIdsThisCollision, tracks, tracks, bcs); + } + } + PROCESS_SWITCH(HfDataCreatorCharmHadPiReduced, processDplusPiDataWithMl, "Process DplusPi without MC info and with ML info", false); + + void processD0PiData(aod::Collisions const& collisions, + CandsD0Filtered const& candsC, + aod::TrackAssoc const& trackIndices, + TracksPidWithSel const& tracks, + aod::BCsWithTimestamps const& bcs) + { + // store configurables needed for B+ workflow + if (!isHfCandBhadConfigFilled) { + rowCandidateConfigBplus(selectionFlagD0.value, selectionFlagD0bar.value, invMassWindowCharmHadPi.value); + isHfCandBhadConfigFilled = true; + } + + // handle normalization by the right number of collisions + hfCollisionCounter(collisions.tableSize()); + + for (const auto& collision : collisions) { + auto thisCollId = collision.globalIndex(); + auto candsCThisColl = candsC.sliceBy(candsD0PerCollision, thisCollId); + auto trackIdsThisCollision = trackIndices.sliceBy(trackIndicesPerCollision, thisCollId); + runDataCreation(collision, candsCThisColl, trackIdsThisCollision, tracks, tracks, bcs); + } + } + PROCESS_SWITCH(HfDataCreatorCharmHadPiReduced, processD0PiData, "Process D0Pi without MC info and without ML info", true); + + void processD0PiDataWithMl(aod::Collisions const& collisions, + CandsD0FilteredWithMl const& candsC, + aod::TrackAssoc const& trackIndices, + TracksPidWithSel const& tracks, + aod::BCsWithTimestamps const& bcs) + { + // store configurables needed for B+ workflow + if (!isHfCandBhadConfigFilled) { + rowCandidateConfigBplus(selectionFlagD0.value, selectionFlagD0bar.value, invMassWindowCharmHadPi.value); + isHfCandBhadConfigFilled = true; + } + + // handle normalization by the right number of collisions + hfCollisionCounter(collisions.tableSize()); + + for (const auto& collision : collisions) { + auto thisCollId = collision.globalIndex(); + auto candsCThisColl = candsC.sliceBy(candsD0PerCollisionWithMl, thisCollId); + auto trackIdsThisCollision = trackIndices.sliceBy(trackIndicesPerCollision, thisCollId); + runDataCreation(collision, candsCThisColl, trackIdsThisCollision, tracks, tracks, bcs); + } + } + PROCESS_SWITCH(HfDataCreatorCharmHadPiReduced, processD0PiDataWithMl, "Process D0Pi without MC info and with ML info", false); + + //////////////////////////////////////////////////////////////////////////////////////////////////// + // PROCESS FUNCTIONS FOR MC + + void processDplusPiMc(aod::Collisions const& collisions, + CandsDplusFiltered const& candsC, + aod::TrackAssoc const& trackIndices, + TracksPidWithSelAndMc const& tracks, + aod::McParticles const& particlesMc, + aod::BCsWithTimestamps const& bcs) + { + // store configurables needed for B0 workflow + if (!isHfCandBhadConfigFilled) { + rowCandidateConfigB0(selectionFlagDplus.value, invMassWindowCharmHadPi.value); + isHfCandBhadConfigFilled = true; + } + + // handle normalization by the right number of collisions + hfCollisionCounter(collisions.tableSize()); + + for (const auto& collision : collisions) { + auto thisCollId = collision.globalIndex(); + auto candsCThisColl = candsC.sliceBy(candsDplusPerCollision, thisCollId); + auto trackIdsThisCollision = trackIndices.sliceBy(trackIndicesPerCollision, thisCollId); + runDataCreation(collision, candsCThisColl, trackIdsThisCollision, tracks, particlesMc, bcs); + } + runMcGen(particlesMc); + } + PROCESS_SWITCH(HfDataCreatorCharmHadPiReduced, processDplusPiMc, "Process DplusPi with MC info and without ML info", false); + + void processDplusPiMcWithMl(aod::Collisions const& collisions, + CandsDplusFilteredWithMl const& candsC, + aod::TrackAssoc const& trackIndices, + TracksPidWithSelAndMc const& tracks, + aod::McParticles const& particlesMc, + aod::BCsWithTimestamps const& bcs) + { + // store configurables needed for B0 workflow + if (!isHfCandBhadConfigFilled) { + rowCandidateConfigB0(selectionFlagDplus.value, invMassWindowCharmHadPi.value); + isHfCandBhadConfigFilled = true; + } + + // handle normalization by the right number of collisions + hfCollisionCounter(collisions.tableSize()); + + for (const auto& collision : collisions) { + auto thisCollId = collision.globalIndex(); + auto candsCThisColl = candsC.sliceBy(candsDplusPerCollisionWithMl, thisCollId); + auto trackIdsThisCollision = trackIndices.sliceBy(trackIndicesPerCollision, thisCollId); + runDataCreation(collision, candsCThisColl, trackIdsThisCollision, tracks, particlesMc, bcs); + } + runMcGen(particlesMc); + } + PROCESS_SWITCH(HfDataCreatorCharmHadPiReduced, processDplusPiMcWithMl, "Process DplusPi with MC info and with ML info", false); + + void processD0PiMc(aod::Collisions const& collisions, + CandsD0Filtered const& candsC, + aod::TrackAssoc const& trackIndices, + TracksPidWithSelAndMc const& tracks, + aod::McParticles const& particlesMc, + aod::BCsWithTimestamps const& bcs) + { + // store configurables needed for B+ workflow + if (!isHfCandBhadConfigFilled) { + rowCandidateConfigBplus(selectionFlagD0.value, selectionFlagD0bar.value, invMassWindowCharmHadPi.value); + isHfCandBhadConfigFilled = true; + } + + // handle normalization by the right number of collisions + hfCollisionCounter(collisions.tableSize()); + + for (const auto& collision : collisions) { + auto thisCollId = collision.globalIndex(); + auto candsCThisColl = candsC.sliceBy(candsD0PerCollision, thisCollId); + auto trackIdsThisCollision = trackIndices.sliceBy(trackIndicesPerCollision, thisCollId); + runDataCreation(collision, candsCThisColl, trackIdsThisCollision, tracks, particlesMc, bcs); + } + runMcGen(particlesMc); + } + PROCESS_SWITCH(HfDataCreatorCharmHadPiReduced, processD0PiMc, "Process D0Pi with MC info and without ML info", false); + + void processD0PiMcWithMl(aod::Collisions const& collisions, + CandsD0FilteredWithMl const& candsC, + aod::TrackAssoc const& trackIndices, + TracksPidWithSelAndMc const& tracks, + aod::McParticles const& particlesMc, + aod::BCsWithTimestamps const& bcs) + { + // store configurables needed for B+ workflow + if (!isHfCandBhadConfigFilled) { + rowCandidateConfigBplus(selectionFlagD0.value, selectionFlagD0bar.value, invMassWindowCharmHadPi.value); + isHfCandBhadConfigFilled = true; + } + + // handle normalization by the right number of collisions + hfCollisionCounter(collisions.tableSize()); + + for (const auto& collision : collisions) { + auto thisCollId = collision.globalIndex(); + auto candsCThisColl = candsC.sliceBy(candsD0PerCollisionWithMl, thisCollId); + auto trackIdsThisCollision = trackIndices.sliceBy(trackIndicesPerCollision, thisCollId); + runDataCreation(collision, candsCThisColl, trackIdsThisCollision, tracks, particlesMc, bcs); + } + runMcGen(particlesMc); + } + PROCESS_SWITCH(HfDataCreatorCharmHadPiReduced, processD0PiMcWithMl, "Process D0Pi with MC info and with ML info", false); +}; // struct + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{adaptAnalysisTask(cfgc)}; +} diff --git a/PWGHF/D2H/TableProducer/dataCreatorD0PiReduced.cxx b/PWGHF/D2H/TableProducer/dataCreatorD0PiReduced.cxx deleted file mode 100644 index 9f501603cd3..00000000000 --- a/PWGHF/D2H/TableProducer/dataCreatorD0PiReduced.cxx +++ /dev/null @@ -1,568 +0,0 @@ -// 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 dataCreatorD0PiReduced.cxx -/// \brief Creation of D0-Pi pairs -/// -/// \author Antonio Palasciano , Università degli Studi di Bari -/// \author Fabrizio Grosa , CERN - -#include - -#include "CommonConstants/PhysicsConstants.h" -#include "DCAFitter/DCAFitterN.h" -#include "Framework/AnalysisTask.h" -#include "Framework/O2DatabasePDGPlugin.h" -#include "Framework/runDataProcessing.h" -#include "ReconstructionDataFormats/DCA.h" - -#include "Common/Core/trackUtilities.h" -#include "Common/DataModel/CollisionAssociationTables.h" - -#include "PWGHF/Core/HfHelper.h" -#include "PWGHF/DataModel/CandidateReconstructionTables.h" -#include "PWGHF/DataModel/CandidateSelectionTables.h" -#include "PWGHF/Utils/utilsBfieldCCDB.h" -#include "PWGHF/D2H/DataModel/ReducedDataModel.h" - -using namespace o2; -using namespace o2::analysis; -using namespace o2::aod; -using namespace o2::constants::physics; -using namespace o2::framework; -using namespace o2::framework::expressions; - -// event types -enum Event : uint8_t { - Processed = 0, - NoD0PiSelected, - D0PiSelected, - kNEvent -}; - -/// Creation of D0-Pi pairs -struct HfDataCreatorD0PiReduced { - // Produces AOD tables to store track information - Produces hfReducedCollision; - Produces hfCollisionCounter; - Produces hfTrackPion; - Produces hfTrackCovPion; - Produces hfTrackPidPion; - Produces hfCand2Prong; - Produces hfCand2ProngCov; - Produces hfCand2ProngMl; - Produces rowCandidateConfig; - Produces rowHfD0PiMcRecReduced; - Produces rowHfBpMcGenReduced; - - // vertexing - // Configurable bz{"bz", 5., "magnetic field"}; - Configurable propagateToPCA{"propagateToPCA", true, "create tracks version propagated to PCA"}; - Configurable useAbsDCA{"useAbsDCA", false, "Minimise abs. distance rather than chi2"}; - Configurable useWeightedFinalPCA{"useWeightedFinalPCA", false, "Recalculate vertex position using track covariances, effective only if useAbsDCA is true"}; - Configurable maxR{"maxR", 200., "reject PCA's above this radius"}; - Configurable maxDZIni{"maxDZIni", 4., "reject (if>0) PCA candidate if tracks DZ exceeds threshold"}; - Configurable minParamChange{"minParamChange", 1.e-3, "stop iterations if largest change of any B+ is smaller than this"}; - Configurable minRelChi2Change{"minRelChi2Change", 0.9, "stop iterations is chi2/chi2old > this"}; - // selection - Configurable usePionIsGlobalTrackWoDCA{"usePionIsGlobalTrackWoDCA", true, "check isGlobalTrackWoDCA status for pions, for Run3 studies"}; - Configurable ptPionMin{"ptPionMin", 0.5, "minimum pion pT threshold (GeV/c)"}; - Configurable> binsPtPion{"binsPtPion", std::vector{hf_cuts_single_track::vecBinsPtTrack}, "track pT bin limits for pion DCA XY pT-dependent cut"}; - Configurable> cutsTrackPionDCA{"cutsTrackPionDCA", {hf_cuts_single_track::cutsTrack[0], hf_cuts_single_track::nBinsPtTrack, hf_cuts_single_track::nCutVarsTrack, hf_cuts_single_track::labelsPtTrack, hf_cuts_single_track::labelsCutVarTrack}, "Single-track selections per pT bin for pions"}; - Configurable invMassWindowD0Pi{"invMassWindowD0Pi", 0.3, "invariant-mass window for D0Pi pair preselections (GeV/c2)"}; - Configurable selectionFlagD0{"selectionFlagD0", 1, "Selection Flag for D0"}; - Configurable selectionFlagD0bar{"selectionFlagD0bar", 1, "Selection Flag for D0bar"}; - // magnetic field setting from CCDB - Configurable ccdbUrl{"ccdbUrl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; - Configurable ccdbPathGrpMag{"ccdbPathGrpMag", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object (Run 3)"}; - - HfHelper hfHelper; - - // CCDB service - Service ccdb; - o2::base::Propagator::MatCorrType noMatCorr = o2::base::Propagator::MatCorrType::USEMatCorrNONE; - int runNumber; - - // O2DatabasePDG service - Service pdg; - - double massPi{0.}; - double massD0{0.}; - double massBplus{0.}; - double invMass2D0PiMin{0.}; - double invMass2D0PiMax{0.}; - double bz{0.}; - - bool isHfCandBplusConfigFilled = false; - - // Fitter to redo D0-vertex to get extrapolated daughter tracks (2-prong vertex filter) - o2::vertexing::DCAFitterN<2> df2; - - using TracksPidAll = soa::Join; - using TracksPIDWithSel = soa::Join; - using TracksPIDWithSelAndMc = soa::Join; - using CandsDFiltered = soa::Filtered>; - using CandsDFilteredWithMl = soa::Filtered>; - - Filter filterSelectCandidates = (aod::hf_sel_candidate_d0::isSelD0 >= selectionFlagD0 || aod::hf_sel_candidate_d0::isSelD0bar >= selectionFlagD0bar); - - Preslice candsDPerCollision = aod::track_association::collisionId; - Preslice candsDPerCollisionWithMl = aod::track_association::collisionId; - Preslice trackIndicesPerCollision = aod::track_association::collisionId; - - HistogramRegistry registry{"registry"}; - - void init(InitContext const&) - { - // invariant-mass window cut - massPi = MassPiPlus; - massD0 = MassD0; - massBplus = MassBPlus; - invMass2D0PiMin = (massBplus - invMassWindowD0Pi) * (massBplus - invMassWindowD0Pi); - invMass2D0PiMax = (massBplus + invMassWindowD0Pi) * (massBplus + invMassWindowD0Pi); - - // Initialize fitter - df2.setPropagateToPCA(propagateToPCA); - df2.setMaxR(maxR); - df2.setMaxDZIni(maxDZIni); - df2.setMinParamChange(minParamChange); - df2.setMinRelChi2Change(minRelChi2Change); - df2.setUseAbsDCA(useAbsDCA); - df2.setWeightedFinalPCA(useWeightedFinalPCA); - df2.setMatCorrType(noMatCorr); - - // Configure CCDB access - ccdb->setURL(ccdbUrl); - ccdb->setCaching(true); - ccdb->setLocalObjectValidityChecking(); - runNumber = 0; - - // histograms - constexpr int kNBinsEvents = kNEvent; - std::string labels[kNBinsEvents]; - labels[Event::Processed] = "processed"; - labels[Event::NoD0PiSelected] = "without D0Pi pairs"; - labels[Event::D0PiSelected] = "with D0Pi pairs"; - static const AxisSpec axisEvents = {kNBinsEvents, 0.5, kNBinsEvents + 0.5, ""}; - registry.add("hEvents", "Events;;entries", HistType::kTH1F, {axisEvents}); - for (int iBin = 0; iBin < kNBinsEvents; iBin++) { - registry.get(HIST("hEvents"))->GetXaxis()->SetBinLabel(iBin + 1, labels[iBin].data()); - } - - registry.add("hMassD0ToKPi", "D^{0}} candidates;inv. mass (K^{#minus} #pi^{#plus}) (GeV/#it{c}^{2});entries", {HistType::kTH1F, {{500, 0., 5.}}}); - registry.add("hPtD0", "D^{0} candidates;D^{0} candidate #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{100, 0., 10.}}}); - registry.add("hPtPion", "#pi^{#plus} candidates;#pi^{#plus} candidate #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{100, 0., 10.}}}); - registry.add("hCPAD0", "D^{0} candidates;D^{0} cosine of pointing angle;entries", {HistType::kTH1F, {{110, -1.1, 1.1}}}); - } - - /// Pion selection (D0 Pi <-- B+) - /// \param trackPion is a track with the pion hypothesis - /// \param trackParCovPion is the track parametrisation of the pion - /// \param dcaPion is the 2-D array with track DCAs of the pion - /// \param track0 is prong0 of selected D0 candidate - /// \param track1 is prong1 of selected D0 candidate - /// \param candD0 is the D0 candidate - /// \return true if trackPion passes all cuts - template - bool isPionSelected(const T1& trackPion, const T2& trackParCovPion, const T3& dcaPion, const T1& track0, const T1& track1, const T4& candD0) - { - // check isGlobalTrackWoDCA status for pions if wanted - if (usePionIsGlobalTrackWoDCA && !trackPion.isGlobalTrackWoDCA()) { - return false; - } - // minimum pT selection - if (trackParCovPion.getPt() < ptPionMin || !isSelectedTrackDCA(trackParCovPion, dcaPion)) { - return false; - } - // reject pion not compatible with D0/D0bar hypothesis - if (!((candD0.isSelD0() >= selectionFlagD0 && trackPion.sign() < 0) || (candD0.isSelD0bar() >= selectionFlagD0bar && trackPion.sign() > 0))) { - // LOGF(debug, "D0: %d, D0bar%d, sign: %d", candD0.isSelD0(), candD0.isSelD0bar(), track.sign()); - return false; - } - // reject pions that are D daughters - if (trackPion.globalIndex() == track0.globalIndex() || trackPion.globalIndex() == track1.globalIndex()) { - return false; - } - - return true; - } - - /// Single-track cuts for pions on dcaXY - /// \param trackPar is the track parametrisation - /// \param dca is the 2-D array with track DCAs - /// \return true if track passes all cuts - template - bool isSelectedTrackDCA(const T1& trackPar, const T2& dca) - { - auto pTBinTrack = findBin(binsPtPion, trackPar.getPt()); - if (pTBinTrack == -1) { - return false; - } - - if (std::abs(dca[0]) < cutsTrackPionDCA->get(pTBinTrack, "min_dcaxytoprimary")) { - return false; // minimum DCAxy - } - if (std::abs(dca[0]) > cutsTrackPionDCA->get(pTBinTrack, "max_dcaxytoprimary")) { - return false; // maximum DCAxy - } - return true; - } - - template - void runDataCreation(aod::Collision const& collision, - C const& candsD0, - aod::TrackAssoc const& trackIndices, - T const& tracks, - P const& particlesMc, - aod::BCsWithTimestamps const& bcs) - { - // helpers for ReducedTables filling - int indexHfReducedCollision = hfReducedCollision.lastIndex() + 1; - // std::map where the key is the track.globalIndex() and - // the value is the track index in the table of the selected pions - std::map selectedTracksPion; - bool fillHfReducedCollision = false; - - auto primaryVertex = getPrimaryVertex(collision); - - // Set the magnetic field from ccdb. - // The static instance of the propagator was already modified in the HFTrackIndexSkimCreator, - // but this is not true when running on Run2 data/MC already converted into AO2Ds. - auto bc = collision.bc_as(); - if (runNumber != bc.runNumber()) { - LOG(info) << ">>>>>>>>>>>> Current run number: " << runNumber; - o2::parameters::GRPMagField* grpo = ccdb->getForTimeStamp(ccdbPathGrpMag, bc.timestamp()); - if (grpo == nullptr) { - LOGF(fatal, "Run 3 GRP object (type o2::parameters::GRPMagField) is not available in CCDB for run=%d at timestamp=%llu", bc.runNumber(), bc.timestamp()); - } - o2::base::Propagator::initFieldFromGRP(grpo); - bz = o2::base::Propagator::Instance()->getNominalBz(); - LOG(info) << ">>>>>>>>>>>> Magnetic field: " << bz; - runNumber = bc.runNumber(); - } - df2.setBz(bz); - - auto thisCollId = collision.globalIndex(); - for (const auto& candD0 : candsD0) { - int indexHfCand2Prong = hfCand2Prong.lastIndex() + 1; - bool fillHfCand2Prong = false; - float invMassD0{-1.f}, invMassD0bar{-1.f}; - - if (candD0.isSelD0() >= selectionFlagD0) { - invMassD0 = hfHelper.invMassD0ToPiK(candD0); - registry.fill(HIST("hMassD0ToKPi"), invMassD0); - } - if (candD0.isSelD0bar() >= selectionFlagD0bar) { - invMassD0bar = hfHelper.invMassD0barToKPi(candD0); - registry.fill(HIST("hMassD0ToKPi"), invMassD0bar); - } - registry.fill(HIST("hPtD0"), candD0.pt()); - registry.fill(HIST("hCPAD0"), candD0.cpa()); - - auto track0 = candD0.template prong0_as(); - auto track1 = candD0.template prong1_as(); - auto trackParCov0 = getTrackParCov(track0); - auto trackParCov1 = getTrackParCov(track1); - - std::array pVec0 = {track0.px(), track0.py(), track0.pz()}; - std::array pVec1 = {track1.px(), track1.py(), track1.pz()}; - - auto dca0 = o2::dataformats::DCA(track0.dcaXY(), track0.dcaZ(), track0.cYY(), track0.cZY(), track0.cZZ()); - auto dca1 = o2::dataformats::DCA(track1.dcaXY(), track1.dcaZ(), track1.cYY(), track1.cZY(), track1.cZZ()); - - // repropagate tracks to this collision if needed - if (track0.collisionId() != thisCollId) { - trackParCov0.propagateToDCA(primaryVertex, bz, &dca0); - } - - if (track1.collisionId() != thisCollId) { - trackParCov1.propagateToDCA(primaryVertex, bz, &dca1); - } - - // --------------------------------- - // reconstruct 2-prong secondary vertex (D0) - if (df2.process(trackParCov0, trackParCov1) == 0) { - continue; - } - - const auto& secondaryVertexD0 = df2.getPCACandidate(); - // propagate the 2 prongs to the secondary vertex - trackParCov0.propagateTo(secondaryVertexD0[0], bz); - trackParCov1.propagateTo(secondaryVertexD0[0], bz); - - // update pVec of tracks - df2.getTrack(0).getPxPyPzGlo(pVec0); - df2.getTrack(1).getPxPyPzGlo(pVec1); - - // D0(bar) → π∓ K± - std::array pVecD0 = RecoDecay::pVec(pVec0, pVec1); - auto trackParCovD0 = df2.createParentTrackParCov(); - trackParCovD0.setAbsCharge(0); // to be sure - - for (const auto& trackId : trackIndices) { - auto trackPion = trackId.template track_as(); - - // apply selections on pion tracks - auto trackParCovPion = getTrackParCov(trackPion); - o2::gpu::gpustd::array dcaPion{trackPion.dcaXY(), trackPion.dcaZ()}; - std::array pVecPion = {trackPion.px(), trackPion.py(), trackPion.pz()}; - if (trackPion.collisionId() != thisCollId) { - o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, trackParCovPion, 2.f, noMatCorr, &dcaPion); - getPxPyPz(trackParCovPion, pVecPion); - } - - // apply selections on pion tracks - if (!isPionSelected(trackPion, trackParCovPion, dcaPion, track0, track1, candD0)) { - continue; - } - registry.fill(HIST("hPtPion"), trackParCovPion.getPt()); - // compute invariant mass square and apply selection - auto invMass2D0Pi = RecoDecay::m2(std::array{pVecD0, pVecPion}, std::array{massD0, massPi}); - if ((invMass2D0Pi < invMass2D0PiMin) || (invMass2D0Pi > invMass2D0PiMax)) { - continue; - } - - // fill Pion tracks table - // if information on track already stored, go to next track - if (!selectedTracksPion.count(trackPion.globalIndex())) { - hfTrackPion(trackPion.globalIndex(), indexHfReducedCollision, - trackParCovPion.getX(), trackParCovPion.getAlpha(), - trackParCovPion.getY(), trackParCovPion.getZ(), trackParCovPion.getSnp(), - trackParCovPion.getTgl(), trackParCovPion.getQ2Pt()); - hfTrackCovPion(trackParCovPion.getSigmaY2(), trackParCovPion.getSigmaZY(), trackParCovPion.getSigmaZ2(), - trackParCovPion.getSigmaSnpY(), trackParCovPion.getSigmaSnpZ(), - trackParCovPion.getSigmaSnp2(), trackParCovPion.getSigmaTglY(), trackParCovPion.getSigmaTglZ(), - trackParCovPion.getSigmaTglSnp(), trackParCovPion.getSigmaTgl2(), - trackParCovPion.getSigma1PtY(), trackParCovPion.getSigma1PtZ(), trackParCovPion.getSigma1PtSnp(), - trackParCovPion.getSigma1PtTgl(), trackParCovPion.getSigma1Pt2()); - hfTrackPidPion(trackPion.hasTPC(), trackPion.hasTOF(), - trackPion.tpcNSigmaPi(), trackPion.tofNSigmaPi()); - // add trackPion.globalIndex() to a list - // to keep memory of the pions filled in the table and avoid refilling them if they are paired to another D candidate - // and keep track of their index in hfTrackPion for McRec purposes - selectedTracksPion[trackPion.globalIndex()] = hfTrackPion.lastIndex(); - } - - if constexpr (doMc) { - // we check the MC matching to be stored - auto arrayDaughtersD0 = std::array{track0, track1}; - auto arrayDaughtersBplus = std::array{track0, track1, trackPion}; - int8_t sign{0}; - int8_t flag{0}; - int8_t debug{0}; - // B+ → D0(bar) π+ → (K+ π-) π+ - // Printf("Checking B+ → D0bar π+"); - auto indexRec = RecoDecay::getMatchedMCRec(particlesMc, arrayDaughtersBplus, Pdg::kBPlus, std::array{+kPiPlus, +kKPlus, -kPiPlus}, true, &sign, 2); - auto motherPt = -1.f; - if (indexRec > -1) { - // D0bar → K+ π- - // Printf("Checking D0bar → K+ π-"); - indexRec = RecoDecay::getMatchedMCRec(particlesMc, arrayDaughtersD0, Pdg::kD0, std::array{+kPiPlus, -kKPlus}, true, &sign, 1); - if (indexRec > -1) { - flag = sign * BIT(hf_cand_bplus::DecayType::BplusToD0Pi); - } else { - debug = 1; - LOGF(info, "WARNING: B+ decays in the expected final state but the condition on the intermediate state is not fulfilled"); - } - - auto indexMother = RecoDecay::getMother(particlesMc, trackPion.template mcParticle_as

(), Pdg::kBPlus, true); - if (indexMother >= 0) { - auto particleMother = particlesMc.rawIteratorAt(indexMother); - motherPt = particleMother.pt(); - } - } - - rowHfD0PiMcRecReduced(indexHfCand2Prong, selectedTracksPion[trackPion.globalIndex()], flag, debug, motherPt); - } - fillHfCand2Prong = true; - } // pion loop - if (fillHfCand2Prong) { // fill candD0 table only once per D0 candidate - hfCand2Prong(track0.globalIndex(), track1.globalIndex(), - indexHfReducedCollision, - trackParCovD0.getX(), trackParCovD0.getAlpha(), - trackParCovD0.getY(), trackParCovD0.getZ(), trackParCovD0.getSnp(), - trackParCovD0.getTgl(), trackParCovD0.getQ2Pt(), - candD0.xSecondaryVertex(), candD0.ySecondaryVertex(), candD0.zSecondaryVertex(), invMassD0, invMassD0bar); - hfCand2ProngCov(trackParCovD0.getSigmaY2(), trackParCovD0.getSigmaZY(), trackParCovD0.getSigmaZ2(), - trackParCovD0.getSigmaSnpY(), trackParCovD0.getSigmaSnpZ(), - trackParCovD0.getSigmaSnp2(), trackParCovD0.getSigmaTglY(), trackParCovD0.getSigmaTglZ(), - trackParCovD0.getSigmaTglSnp(), trackParCovD0.getSigmaTgl2(), - trackParCovD0.getSigma1PtY(), trackParCovD0.getSigma1PtZ(), trackParCovD0.getSigma1PtSnp(), - trackParCovD0.getSigma1PtTgl(), trackParCovD0.getSigma1Pt2()); - if constexpr (withMl) { - hfCand2ProngMl(candD0.mlProbD0()[0], candD0.mlProbD0()[1], candD0.mlProbD0()[2]); - } - fillHfReducedCollision = true; - } - } // candsD loop - registry.fill(HIST("hEvents"), 1 + Event::Processed); - if (!fillHfReducedCollision) { - registry.fill(HIST("hEvents"), 1 + Event::NoD0PiSelected); - return; - } - registry.fill(HIST("hEvents"), 1 + Event::D0PiSelected); - // fill collision table if it contains a D0Pi pair a minima - hfReducedCollision(collision.posX(), collision.posY(), collision.posZ(), - collision.covXX(), collision.covXY(), collision.covYY(), - collision.covXZ(), collision.covYZ(), collision.covZZ(), - bz); - } - - void runMcGen(aod::McParticles const& particlesMc) - { - // Match generated particles. - for (const auto& particle : particlesMc) { - int8_t sign{0}; - int8_t flag{0}; - // B+ → D0bar π+ - if (RecoDecay::isMatchedMCGen(particlesMc, particle, Pdg::kBPlus, std::array{static_cast(Pdg::kD0), +kPiPlus}, true)) { - // Match D0bar -> π- K+ - auto candD0MC = particlesMc.rawIteratorAt(particle.daughtersIds().front()); - // Printf("Checking D0bar -> π- K+"); - if (RecoDecay::isMatchedMCGen(particlesMc, candD0MC, static_cast(Pdg::kD0), std::array{+kPiPlus, -kKPlus}, true, &sign)) { - flag = sign * BIT(hf_cand_bplus::DecayType::BplusToD0Pi); - } - } - - // save information for B+ task - if (!TESTBIT(std::abs(flag), hf_cand_bplus::DecayType::BplusToD0Pi)) { - continue; - } - - auto ptParticle = particle.pt(); - auto yParticle = RecoDecay::y(std::array{particle.px(), particle.py(), particle.pz()}, massBplus); - auto etaParticle = particle.eta(); - - std::array ptProngs; - std::array yProngs; - std::array etaProngs; - int counter = 0; - for (const auto& daught : particle.daughters_as()) { - ptProngs[counter] = daught.pt(); - etaProngs[counter] = daught.eta(); - yProngs[counter] = RecoDecay::y(std::array{daught.px(), daught.py(), daught.pz()}, pdg->Mass(daught.pdgCode())); - counter++; - } - rowHfBpMcGenReduced(flag, ptParticle, yParticle, etaParticle, - ptProngs[0], yProngs[0], etaProngs[0], - ptProngs[1], yProngs[1], etaProngs[1]); - } // gen - } - - void processData(aod::Collisions const& collisions, - CandsDFiltered const& candsD0, - aod::TrackAssoc const& trackIndices, - TracksPIDWithSel const& tracks, - aod::BCsWithTimestamps const& bcs) - { - // store configurables needed for B+ workflow - if (!isHfCandBplusConfigFilled) { - rowCandidateConfig(selectionFlagD0.value, selectionFlagD0bar.value, invMassWindowD0Pi.value); - isHfCandBplusConfigFilled = true; - } - - // handle normalization by the right number of collisions - hfCollisionCounter(collisions.tableSize()); - - for (const auto& collision : collisions) { - auto thisCollId = collision.globalIndex(); - auto candsDThisColl = candsD0.sliceBy(candsDPerCollision, thisCollId); - auto trackIdsThisCollision = trackIndices.sliceBy(trackIndicesPerCollision, thisCollId); - runDataCreation(collision, candsDThisColl, trackIdsThisCollision, tracks, tracks, bcs); - } - } - PROCESS_SWITCH(HfDataCreatorD0PiReduced, processData, "Process without MC info and without ML info", true); - - void processDataWithMl(aod::Collisions const& collisions, - CandsDFilteredWithMl const& candsD0, - aod::TrackAssoc const& trackIndices, - TracksPIDWithSel const& tracks, - aod::BCsWithTimestamps const& bcs) - { - // store configurables needed for B+ workflow - if (!isHfCandBplusConfigFilled) { - rowCandidateConfig(selectionFlagD0.value, selectionFlagD0bar.value, invMassWindowD0Pi.value); - isHfCandBplusConfigFilled = true; - } - - // handle normalization by the right number of collisions - hfCollisionCounter(collisions.tableSize()); - - for (const auto& collision : collisions) { - auto thisCollId = collision.globalIndex(); - auto candsDThisColl = candsD0.sliceBy(candsDPerCollision, thisCollId); - auto trackIdsThisCollision = trackIndices.sliceBy(trackIndicesPerCollision, thisCollId); - runDataCreation(collision, candsDThisColl, trackIdsThisCollision, tracks, tracks, bcs); - } - } - PROCESS_SWITCH(HfDataCreatorD0PiReduced, processDataWithMl, "Process without MC info and with ML info", false); - - void processMc(aod::Collisions const& collisions, - CandsDFiltered const& candsD0, - aod::TrackAssoc const& trackIndices, - TracksPIDWithSelAndMc const& tracks, - aod::McParticles const& particlesMc, - aod::BCsWithTimestamps const& bcs) - { - // store configurables needed for B+ workflow - if (!isHfCandBplusConfigFilled) { - rowCandidateConfig(selectionFlagD0.value, selectionFlagD0bar.value, invMassWindowD0Pi.value); - isHfCandBplusConfigFilled = true; - } - - // handle normalization by the right number of collisions - hfCollisionCounter(collisions.tableSize()); - - for (const auto& collision : collisions) { - auto thisCollId = collision.globalIndex(); - auto candsDThisColl = candsD0.sliceBy(candsDPerCollision, thisCollId); - auto trackIdsThisCollision = trackIndices.sliceBy(trackIndicesPerCollision, thisCollId); - runDataCreation(collision, candsDThisColl, trackIdsThisCollision, tracks, particlesMc, bcs); - } - - runMcGen(particlesMc); - } - PROCESS_SWITCH(HfDataCreatorD0PiReduced, processMc, "Process with MC info and without ML info", false); - - void processMcWithMl(aod::Collisions const& collisions, - CandsDFilteredWithMl const& candsD0, - aod::TrackAssoc const& trackIndices, - TracksPIDWithSelAndMc const& tracks, - aod::McParticles const& particlesMc, - aod::BCsWithTimestamps const& bcs) - { - // store configurables needed for B+ workflow - if (!isHfCandBplusConfigFilled) { - rowCandidateConfig(selectionFlagD0.value, selectionFlagD0bar.value, invMassWindowD0Pi.value); - isHfCandBplusConfigFilled = true; - } - - // handle normalization by the right number of collisions - hfCollisionCounter(collisions.tableSize()); - - for (const auto& collision : collisions) { - auto thisCollId = collision.globalIndex(); - auto candsDThisColl = candsD0.sliceBy(candsDPerCollision, thisCollId); - auto trackIdsThisCollision = trackIndices.sliceBy(trackIndicesPerCollision, thisCollId); - runDataCreation(collision, candsDThisColl, trackIdsThisCollision, tracks, particlesMc, bcs); - } - - runMcGen(particlesMc); - } - PROCESS_SWITCH(HfDataCreatorD0PiReduced, processMcWithMl, "Process with MC info and with ML info", false); - -}; // struct - -WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) -{ - return WorkflowSpec{adaptAnalysisTask(cfgc)}; -} diff --git a/PWGHF/D2H/TableProducer/dataCreatorDplusPiReduced.cxx b/PWGHF/D2H/TableProducer/dataCreatorDplusPiReduced.cxx deleted file mode 100644 index 3387853b11a..00000000000 --- a/PWGHF/D2H/TableProducer/dataCreatorDplusPiReduced.cxx +++ /dev/null @@ -1,645 +0,0 @@ -// 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 dataCreatorDplusPiReduced.cxx -/// \brief Creation of Dplus-Pi pairs -/// -/// \author Alexandre Bigot , IPHC Strasbourg -/// \author Fabrizio Grosa , CERN - -#include - -#include "CommonConstants/PhysicsConstants.h" -#include "DCAFitter/DCAFitterN.h" -#include "Framework/AnalysisTask.h" -#include "Framework/O2DatabasePDGPlugin.h" -#include "Framework/runDataProcessing.h" -#include "ReconstructionDataFormats/DCA.h" - -#include "Common/Core/trackUtilities.h" -#include "Common/DataModel/CollisionAssociationTables.h" - -#include "PWGHF/Core/HfHelper.h" -#include "PWGHF/DataModel/CandidateReconstructionTables.h" -#include "PWGHF/DataModel/CandidateSelectionTables.h" -#include "PWGHF/Utils/utilsBfieldCCDB.h" -#include "PWGHF/D2H/DataModel/ReducedDataModel.h" - -using namespace o2; -using namespace o2::analysis; -using namespace o2::aod; -using namespace o2::constants::physics; -using namespace o2::framework; -using namespace o2::framework::expressions; - -// event types -enum Event : uint8_t { - Processed = 0, - NoDPiSelected, - DPiSelected, - kNEvent -}; - -/// Creation of Dplus-Pi pairs -struct HfDataCreatorDplusPiReduced { - // Produces AOD tables to store track information - Produces hfReducedCollision; - Produces hfCollisionCounter; - Produces hfTrackPion; - Produces hfTrackCovPion; - Produces hfTrackPidPion; - Produces hfCand3Prong; - Produces hfCand3ProngCov; - Produces hfCand3ProngMl; - Produces rowCandidateConfig; - Produces rowHfDPiMcRecReduced; - Produces rowHfDPiMcCheckReduced; - Produces rowHfB0McGenReduced; - - // vertexing - // Configurable bz{"bz", 5., "magnetic field"}; - Configurable propagateToPCA{"propagateToPCA", true, "create tracks version propagated to PCA"}; - Configurable useAbsDCA{"useAbsDCA", false, "Minimise abs. distance rather than chi2"}; - Configurable useWeightedFinalPCA{"useWeightedFinalPCA", false, "Recalculate vertex position using track covariances, effective only if useAbsDCA is true"}; - Configurable maxR{"maxR", 200., "reject PCA's above this radius"}; - Configurable maxDZIni{"maxDZIni", 4., "reject (if>0) PCA candidate if tracks DZ exceeds threshold"}; - Configurable minParamChange{"minParamChange", 1.e-3, "stop iterations if largest change of any B0 is smaller than this"}; - Configurable minRelChi2Change{"minRelChi2Change", 0.9, "stop iterations is chi2/chi2old > this"}; - // selection - Configurable usePionIsGlobalTrackWoDCA{"usePionIsGlobalTrackWoDCA", true, "check isGlobalTrackWoDCA status for pions, for Run3 studies"}; - Configurable ptPionMin{"ptPionMin", 0.5, "minimum pion pT threshold (GeV/c)"}; - Configurable> binsPtPion{"binsPtPion", std::vector{hf_cuts_single_track::vecBinsPtTrack}, "track pT bin limits for pion DCA XY pT-dependent cut"}; - Configurable> cutsTrackPionDCA{"cutsTrackPionDCA", {hf_cuts_single_track::cutsTrack[0], hf_cuts_single_track::nBinsPtTrack, hf_cuts_single_track::nCutVarsTrack, hf_cuts_single_track::labelsPtTrack, hf_cuts_single_track::labelsCutVarTrack}, "Single-track selections per pT bin for pions"}; - Configurable invMassWindowDPi{"invMassWindowDPi", 0.3, "invariant-mass window for DPi pair preselections (GeV/c2)"}; - Configurable selectionFlagD{"selectionFlagD", 1, "Selection Flag for D"}; - - // magnetic field setting from CCDB - Configurable ccdbUrl{"ccdbUrl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; - Configurable ccdbPathGrpMag{"ccdbPathGrpMag", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object (Run 3)"}; - - // MC extra - Configurable checkDecayTypeMc{"checkDecayTypeMc", false, "flag to enable MC checks on decay type"}; - - HfHelper hfHelper; - - // CCDB service - Service ccdb; - o2::base::Propagator::MatCorrType noMatCorr = o2::base::Propagator::MatCorrType::USEMatCorrNONE; - int runNumber; - - // O2DatabasePDG service - Service pdg; - - double massPi{0.}; - double massD{0.}; - double massB0{0.}; - double invMass2DPiMin{0.}; - double invMass2DPiMax{0.}; - double bz{0.}; - - bool isHfCandB0ConfigFilled = false; - - // Fitter to redo D-vertex to get extrapolated daughter tracks (3-prong vertex filter) - o2::vertexing::DCAFitterN<3> df3; - - using TracksPidAll = soa::Join; - using TracksPIDWithSel = soa::Join; - using TracksPIDWithSelAndMc = soa::Join; - using CandsDFiltered = soa::Filtered>; - using CandsDFilteredWithMl = soa::Filtered>; - - Filter filterSelectCandidates = (aod::hf_sel_candidate_dplus::isSelDplusToPiKPi >= selectionFlagD); - - Preslice candsDPerCollision = aod::track_association::collisionId; - Preslice candsDPerCollisionWithMl = aod::track_association::collisionId; - Preslice trackIndicesPerCollision = aod::track_association::collisionId; - - HistogramRegistry registry{"registry"}; - - void init(InitContext const&) - { - // invariant-mass window cut - massPi = MassPiPlus; - massD = MassDMinus; - massB0 = MassB0; - invMass2DPiMin = (massB0 - invMassWindowDPi) * (massB0 - invMassWindowDPi); - invMass2DPiMax = (massB0 + invMassWindowDPi) * (massB0 + invMassWindowDPi); - - // Initialize fitter - df3.setPropagateToPCA(propagateToPCA); - df3.setMaxR(maxR); - df3.setMaxDZIni(maxDZIni); - df3.setMinParamChange(minParamChange); - df3.setMinRelChi2Change(minRelChi2Change); - df3.setUseAbsDCA(useAbsDCA); - df3.setWeightedFinalPCA(useWeightedFinalPCA); - df3.setMatCorrType(noMatCorr); - - // Configure CCDB access - ccdb->setURL(ccdbUrl); - ccdb->setCaching(true); - ccdb->setLocalObjectValidityChecking(); - runNumber = 0; - - // histograms - constexpr int kNBinsEvents = kNEvent; - std::string labels[kNBinsEvents]; - labels[Event::Processed] = "processed"; - labels[Event::NoDPiSelected] = "without DPi pairs"; - labels[Event::DPiSelected] = "with DPi pairs"; - static const AxisSpec axisEvents = {kNBinsEvents, 0.5, kNBinsEvents + 0.5, ""}; - registry.add("hEvents", "Events;;entries", HistType::kTH1F, {axisEvents}); - for (int iBin = 0; iBin < kNBinsEvents; iBin++) { - registry.get(HIST("hEvents"))->GetXaxis()->SetBinLabel(iBin + 1, labels[iBin].data()); - } - - registry.add("hMassDToPiKPi", "D^{#minus} candidates;inv. mass (p^{#minus} K^{#plus} #pi^{#minus}) (GeV/#it{c}^{2});entries", {HistType::kTH1F, {{500, 0., 5.}}}); - registry.add("hPtD", "D^{#minus} candidates;D^{#minus} candidate #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{100, 0., 10.}}}); - registry.add("hPtPion", "#pi^{#plus} candidates;#pi^{#plus} candidate #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{100, 0., 10.}}}); - registry.add("hCPAD", "D^{#minus} candidates;D^{#minus} cosine of pointing angle;entries", {HistType::kTH1F, {{110, -1.1, 1.1}}}); - } - - /// Pion selection (D Pi <-- B0) - /// \param trackPion is a track with the pion hypothesis - /// \param trackParCovPion is the track parametrisation of the pion - /// \param dcaPion is the 2-D array with track DCAs of the pion - /// \param track0 is prong0 of selected D candidate - /// \param track1 is prong1 of selected D candidate - /// \param track2 is prong2 of selected D candidate - /// \return true if trackPion passes all cuts - template - bool isPionSelected(const T1& trackPion, const T2& trackParCovPion, const T3& dcaPion, const T1& track0, const T1& track1, const T1& track2) - { - // check isGlobalTrackWoDCA status for pions if wanted - if (usePionIsGlobalTrackWoDCA && !trackPion.isGlobalTrackWoDCA()) { - return false; - } - // minimum pT selection - if (trackParCovPion.getPt() < ptPionMin || !isSelectedTrackDCA(trackParCovPion, dcaPion)) { - return false; - } - // reject pions that are D daughters - if (trackPion.globalIndex() == track0.globalIndex() || trackPion.globalIndex() == track1.globalIndex() || trackPion.globalIndex() == track2.globalIndex()) { - return false; - } - // reject pi D with same sign as D - if (trackPion.sign() * track0.sign() > 0) { - return false; - } - return true; - } - - /// Single-track cuts for pions on dcaXY - /// \param trackPar is the track parametrisation - /// \param dca is the 2-D array with track DCAs - /// \return true if track passes all cuts - template - bool isSelectedTrackDCA(const T1& trackPar, const T2& dca) - { - auto pTBinTrack = findBin(binsPtPion, trackPar.getPt()); - if (pTBinTrack == -1) { - return false; - } - - if (std::abs(dca[0]) < cutsTrackPionDCA->get(pTBinTrack, "min_dcaxytoprimary")) { - return false; // minimum DCAxy - } - if (std::abs(dca[0]) > cutsTrackPionDCA->get(pTBinTrack, "max_dcaxytoprimary")) { - return false; // maximum DCAxy - } - return true; - } - - template - void runDataCreation(aod::Collision const& collision, - C const& candsD, - aod::TrackAssoc const& trackIndices, - T const& tracks, - P const& particlesMc, - aod::BCsWithTimestamps const& bcs) - { - // helpers for ReducedTables filling - int indexHfReducedCollision = hfReducedCollision.lastIndex() + 1; - // std::map where the key is the track.globalIndex() and - // the value is the track index in the table of the selected pions - std::map selectedTracksPion; - bool fillHfReducedCollision = false; - - auto primaryVertex = getPrimaryVertex(collision); - - // Set the magnetic field from ccdb. - // The static instance of the propagator was already modified in the HFTrackIndexSkimCreator, - // but this is not true when running on Run2 data/MC already converted into AO2Ds. - auto bc = collision.bc_as(); - if (runNumber != bc.runNumber()) { - LOG(info) << ">>>>>>>>>>>> Current run number: " << runNumber; - o2::parameters::GRPMagField* grpo = ccdb->getForTimeStamp(ccdbPathGrpMag, bc.timestamp()); - if (grpo == nullptr) { - LOGF(fatal, "Run 3 GRP object (type o2::parameters::GRPMagField) is not available in CCDB for run=%d at timestamp=%llu", bc.runNumber(), bc.timestamp()); - } - o2::base::Propagator::initFieldFromGRP(grpo); - bz = o2::base::Propagator::Instance()->getNominalBz(); - LOG(info) << ">>>>>>>>>>>> Magnetic field: " << bz; - runNumber = bc.runNumber(); - } - df3.setBz(bz); - - auto thisCollId = collision.globalIndex(); - for (const auto& candD : candsD) { - int indexHfCand3Prong = hfCand3Prong.lastIndex() + 1; - bool fillHfCand3Prong = false; - float invMassD = hfHelper.invMassDplusToPiKPi(candD); - - registry.fill(HIST("hMassDToPiKPi"), invMassD); - registry.fill(HIST("hPtD"), candD.pt()); - registry.fill(HIST("hCPAD"), candD.cpa()); - - // track0 <-> pi, track1 <-> K, track2 <-> pi - auto track0 = candD.template prong0_as(); - auto track1 = candD.template prong1_as(); - auto track2 = candD.template prong2_as(); - auto trackParCov0 = getTrackParCov(track0); - auto trackParCov1 = getTrackParCov(track1); - auto trackParCov2 = getTrackParCov(track2); - - std::array pVec0 = {track0.px(), track0.py(), track0.pz()}; - std::array pVec1 = {track1.px(), track1.py(), track1.pz()}; - std::array pVec2 = {track2.px(), track2.py(), track2.pz()}; - - auto dca0 = o2::dataformats::DCA(track0.dcaXY(), track0.dcaZ(), track0.cYY(), track0.cZY(), track0.cZZ()); - auto dca1 = o2::dataformats::DCA(track1.dcaXY(), track1.dcaZ(), track1.cYY(), track1.cZY(), track1.cZZ()); - auto dca2 = o2::dataformats::DCA(track2.dcaXY(), track2.dcaZ(), track2.cYY(), track2.cZY(), track2.cZZ()); - - // repropagate tracks to this collision if needed - if (track0.collisionId() != thisCollId) { - trackParCov0.propagateToDCA(primaryVertex, bz, &dca0); - } - - if (track1.collisionId() != thisCollId) { - trackParCov1.propagateToDCA(primaryVertex, bz, &dca1); - } - - if (track2.collisionId() != thisCollId) { - trackParCov2.propagateToDCA(primaryVertex, bz, &dca2); - } - - // --------------------------------- - // reconstruct 3-prong secondary vertex (D±) - if (df3.process(trackParCov0, trackParCov1, trackParCov2) == 0) { - continue; - } - - const auto& secondaryVertexD = df3.getPCACandidate(); - // propagate the 3 prongs to the secondary vertex - trackParCov0.propagateTo(secondaryVertexD[0], bz); - trackParCov1.propagateTo(secondaryVertexD[0], bz); - trackParCov2.propagateTo(secondaryVertexD[0], bz); - - // update pVec of tracks - df3.getTrack(0).getPxPyPzGlo(pVec0); - df3.getTrack(1).getPxPyPzGlo(pVec1); - df3.getTrack(2).getPxPyPzGlo(pVec2); - - // D∓ → π∓ K± π∓ - std::array pVecD = RecoDecay::pVec(pVec0, pVec1, pVec2); - auto trackParCovD = df3.createParentTrackParCov(); - int charge = (track1.signed1Pt() > 0) ? -1 : 1; - trackParCovD.setAbsCharge(charge); // to be sure - - for (const auto& trackId : trackIndices) { - auto trackPion = trackId.template track_as(); - - // apply selections on pion tracks - auto trackParCovPion = getTrackParCov(trackPion); - o2::gpu::gpustd::array dcaPion{trackPion.dcaXY(), trackPion.dcaZ()}; - std::array pVecPion = {trackPion.px(), trackPion.py(), trackPion.pz()}; - if (trackPion.collisionId() != thisCollId) { - o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, trackParCovPion, 2.f, noMatCorr, &dcaPion); - getPxPyPz(trackParCovPion, pVecPion); - } - - // apply selections on pion tracks - if (!isPionSelected(trackPion, trackParCovPion, dcaPion, track0, track1, track2)) { - continue; - } - registry.fill(HIST("hPtPion"), trackParCovPion.getPt()); - // compute invariant mass square and apply selection - auto invMass2DPi = RecoDecay::m2(std::array{pVecD, pVecPion}, std::array{massD, massPi}); - if ((invMass2DPi < invMass2DPiMin) || (invMass2DPi > invMass2DPiMax)) { - continue; - } - - // fill Pion tracks table - // if information on track already stored, go to next track - if (!selectedTracksPion.count(trackPion.globalIndex())) { - hfTrackPion(trackPion.globalIndex(), indexHfReducedCollision, - trackParCovPion.getX(), trackParCovPion.getAlpha(), - trackParCovPion.getY(), trackParCovPion.getZ(), trackParCovPion.getSnp(), - trackParCovPion.getTgl(), trackParCovPion.getQ2Pt()); - hfTrackCovPion(trackParCovPion.getSigmaY2(), trackParCovPion.getSigmaZY(), trackParCovPion.getSigmaZ2(), - trackParCovPion.getSigmaSnpY(), trackParCovPion.getSigmaSnpZ(), - trackParCovPion.getSigmaSnp2(), trackParCovPion.getSigmaTglY(), trackParCovPion.getSigmaTglZ(), - trackParCovPion.getSigmaTglSnp(), trackParCovPion.getSigmaTgl2(), - trackParCovPion.getSigma1PtY(), trackParCovPion.getSigma1PtZ(), trackParCovPion.getSigma1PtSnp(), - trackParCovPion.getSigma1PtTgl(), trackParCovPion.getSigma1Pt2()); - hfTrackPidPion(trackPion.hasTPC(), trackPion.hasTOF(), - trackPion.tpcNSigmaPi(), trackPion.tofNSigmaPi()); - // add trackPion.globalIndex() to a list - // to keep memory of the pions filled in the table and avoid refilling them if they are paired to another D candidate - // and keep track of their index in hfTrackPion for McRec purposes - selectedTracksPion[trackPion.globalIndex()] = hfTrackPion.lastIndex(); - } - - if constexpr (doMc) { - // we check the MC matching to be stored - auto arrayDaughtersD = std::array{track0, track1, track2}; - auto arrayDaughtersB0 = std::array{track0, track1, track2, trackPion}; - int8_t sign{0}; - int8_t flag{0}; - int8_t debug{0}; - int pdgCodeBeautyMother{-1}; - int pdgCodeCharmMother{-1}; - int pdgCodeProng0{0}; - int pdgCodeProng1{0}; - int pdgCodeProng2{0}; - int pdgCodeProng3{0}; - // B0 → D- π+ → (π- K+ π-) π+ - auto indexRec = RecoDecay::getMatchedMCRec(particlesMc, arrayDaughtersB0, Pdg::kB0, std::array{-kPiPlus, +kKPlus, -kPiPlus, +kPiPlus}, true, &sign, 3); - auto motherPt = -1.f; - if (indexRec > -1) { - // D- → π- K+ π- - // Printf("Checking D- → π- K+ π-"); - indexRec = RecoDecay::getMatchedMCRec(particlesMc, arrayDaughtersD, Pdg::kDMinus, std::array{-kPiPlus, +kKPlus, -kPiPlus}, true, &sign, 2); - if (indexRec > -1) { - flag = sign * BIT(hf_cand_b0::DecayTypeMc::B0ToDplusPiToPiKPiPi); - } else { - debug = 1; - LOGF(debug, "B0 decays in the expected final state but the condition on the intermediate state is not fulfilled"); - } - - auto indexMother = RecoDecay::getMother(particlesMc, trackPion.template mcParticle_as

(), Pdg::kB0, true); - if (indexMother >= 0) { - auto particleMother = particlesMc.rawIteratorAt(indexMother); - motherPt = particleMother.pt(); - } - } - // B0 → Ds- π+ → (K- K+ π-) π+ - if (!flag && checkDecayTypeMc) { - indexRec = RecoDecay::getMatchedMCRec(particlesMc, arrayDaughtersB0, Pdg::kB0, std::array{-kKPlus, +kKPlus, -kPiPlus, +kPiPlus}, true, &sign, 3); - if (indexRec > -1) { - // Ds- → K- K+ π- - indexRec = RecoDecay::getMatchedMCRec(particlesMc, arrayDaughtersD, -Pdg::kDS, std::array{-kKPlus, +kKPlus, -kPiPlus}, true, &sign, 2); - if (indexRec > -1) { - flag = sign * BIT(hf_cand_b0::DecayTypeMc::B0ToDsPiToKKPiPi); - } - } - } - // Partly reconstructed decays, i.e. the 4 prongs have a common b-hadron ancestor - // convention: final state particles are prong0,1,2,3 - if (!flag && checkDecayTypeMc) { - auto particleProng0 = arrayDaughtersB0[0].mcParticle(); - auto particleProng1 = arrayDaughtersB0[1].mcParticle(); - auto particleProng2 = arrayDaughtersB0[2].mcParticle(); - auto particleProng3 = arrayDaughtersB0[3].mcParticle(); - // b-hadron hypothesis - std::array bHadronMotherHypos = {Pdg::kB0, Pdg::kBS, Pdg::kLambdaB0}; - // c-hadron hypothesis - std::array cHadronMotherHypos = {Pdg::kDPlus, Pdg::kDS, Pdg::kDStar}; - - 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); - int index3Mother = RecoDecay::getMother(particlesMc, particleProng3, bHadronMotherHypo, true); - - // look for common b-hadron ancestor - if (index0Mother > -1 && index1Mother > -1 && index2Mother > -1 && index3Mother > -1) { - if (index0Mother == index1Mother && index1Mother == index2Mother && index2Mother == index3Mother) { - flag = BIT(hf_cand_b0::DecayTypeMc::PartlyRecoDecay); - pdgCodeBeautyMother = particlesMc.rawIteratorAt(index0Mother).pdgCode(); - pdgCodeCharmMother = 0; - pdgCodeProng0 = particleProng0.pdgCode(); - pdgCodeProng1 = particleProng1.pdgCode(); - pdgCodeProng2 = particleProng2.pdgCode(); - pdgCodeProng3 = particleProng3.pdgCode(); - // look for common c-hadron mother among prongs 0, 1 and 2 - for (const auto& cHadronMotherHypo : cHadronMotherHypos) { - int8_t depthMax = 2; - if (cHadronMotherHypo == Pdg::kDStar) { // to include D* -> D π0/γ and D* -> D0 π - depthMax += 1; - } - int index0CharmMother = RecoDecay::getMother(particlesMc, particleProng0, cHadronMotherHypo, true, &sign, depthMax); - int index1CharmMother = RecoDecay::getMother(particlesMc, particleProng1, cHadronMotherHypo, true, &sign, depthMax); - int index2CharmMother = RecoDecay::getMother(particlesMc, particleProng2, cHadronMotherHypo, true, &sign, depthMax); - if (index0CharmMother > -1 && index1CharmMother > -1 && index2CharmMother > -1) { - if (index0CharmMother == index1CharmMother && index1CharmMother == index2CharmMother) { - // pdgCodeCharmMother = - // Pdg::kDPlus (if D+ is the mother and does not come from D*+) - // Pdg::kDPlus + Pdg::kDStar (if D+ is the mother and D*+ -> D+ π0/γ) - // Pdg::kDStar (if D*+ is the mother and D*+ -> D0 π+) - // Pdg::kDS (if Ds is the mother) - pdgCodeCharmMother += std::abs(particlesMc.rawIteratorAt(index0CharmMother).pdgCode()); - } - } - } - break; - } - } - } - } - - rowHfDPiMcRecReduced(indexHfCand3Prong, selectedTracksPion[trackPion.globalIndex()], flag, debug, motherPt); - if (checkDecayTypeMc) { - rowHfDPiMcCheckReduced(pdgCodeBeautyMother, pdgCodeCharmMother, pdgCodeProng0, pdgCodeProng1, pdgCodeProng2, pdgCodeProng3); - } - } - fillHfCand3Prong = true; - } // pion loop - if (fillHfCand3Prong) { // fill candDplus table only once per D candidate - hfCand3Prong(track0.globalIndex(), track1.globalIndex(), track2.globalIndex(), - indexHfReducedCollision, - trackParCovD.getX(), trackParCovD.getAlpha(), - trackParCovD.getY(), trackParCovD.getZ(), trackParCovD.getSnp(), - trackParCovD.getTgl(), trackParCovD.getQ2Pt(), - candD.xSecondaryVertex(), candD.ySecondaryVertex(), candD.zSecondaryVertex(), invMassD); - hfCand3ProngCov(trackParCovD.getSigmaY2(), trackParCovD.getSigmaZY(), trackParCovD.getSigmaZ2(), - trackParCovD.getSigmaSnpY(), trackParCovD.getSigmaSnpZ(), - trackParCovD.getSigmaSnp2(), trackParCovD.getSigmaTglY(), trackParCovD.getSigmaTglZ(), - trackParCovD.getSigmaTglSnp(), trackParCovD.getSigmaTgl2(), - trackParCovD.getSigma1PtY(), trackParCovD.getSigma1PtZ(), trackParCovD.getSigma1PtSnp(), - trackParCovD.getSigma1PtTgl(), trackParCovD.getSigma1Pt2()); - if constexpr (withMl) { - hfCand3ProngMl(candD.mlProbDplusToPiKPi()[0], candD.mlProbDplusToPiKPi()[1], candD.mlProbDplusToPiKPi()[2]); - } - fillHfReducedCollision = true; - } - } // candsD loop - - registry.fill(HIST("hEvents"), 1 + Event::Processed); - if (!fillHfReducedCollision) { - registry.fill(HIST("hEvents"), 1 + Event::NoDPiSelected); - return; - } - registry.fill(HIST("hEvents"), 1 + Event::DPiSelected); - // fill collision table if it contains a DPi pair a minima - hfReducedCollision(collision.posX(), collision.posY(), collision.posZ(), - collision.covXX(), collision.covXY(), collision.covYY(), - collision.covXZ(), collision.covYZ(), collision.covZZ(), - bz); - } - - void runMcGen(aod::McParticles const& particlesMc) - { - // Match generated particles. - for (const auto& particle : particlesMc) { - int8_t sign{0}; - int8_t flag{0}; - // B0 → D- π+ - if (RecoDecay::isMatchedMCGen(particlesMc, particle, Pdg::kB0, std::array{-static_cast(Pdg::kDPlus), +kPiPlus}, true)) { - // Match D- -> π- K+ π- - auto candDMC = particlesMc.rawIteratorAt(particle.daughtersIds().front()); - // Printf("Checking D- -> π- K+ π-"); - if (RecoDecay::isMatchedMCGen(particlesMc, candDMC, -static_cast(Pdg::kDPlus), std::array{-kPiPlus, +kKPlus, -kPiPlus}, true, &sign)) { - flag = sign * BIT(hf_cand_b0::DecayType::B0ToDPi); - } - } - - // save information for B0 task - if (!TESTBIT(std::abs(flag), hf_cand_b0::DecayType::B0ToDPi)) { - continue; - } - - auto ptParticle = particle.pt(); - auto yParticle = RecoDecay::y(std::array{particle.px(), particle.py(), particle.pz()}, massB0); - auto etaParticle = particle.eta(); - - std::array ptProngs; - std::array yProngs; - std::array etaProngs; - int counter = 0; - for (const auto& daught : particle.daughters_as()) { - ptProngs[counter] = daught.pt(); - etaProngs[counter] = daught.eta(); - yProngs[counter] = RecoDecay::y(std::array{daught.px(), daught.py(), daught.pz()}, pdg->Mass(daught.pdgCode())); - counter++; - } - rowHfB0McGenReduced(flag, ptParticle, yParticle, etaParticle, - ptProngs[0], yProngs[0], etaProngs[0], - ptProngs[1], yProngs[1], etaProngs[1]); - } // gen - } - - void processData(aod::Collisions const& collisions, - CandsDFiltered const& candsD, - aod::TrackAssoc const& trackIndices, - TracksPIDWithSel const& tracks, - aod::BCsWithTimestamps const& bcs) - { - // store configurables needed for B0 workflow - if (!isHfCandB0ConfigFilled) { - rowCandidateConfig(selectionFlagD.value, invMassWindowDPi.value); - isHfCandB0ConfigFilled = true; - } - - // handle normalization by the right number of collisions - hfCollisionCounter(collisions.tableSize()); - - for (const auto& collision : collisions) { - auto thisCollId = collision.globalIndex(); - auto candsDThisColl = candsD.sliceBy(candsDPerCollision, thisCollId); - auto trackIdsThisCollision = trackIndices.sliceBy(trackIndicesPerCollision, thisCollId); - runDataCreation(collision, candsDThisColl, trackIdsThisCollision, tracks, tracks, bcs); - } - } - PROCESS_SWITCH(HfDataCreatorDplusPiReduced, processData, "Process without MC info and without ML info", true); - - void processDataWithMl(aod::Collisions const& collisions, - CandsDFilteredWithMl const& candsD, - aod::TrackAssoc const& trackIndices, - TracksPIDWithSel const& tracks, - aod::BCsWithTimestamps const& bcs) - { - // store configurables needed for B0 workflow - if (!isHfCandB0ConfigFilled) { - rowCandidateConfig(selectionFlagD.value, invMassWindowDPi.value); - isHfCandB0ConfigFilled = true; - } - - // handle normalization by the right number of collisions - hfCollisionCounter(collisions.tableSize()); - - for (const auto& collision : collisions) { - auto thisCollId = collision.globalIndex(); - auto candsDThisColl = candsD.sliceBy(candsDPerCollisionWithMl, thisCollId); - auto trackIdsThisCollision = trackIndices.sliceBy(trackIndicesPerCollision, thisCollId); - runDataCreation(collision, candsDThisColl, trackIdsThisCollision, tracks, tracks, bcs); - } - } - PROCESS_SWITCH(HfDataCreatorDplusPiReduced, processDataWithMl, "Process without MC info and with ML info", false); - - void processMc(aod::Collisions const& collisions, - CandsDFiltered const& candsD, - aod::TrackAssoc const& trackIndices, - TracksPIDWithSelAndMc const& tracks, - aod::McParticles const& particlesMc, - aod::BCsWithTimestamps const& bcs) - { - // store configurables needed for B0 workflow - if (!isHfCandB0ConfigFilled) { - rowCandidateConfig(selectionFlagD.value, invMassWindowDPi.value); - isHfCandB0ConfigFilled = true; - } - - // handle normalization by the right number of collisions - hfCollisionCounter(collisions.tableSize()); - - for (const auto& collision : collisions) { - auto thisCollId = collision.globalIndex(); - auto candsDThisColl = candsD.sliceBy(candsDPerCollision, thisCollId); - auto trackIdsThisCollision = trackIndices.sliceBy(trackIndicesPerCollision, thisCollId); - runDataCreation(collision, candsDThisColl, trackIdsThisCollision, tracks, particlesMc, bcs); - } - runMcGen(particlesMc); - } - PROCESS_SWITCH(HfDataCreatorDplusPiReduced, processMc, "Process with MC info and without ML info", false); - - void processMcWithMl(aod::Collisions const& collisions, - CandsDFilteredWithMl const& candsD, - aod::TrackAssoc const& trackIndices, - TracksPIDWithSelAndMc const& tracks, - aod::McParticles const& particlesMc, - aod::BCsWithTimestamps const& bcs) - { - // store configurables needed for B0 workflow - if (!isHfCandB0ConfigFilled) { - rowCandidateConfig(selectionFlagD.value, invMassWindowDPi.value); - isHfCandB0ConfigFilled = true; - } - - // handle normalization by the right number of collisions - hfCollisionCounter(collisions.tableSize()); - - for (const auto& collision : collisions) { - auto thisCollId = collision.globalIndex(); - auto candsDThisColl = candsD.sliceBy(candsDPerCollisionWithMl, thisCollId); - auto trackIdsThisCollision = trackIndices.sliceBy(trackIndicesPerCollision, thisCollId); - runDataCreation(collision, candsDThisColl, trackIdsThisCollision, tracks, particlesMc, bcs); - } - runMcGen(particlesMc); - } - PROCESS_SWITCH(HfDataCreatorDplusPiReduced, processMcWithMl, "Process with MC info and with ML info", false); -}; // struct - -WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) -{ - return WorkflowSpec{adaptAnalysisTask(cfgc)}; -} From e1c5091f7c5cfe4650799f42a944c4310e0551fe Mon Sep 17 00:00:00 2001 From: Fabrizio Grosa Date: Thu, 29 Feb 2024 15:50:40 +0100 Subject: [PATCH 2/5] Reorder init and remove unwanted columns in collision table for Ds resonances --- .../dataCreatorCharmHadPiReduced.cxx | 116 +++++++++++------- .../TableProducer/dataCreatorDV0Reduced.cxx | 7 +- 2 files changed, 74 insertions(+), 49 deletions(-) diff --git a/PWGHF/D2H/TableProducer/dataCreatorCharmHadPiReduced.cxx b/PWGHF/D2H/TableProducer/dataCreatorCharmHadPiReduced.cxx index 6b8d018e874..be02117f4fc 100644 --- a/PWGHF/D2H/TableProducer/dataCreatorCharmHadPiReduced.cxx +++ b/PWGHF/D2H/TableProducer/dataCreatorCharmHadPiReduced.cxx @@ -152,6 +152,50 @@ struct HfDataCreatorCharmHadPiReduced { void init(InitContext const&) { + std::array doProcess = {doprocessDplusPiData, doprocessDplusPiDataWithMl, doprocessDplusPiMc, doprocessDplusPiMcWithMl, doprocessD0PiData, doprocessD0PiDataWithMl, doprocessD0PiMc, doprocessD0PiMcWithMl}; + if (std::accumulate(doProcess.begin(), doProcess.end(), 0) != 0) { + LOGP(fatal, "One and only one process function can be enabled at a time, please fix your configuration!"); + } + + // invariant-mass window cut + massPi = MassPiPlus; + if (doprocessDplusPiData || doprocessDplusPiDataWithMl || doprocessDplusPiMc || doprocessDplusPiMcWithMl) { + massC = MassDMinus; + massB = MassB0; + } else if (doprocessD0PiData || doprocessD0PiDataWithMl || doprocessD0PiMc || doprocessD0PiMcWithMl) { + massC = MassD0; + massB = MassBPlus; + } + invMass2ChHadPiMin = (massB - invMassWindowCharmHadPi) * (massB - invMassWindowCharmHadPi); + invMass2ChHadPiMax = (massB + invMassWindowCharmHadPi) * (massB + invMassWindowCharmHadPi); + + // Initialize fitter + if (doprocessDplusPiData || doprocessDplusPiDataWithMl || doprocessDplusPiMc || doprocessDplusPiMcWithMl) { + df3.setPropagateToPCA(propagateToPCA); + df3.setMaxR(maxR); + df3.setMaxDZIni(maxDZIni); + df3.setMinParamChange(minParamChange); + df3.setMinRelChi2Change(minRelChi2Change); + df3.setUseAbsDCA(useAbsDCA); + df3.setWeightedFinalPCA(useWeightedFinalPCA); + df3.setMatCorrType(noMatCorr); + } else if (doprocessD0PiData || doprocessD0PiDataWithMl || doprocessD0PiMc || doprocessD0PiMcWithMl) { + df2.setPropagateToPCA(propagateToPCA); + df2.setMaxR(maxR); + df2.setMaxDZIni(maxDZIni); + df2.setMinParamChange(minParamChange); + df2.setMinRelChi2Change(minRelChi2Change); + df2.setUseAbsDCA(useAbsDCA); + df2.setWeightedFinalPCA(useWeightedFinalPCA); + df2.setMatCorrType(noMatCorr); + } + + // Configure CCDB access + ccdb->setURL(ccdbUrl); + ccdb->setCaching(true); + ccdb->setLocalObjectValidityChecking(); + runNumber = 0; + // histograms constexpr int kNBinsEvents = kNEvent; std::string labels[kNBinsEvents]; @@ -164,48 +208,30 @@ struct HfDataCreatorCharmHadPiReduced { registry.get(HIST("hEvents"))->GetXaxis()->SetBinLabel(iBin + 1, labels[iBin].data()); } - registry.add("hMassCharmHad0", "D^{#minus} candidates;inv. mass (p^{#minus} K^{#plus} #pi^{#minus}) (GeV/#it{c}^{2});entries", {HistType::kTH1F, {{500, 0., 5.}}}); - registry.add("hMassCharmHad1", "D^{#minus} candidates;inv. mass (p^{#minus} K^{#plus} #pi^{#minus}) (GeV/#it{c}^{2});entries", {HistType::kTH1F, {{500, 0., 5.}}}); - registry.add("hPtCharmHad", "D^{#minus} candidates;D^{#minus} candidate #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{100, 0., 10.}}}); - registry.add("hPtPion", "#pi^{#plus} candidates;#pi^{#plus} candidate #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{100, 0., 10.}}}); - registry.add("hCpaCharmHad", "D^{#minus} candidates;D^{#minus} cosine of pointing angle;entries", {HistType::kTH1F, {{110, -1.1, 1.1}}}); - - // Initialize fitter - df2.setPropagateToPCA(propagateToPCA); - df2.setMaxR(maxR); - df2.setMaxDZIni(maxDZIni); - df2.setMinParamChange(minParamChange); - df2.setMinRelChi2Change(minRelChi2Change); - df2.setUseAbsDCA(useAbsDCA); - df2.setWeightedFinalPCA(useWeightedFinalPCA); - df2.setMatCorrType(noMatCorr); - - df3.setPropagateToPCA(propagateToPCA); - df3.setMaxR(maxR); - df3.setMaxDZIni(maxDZIni); - df3.setMinParamChange(minParamChange); - df3.setMinRelChi2Change(minRelChi2Change); - df3.setUseAbsDCA(useAbsDCA); - df3.setWeightedFinalPCA(useWeightedFinalPCA); - df3.setMatCorrType(noMatCorr); - - // Configure CCDB access - ccdb->setURL(ccdbUrl); - ccdb->setCaching(true); - ccdb->setLocalObjectValidityChecking(); - runNumber = 0; - - // invariant-mass window cut - massPi = MassPiPlus; + std::string charmHadInvMassTitle = ""; + std::string charmHadTitle0 = ""; + std::string charmHadTitle1 = ""; + std::string histMassTitle0 = ""; + std::string histMassTitle1 = ""; if (doprocessDplusPiData || doprocessDplusPiDataWithMl || doprocessDplusPiMc || doprocessDplusPiMcWithMl) { - massC = MassDMinus; - massB = MassB0; + charmHadTitle0 = "D^{#plus}"; + histMassTitle0 = "Dplus"; + charmHadInvMassTitle = "#it{M}(K#pi#pi)"; } else if (doprocessD0PiData || doprocessD0PiDataWithMl || doprocessD0PiMc || doprocessD0PiMcWithMl) { - massC = MassD0; - massB = MassBPlus; + charmHadTitle0 = "D^{0}"; + charmHadTitle1 = "#overline{D}^{0}"; + histMassTitle0 = "D0"; + histMassTitle1 = "D0bar"; + charmHadInvMassTitle = "#it{M}(K#pi)"; } - invMass2ChHadPiMin = (massB - invMassWindowCharmHadPi) * (massB - invMassWindowCharmHadPi); - invMass2ChHadPiMax = (massB + invMassWindowCharmHadPi) * (massB + invMassWindowCharmHadPi); + + registry.add(Form("hMass%s", histMassTitle0.data()), Form("%s candidates; %s (GeV/#it{c}^{2});entries", charmHadTitle0.data(), charmHadInvMassTitle.data()), {HistType::kTH1F, {{500, 0., 5.}}}); + if (doprocessD0PiData || doprocessD0PiDataWithMl || doprocessD0PiMc || doprocessD0PiMcWithMl) { + registry.add(Form("hMass%s", histMassTitle1.data()), Form("%s candidates; %s (GeV/#it{c}^{2});entries", charmHadTitle1.data(), charmHadInvMassTitle.data()), {HistType::kTH1F, {{500, 0., 5.}}}); + } + registry.add(Form("hPt%s", histMassTitle0.data()), Form("%s candidates candidates;%s candidate #it{p}_{T} (GeV/#it{c});entries", charmHadTitle0.data(), charmHadTitle0.data()), {HistType::kTH1F, {{100, 0., 10.}}}); + registry.add("hPtPion", "#pi^{#plus} candidates;#pi^{#plus} candidate #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{100, 0., 10.}}}); + registry.add(Form("hCpa%s", histMassTitle0.data()), Form("%s candidates;%s cosine of pointing angle;entries", charmHadTitle0.data(), charmHadTitle0.data()), {HistType::kTH1F, {{110, -1.1, 1.1}}}); } /// Pion selection (D Pi <-- B0) @@ -435,20 +461,22 @@ struct HfDataCreatorCharmHadPiReduced { if constexpr (decChannel == DecayChannel::B0ToDminusPi) { indexHfCandCharm = hfCand3Prong.lastIndex() + 1; invMassC0 = hfHelper.invMassDplusToPiKPi(candC); - registry.fill(HIST("hMassCharmHad"), invMassC0); + registry.fill(HIST("hMassDplus"), invMassC0); + registry.fill(HIST("hPtDplus"), candC.pt()); + registry.fill(HIST("hCpaDplus"), candC.cpa()); } else if constexpr (decChannel == DecayChannel::BplusToD0barPi){ indexHfCandCharm = hfCand2Prong.lastIndex() + 1; if (candC.isSelD0() >= selectionFlagD0) { invMassC0 = hfHelper.invMassD0ToPiK(candC); - registry.fill(HIST("hMassCharmHad0"), invMassC0); + registry.fill(HIST("hMassD0"), invMassC0); } if (candC.isSelD0bar() >= selectionFlagD0bar) { invMassC1 = hfHelper.invMassD0barToKPi(candC); - registry.fill(HIST("hMassCharmHad1"), invMassC1); + registry.fill(HIST("hMassD0bar"), invMassC1); } + registry.fill(HIST("hPtD0"), candC.pt()); + registry.fill(HIST("hCpaD0"), candC.cpa()); } - registry.fill(HIST("hPtCharmHad"), candC.pt()); - registry.fill(HIST("hCpaCharmHad"), candC.cpa()); bool fillHfCandCharm = false; diff --git a/PWGHF/D2H/TableProducer/dataCreatorDV0Reduced.cxx b/PWGHF/D2H/TableProducer/dataCreatorDV0Reduced.cxx index 8d6fd7c947c..908b9ac4e1e 100644 --- a/PWGHF/D2H/TableProducer/dataCreatorDV0Reduced.cxx +++ b/PWGHF/D2H/TableProducer/dataCreatorDV0Reduced.cxx @@ -70,7 +70,7 @@ enum TypeD : uint8_t { /// Creation of D-V0 pairs struct HfDataCreatorDV0Reduced { // Produces AOD tables to store track information - Produces hfReducedCollision; // Defined in PWGLF/DataModel/LFStrangenessTables.h + Produces hfReducedCollision; // Defined in PWGHF/D2H/DataModel/ReducedDataModel.h Produces hfCollisionCounter; // Defined in PWGHF/D2H/DataModel/ReducedDataModel.h Produces hfCandV0; // Defined in PWGHF/D2H/DataModel/ReducedDataModel.h @@ -373,10 +373,7 @@ struct HfDataCreatorDV0Reduced { } registry.fill(HIST("hEvents"), 1 + Event::DV0Selected); // fill collision table if it contains a DPi pair a minima - hfReducedCollision(collision.posX(), collision.posY(), collision.posZ(), - collision.covXX(), collision.covXY(), collision.covYY(), - collision.covXZ(), collision.covYZ(), collision.covZZ(), - 0); + hfReducedCollision(collision.posX(), collision.posY(), collision.posZ()); } // run data creation void processDplusV0(aod::Collisions const& collisions, From 4aed8a657367b9709e26cc65044f440e1d3d463a Mon Sep 17 00:00:00 2001 From: ALICE Action Bot Date: Thu, 29 Feb 2024 14:51:28 +0000 Subject: [PATCH 3/5] Please consider the following formatting changes --- .../dataCreatorCharmHadPiReduced.cxx | 44 +++++++++---------- .../TableProducer/dataCreatorDV0Reduced.cxx | 2 +- 2 files changed, 23 insertions(+), 23 deletions(-) diff --git a/PWGHF/D2H/TableProducer/dataCreatorCharmHadPiReduced.cxx b/PWGHF/D2H/TableProducer/dataCreatorCharmHadPiReduced.cxx index be02117f4fc..91e75b8ef58 100644 --- a/PWGHF/D2H/TableProducer/dataCreatorCharmHadPiReduced.cxx +++ b/PWGHF/D2H/TableProducer/dataCreatorCharmHadPiReduced.cxx @@ -291,7 +291,8 @@ struct HfDataCreatorCharmHadPiReduced { void fillMcRecoInfo(const PParticles& particlesMc, const std::vector& vecDaughtersB, int& indexHfCandCharm, - std::map selectedTracksPion) { + std::map selectedTracksPion) + { // we check the MC matching to be stored int8_t sign{0}; @@ -395,8 +396,7 @@ struct HfDataCreatorCharmHadPiReduced { } rowHfDPiMcCheckReduced(pdgCodeBeautyMother, pdgCodeCharmMother, pdgCodeProng0, pdgCodeProng1, pdgCodeProng2, pdgCodeProng3); } - } - else if constexpr (decChannel == DecayChannel::BplusToD0barPi) { + } else if constexpr (decChannel == DecayChannel::BplusToD0barPi) { // 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) { @@ -464,7 +464,7 @@ struct HfDataCreatorCharmHadPiReduced { registry.fill(HIST("hMassDplus"), invMassC0); registry.fill(HIST("hPtDplus"), candC.pt()); registry.fill(HIST("hCpaDplus"), candC.cpa()); - } else if constexpr (decChannel == DecayChannel::BplusToD0barPi){ + } else if constexpr (decChannel == DecayChannel::BplusToD0barPi) { indexHfCandCharm = hfCand2Prong.lastIndex() + 1; if (candC.isSelD0() >= selectionFlagD0) { invMassC0 = hfHelper.invMassD0ToPiK(candC); @@ -529,7 +529,7 @@ struct HfDataCreatorCharmHadPiReduced { df3.getTrack(1).getPxPyPzGlo(pVec1); df3.getTrack(2).getPxPyPzGlo(pVec2); pVecCharm = RecoDecay::pVec(pVec0, pVec1, pVec2); - trackParCovCharmHad.setAbsCharge(charmHadDauTracks[1].sign()); // to be sure + trackParCovCharmHad.setAbsCharge(charmHadDauTracks[1].sign()); // to be sure } else if constexpr (decChannel == DecayChannel::BplusToD0barPi) { // D0(bar) → K± π∓ if (df2.process(trackParCov0, trackParCov1) == 0) { continue; @@ -571,7 +571,7 @@ struct HfDataCreatorCharmHadPiReduced { if (!isPionSelected(trackPion, trackParCovPion, dcaPion, charmHadDauTracks)) { continue; } - + registry.fill(HIST("hPtPion"), trackParCovPion.getPt()); // compute invariant mass square and apply selection auto invMass2DPi = RecoDecay::m2(std::array{pVecCharm, pVecPion}, std::array{massC, massPi}); @@ -607,15 +607,15 @@ struct HfDataCreatorCharmHadPiReduced { fillMcRecoInfo(particlesMc, beautyHadDauTracks, indexHfCandCharm, selectedTracksPion); } fillHfCandCharm = true; - } // pion loop - if (fillHfCandCharm) { // fill candCplus table only once per D candidate + } // pion loop + if (fillHfCandCharm) { // fill candCplus table only once per D candidate if constexpr (decChannel == DecayChannel::B0ToDminusPi) { // D∓ → π∓ K± π∓ hfCand3Prong(charmHadDauTracks[0].globalIndex(), charmHadDauTracks[1].globalIndex(), charmHadDauTracks[2].globalIndex(), - indexHfReducedCollision, - trackParCovCharmHad.getX(), trackParCovCharmHad.getAlpha(), - trackParCovCharmHad.getY(), trackParCovCharmHad.getZ(), trackParCovCharmHad.getSnp(), - trackParCovCharmHad.getTgl(), trackParCovCharmHad.getQ2Pt(), - candC.xSecondaryVertex(), candC.ySecondaryVertex(), candC.zSecondaryVertex(), invMassC0); + indexHfReducedCollision, + trackParCovCharmHad.getX(), trackParCovCharmHad.getAlpha(), + trackParCovCharmHad.getY(), trackParCovCharmHad.getZ(), trackParCovCharmHad.getSnp(), + trackParCovCharmHad.getTgl(), trackParCovCharmHad.getQ2Pt(), + candC.xSecondaryVertex(), candC.ySecondaryVertex(), candC.zSecondaryVertex(), invMassC0); hfCand3ProngCov(trackParCovCharmHad.getSigmaY2(), trackParCovCharmHad.getSigmaZY(), trackParCovCharmHad.getSigmaZ2(), trackParCovCharmHad.getSigmaSnpY(), trackParCovCharmHad.getSigmaSnpZ(), trackParCovCharmHad.getSigmaSnp2(), trackParCovCharmHad.getSigmaTglY(), trackParCovCharmHad.getSigmaTglZ(), @@ -627,11 +627,11 @@ struct HfDataCreatorCharmHadPiReduced { } } else if constexpr (decChannel == DecayChannel::BplusToD0barPi) { // D0(bar) → K± π∓ hfCand2Prong(charmHadDauTracks[0].globalIndex(), charmHadDauTracks[1].globalIndex(), - indexHfReducedCollision, - trackParCovCharmHad.getX(), trackParCovCharmHad.getAlpha(), - trackParCovCharmHad.getY(), trackParCovCharmHad.getZ(), trackParCovCharmHad.getSnp(), - trackParCovCharmHad.getTgl(), trackParCovCharmHad.getQ2Pt(), - candC.xSecondaryVertex(), candC.ySecondaryVertex(), candC.zSecondaryVertex(), invMassC0, invMassC1); + indexHfReducedCollision, + trackParCovCharmHad.getX(), trackParCovCharmHad.getAlpha(), + trackParCovCharmHad.getY(), trackParCovCharmHad.getZ(), trackParCovCharmHad.getSnp(), + trackParCovCharmHad.getTgl(), trackParCovCharmHad.getQ2Pt(), + candC.xSecondaryVertex(), candC.ySecondaryVertex(), candC.zSecondaryVertex(), invMassC0, invMassC1); hfCand2ProngCov(trackParCovCharmHad.getSigmaY2(), trackParCovCharmHad.getSigmaZY(), trackParCovCharmHad.getSigmaZ2(), trackParCovCharmHad.getSigmaSnpY(), trackParCovCharmHad.getSigmaSnpZ(), trackParCovCharmHad.getSigmaSnp2(), trackParCovCharmHad.getSigmaTglY(), trackParCovCharmHad.getSigmaTglZ(), @@ -641,10 +641,10 @@ struct HfDataCreatorCharmHadPiReduced { if constexpr (withMl) { std::array bdtScores = {-1.f, -1.f, -1.f, -1.f, -1.f, -1.f}; if (candC.mlProbD0().size() == 3) { - std::copy(candC.mlProbD0().begin(), candC.mlProbD0().end(), bdtScores.begin()); + std::copy(candC.mlProbD0().begin(), candC.mlProbD0().end(), bdtScores.begin()); } if (candC.mlProbD0bar().size() == 3) { - std::copy(candC.mlProbD0bar().begin(), candC.mlProbD0bar().end(), bdtScores.begin() + 3); + std::copy(candC.mlProbD0bar().begin(), candC.mlProbD0bar().end(), bdtScores.begin() + 3); } hfCand2ProngMl(bdtScores[0], bdtScores[1], bdtScores[2], bdtScores[3], bdtScores[4], bdtScores[5]); @@ -667,7 +667,7 @@ struct HfDataCreatorCharmHadPiReduced { bz); } - template + template void runMcGen(aod::McParticles const& particlesMc) { // Match generated particles. @@ -708,7 +708,7 @@ struct HfDataCreatorCharmHadPiReduced { ptProngs[0], yProngs[0], etaProngs[0], ptProngs[1], yProngs[1], etaProngs[1]); } else if constexpr (decayChannel == DecayChannel::BplusToD0barPi) { - // B+ → D0bar π+ + // B+ → D0bar π+ if (RecoDecay::isMatchedMCGen(particlesMc, particle, Pdg::kBPlus, std::array{static_cast(Pdg::kD0), +kPiPlus}, true)) { // Match D0bar -> π- K+ auto candD0MC = particlesMc.rawIteratorAt(particle.daughtersIds().front()); diff --git a/PWGHF/D2H/TableProducer/dataCreatorDV0Reduced.cxx b/PWGHF/D2H/TableProducer/dataCreatorDV0Reduced.cxx index 908b9ac4e1e..cc89b6beb5d 100644 --- a/PWGHF/D2H/TableProducer/dataCreatorDV0Reduced.cxx +++ b/PWGHF/D2H/TableProducer/dataCreatorDV0Reduced.cxx @@ -70,7 +70,7 @@ enum TypeD : uint8_t { /// Creation of D-V0 pairs struct HfDataCreatorDV0Reduced { // Produces AOD tables to store track information - Produces hfReducedCollision; // Defined in PWGHF/D2H/DataModel/ReducedDataModel.h + Produces hfReducedCollision; // Defined in PWGHF/D2H/DataModel/ReducedDataModel.h Produces hfCollisionCounter; // Defined in PWGHF/D2H/DataModel/ReducedDataModel.h Produces hfCandV0; // Defined in PWGHF/D2H/DataModel/ReducedDataModel.h From 014246c52c7b6bbf602986707804f71e865d5b23 Mon Sep 17 00:00:00 2001 From: Fabrizio Grosa Date: Thu, 29 Feb 2024 17:12:21 +0100 Subject: [PATCH 4/5] Fix typo --- PWGHF/D2H/TableProducer/dataCreatorCharmHadPiReduced.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/PWGHF/D2H/TableProducer/dataCreatorCharmHadPiReduced.cxx b/PWGHF/D2H/TableProducer/dataCreatorCharmHadPiReduced.cxx index 91e75b8ef58..7a7a90fdce0 100644 --- a/PWGHF/D2H/TableProducer/dataCreatorCharmHadPiReduced.cxx +++ b/PWGHF/D2H/TableProducer/dataCreatorCharmHadPiReduced.cxx @@ -153,7 +153,7 @@ struct HfDataCreatorCharmHadPiReduced { void init(InitContext const&) { std::array doProcess = {doprocessDplusPiData, doprocessDplusPiDataWithMl, doprocessDplusPiMc, doprocessDplusPiMcWithMl, doprocessD0PiData, doprocessD0PiDataWithMl, doprocessD0PiMc, doprocessD0PiMcWithMl}; - if (std::accumulate(doProcess.begin(), doProcess.end(), 0) != 0) { + if (std::accumulate(doProcess.begin(), doProcess.end(), 0) != 1) { LOGP(fatal, "One and only one process function can be enabled at a time, please fix your configuration!"); } @@ -817,7 +817,7 @@ struct HfDataCreatorCharmHadPiReduced { runDataCreation(collision, candsCThisColl, trackIdsThisCollision, tracks, tracks, bcs); } } - PROCESS_SWITCH(HfDataCreatorCharmHadPiReduced, processD0PiData, "Process D0Pi without MC info and without ML info", true); + PROCESS_SWITCH(HfDataCreatorCharmHadPiReduced, processD0PiData, "Process D0Pi without MC info and without ML info", false); void processD0PiDataWithMl(aod::Collisions const& collisions, CandsD0FilteredWithMl const& candsC, From 9cdb110b598511ffadcfce2e462d8f33939bd435 Mon Sep 17 00:00:00 2001 From: Fabrizio Grosa Date: Thu, 29 Feb 2024 18:03:21 +0100 Subject: [PATCH 5/5] Fix std vector copy --- PWGHF/D2H/TableProducer/dataCreatorCharmHadPiReduced.cxx | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/PWGHF/D2H/TableProducer/dataCreatorCharmHadPiReduced.cxx b/PWGHF/D2H/TableProducer/dataCreatorCharmHadPiReduced.cxx index 7a7a90fdce0..70d66c2b68e 100644 --- a/PWGHF/D2H/TableProducer/dataCreatorCharmHadPiReduced.cxx +++ b/PWGHF/D2H/TableProducer/dataCreatorCharmHadPiReduced.cxx @@ -602,7 +602,9 @@ struct HfDataCreatorCharmHadPiReduced { if constexpr (doMc) { std::vector beautyHadDauTracks{}; - std::copy(charmHadDauTracks.begin(), charmHadDauTracks.end(), beautyHadDauTracks.begin()); + for (const auto& track : charmHadDauTracks) { + beautyHadDauTracks.push_back(track); + } beautyHadDauTracks.push_back(trackPion); fillMcRecoInfo(particlesMc, beautyHadDauTracks, indexHfCandCharm, selectedTracksPion); }