From e88eb897b3ad3154cb2815dbf19a2b46e95e2a57 Mon Sep 17 00:00:00 2001 From: samrangy Date: Fri, 13 Oct 2023 12:14:01 +0530 Subject: [PATCH 01/11] type conversion for PDG value --- PWGHF/HFC/TableProducer/correlatorD0Hadrons.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PWGHF/HFC/TableProducer/correlatorD0Hadrons.cxx b/PWGHF/HFC/TableProducer/correlatorD0Hadrons.cxx index 8e21d248eff..1beac83f158 100644 --- a/PWGHF/HFC/TableProducer/correlatorD0Hadrons.cxx +++ b/PWGHF/HFC/TableProducer/correlatorD0Hadrons.cxx @@ -199,7 +199,7 @@ struct HfCorrelatorD0Hadrons { (aod::track::dcaZ > static_cast(-dcaZTrackMax)) && (aod::track::dcaZ < static_cast(dcaZTrackMax)); Filter d0Filter = (aod::hf_sel_candidate_d0::isSelD0 >= 1) || (aod::hf_sel_candidate_d0::isSelD0bar >= 1); Filter collisionFilterGen = aod::hf_selection_dmeson_collision::dmesonSel == true; - Filter particlesFilter = nabs(aod::mcparticle::pdgCode) == pdg::Code::kD0 || ((aod::mcparticle::flags & (uint8_t)o2::aod::mcparticle::enums::PhysicalPrimary) == (uint8_t)o2::aod::mcparticle::enums::PhysicalPrimary); + Filter particlesFilter = nabs(aod::mcparticle::pdgCode) == static_cast(pdg::Code::kD0) || ((aod::mcparticle::flags & (uint8_t)o2::aod::mcparticle::enums::PhysicalPrimary) == (uint8_t)o2::aod::mcparticle::enums::PhysicalPrimary); HistogramRegistry registry{ "registry", From 5a65fdfced99a5c319d892fccde2786b06a0b935 Mon Sep 17 00:00:00 2001 From: samrangy Date: Sat, 20 Jan 2024 13:28:20 +0530 Subject: [PATCH 02/11] extending mass axis range and modifying SB values --- PWGHF/HFC/TableProducer/correlatorD0Hadrons.cxx | 6 +++--- PWGHF/HFC/Tasks/taskCorrelationD0Hadrons.cxx | 16 ++++++++-------- 2 files changed, 11 insertions(+), 11 deletions(-) diff --git a/PWGHF/HFC/TableProducer/correlatorD0Hadrons.cxx b/PWGHF/HFC/TableProducer/correlatorD0Hadrons.cxx index 7a2aa274756..681ef3ef042 100644 --- a/PWGHF/HFC/TableProducer/correlatorD0Hadrons.cxx +++ b/PWGHF/HFC/TableProducer/correlatorD0Hadrons.cxx @@ -49,9 +49,9 @@ const double efficiencyDmesonDefault[nPtBinsMassAndEfficiency] = {}; auto vecEfficiencyDmeson = std::vector{efficiencyDmesonDefault, efficiencyDmesonDefault + nPtBinsMassAndEfficiency}; // histogram binning definition -const int massAxisNBins = 120; -const double massAxisMin = 1.5848; -const double massAxisMax = 2.1848; +const int massAxisNBins = 200; +const double massAxisMin = 1.3848; +const double massAxisMax = 2.3848; const int phiAxisNBins = 32; const double phiAxisMin = 0.; const double phiAxisMax = o2::constants::math::TwoPI; diff --git a/PWGHF/HFC/Tasks/taskCorrelationD0Hadrons.cxx b/PWGHF/HFC/Tasks/taskCorrelationD0Hadrons.cxx index 8b23113b392..1a47f61d8fc 100644 --- a/PWGHF/HFC/Tasks/taskCorrelationD0Hadrons.cxx +++ b/PWGHF/HFC/Tasks/taskCorrelationD0Hadrons.cxx @@ -54,15 +54,15 @@ AxisSpec axisPtHadron = {11, 0., 11., ""}; AxisSpec axisPoolBin = {9, 0., 9., ""}; // definition of vectors for standard ptbin and invariant mass configurables -const int nPtBinsCorrelations = 8; -const double pTBinsCorrelations[nPtBinsCorrelations + 1] = {0., 2., 4., 6., 8., 12., 16., 24., 99.}; +const int nPtBinsCorrelations = 12; +const double pTBinsCorrelations[nPtBinsCorrelations + 1] = {0., 1., 2., 3., 4., 5., 6., 7., 8., 12., 16., 24., 99.}; auto vecPtBinsCorrelations = std::vector{pTBinsCorrelations, pTBinsCorrelations + nPtBinsCorrelations + 1}; -const double signalRegionLeftDefault[nPtBinsCorrelations] = {1.810, 1.810, 1.810, 1.810, 1.810, 1.810, 1.810, 1.810}; -const double signalRegionRightDefault[nPtBinsCorrelations] = {1.922, 1.922, 1.922, 1.922, 1.922, 1.922, 1.922, 1.922}; -const double sidebandLeftInnerDefault[nPtBinsCorrelations] = {1.754, 1.754, 1.754, 1.754, 1.754, 1.754, 1.754, 1.754}; -const double sidebandLeftOuterDefault[nPtBinsCorrelations] = {1.642, 1.642, 1.642, 1.642, 1.642, 1.642, 1.642, 1.642}; -const double sidebandRightInnerDefault[nPtBinsCorrelations] = {1.978, 1.978, 1.978, 1.978, 1.978, 1.978, 1.978, 1.978}; -const double sidebandRightOuterDefault[nPtBinsCorrelations] = {2.090, 2.090, 2.090, 2.090, 2.090, 2.090, 2.090, 2.090}; +const double signalRegionLeftDefault[nPtBinsCorrelations] = {1.7948, 1.8198, 1.8198, 1.8148, 1.8148, 1.8048, 1.8048, 1.7948, 1.7948, 1.7898, 1.7848, 1.7598}; +const double signalRegionRightDefault[nPtBinsCorrelations] = {1.9098, 1.8998, 1.9048, 1.9048, 1.9148, 1.9248, 1.9298, 1.9348, 1.9398, 1.9298, 1.9398, 1.9198}; +const double sidebandLeftInnerDefault[nPtBinsCorrelations] = {1.7398, 1.7748, 1.7798, 1.7698, 1.7648, 1.7448, 1.7448, 1.7198, 1.7198, 1.7198, 1.7048, 1.6798}; +const double sidebandLeftOuterDefault[nPtBinsCorrelations] = {1.6298, 1.6898, 1.6948, 1.6748, 1.6648, 1.6248, 1.6198, 1.5748, 1.5748, 1.5798, 1.5448, 1.5198}; +const double sidebandRightInnerDefault[nPtBinsCorrelations] = {1.9648, 1.9448, 1.9448, 1.9548, 1.9648, 1.9848, 1.9948, 2.0098, 2.0148, 1.9998, 2.0248, 1.9998}; +const double sidebandRightOuterDefault[nPtBinsCorrelations] = {2.0748, 2.0248, 2.0298, 2.0448, 2.0648, 2.1048, 2.1148, 2.1548, 2.1648, 2.1398, 2.1848, 2.1598}; auto vecsignalRegionLeft = std::vector{signalRegionLeftDefault, signalRegionLeftDefault + nPtBinsCorrelations}; auto vecsignalRegionRight = std::vector{signalRegionRightDefault, signalRegionRightDefault + nPtBinsCorrelations}; auto vecSidebandLeftInner = std::vector{sidebandLeftInnerDefault, sidebandLeftInnerDefault + nPtBinsCorrelations}; From dc85814565e69bd67f789ad57008335f3b46e741 Mon Sep 17 00:00:00 2001 From: samrangy Date: Sat, 20 Jan 2024 14:40:06 +0530 Subject: [PATCH 03/11] resolving conflicts --- PWGHF/HFC/TableProducer/correlatorD0Hadrons.cxx | 1 - 1 file changed, 1 deletion(-) diff --git a/PWGHF/HFC/TableProducer/correlatorD0Hadrons.cxx b/PWGHF/HFC/TableProducer/correlatorD0Hadrons.cxx index 5c18e92b1ac..681ef3ef042 100644 --- a/PWGHF/HFC/TableProducer/correlatorD0Hadrons.cxx +++ b/PWGHF/HFC/TableProducer/correlatorD0Hadrons.cxx @@ -204,7 +204,6 @@ struct HfCorrelatorD0Hadrons { Filter d0Filter = (aod::hf_sel_candidate_d0::isSelD0 >= 1) || (aod::hf_sel_candidate_d0::isSelD0bar >= 1); Filter collisionFilterGen = aod::hf_selection_dmeson_collision::dmesonSel == true; Filter particlesFilter = nabs(aod::mcparticle::pdgCode) == static_cast(Pdg::kD0) || ((aod::mcparticle::flags & (uint8_t)o2::aod::mcparticle::enums::PhysicalPrimary) == (uint8_t)o2::aod::mcparticle::enums::PhysicalPrimary); - Filter particlesFilter = nabs(aod::mcparticle::pdgCode) == static_cast(pdg::Code::kD0) || ((aod::mcparticle::flags & (uint8_t)o2::aod::mcparticle::enums::PhysicalPrimary) == (uint8_t)o2::aod::mcparticle::enums::PhysicalPrimary); HistogramRegistry registry{ "registry", From 7127a94b568ffaa3d3a532da2f25a5e25559d7ff Mon Sep 17 00:00:00 2001 From: ALICE Action Bot Date: Sat, 20 Jan 2024 09:17:34 +0000 Subject: [PATCH 04/11] MegaLinter fixes --- PWGHF/HFC/TableProducer/correlatorD0Hadrons.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/PWGHF/HFC/TableProducer/correlatorD0Hadrons.cxx b/PWGHF/HFC/TableProducer/correlatorD0Hadrons.cxx index 681ef3ef042..da890154128 100644 --- a/PWGHF/HFC/TableProducer/correlatorD0Hadrons.cxx +++ b/PWGHF/HFC/TableProducer/correlatorD0Hadrons.cxx @@ -458,7 +458,7 @@ struct HfCorrelatorD0Hadrons { registry.fill(HIST("hSelectionStatusRec"), candidate1.isSelD0bar() + (candidate1.isSelD0() * 2)); } // fill invariant mass plots from D0/D0bar signal and background candidates - if (candidate1.isSelD0() >= selectionFlagD0) { // only reco as D0 + if (candidate1.isSelD0() >= selectionFlagD0) { // only reco as D0 if (candidate1.flagMcMatchRec() == 1 << aod::hf_cand_2prong::DecayType::D0ToPiK) { // also matched as D0 registry.fill(HIST("hMassD0RecSig"), hfHelper.invMassD0ToPiK(candidate1), candidate1.pt(), efficiencyWeight); } else if (candidate1.flagMcMatchRec() == -(1 << aod::hf_cand_2prong::DecayType::D0ToPiK)) { @@ -467,7 +467,7 @@ struct HfCorrelatorD0Hadrons { registry.fill(HIST("hMassD0RecBg"), hfHelper.invMassD0ToPiK(candidate1), candidate1.pt(), efficiencyWeight); } } - if (candidate1.isSelD0bar() >= selectionFlagD0bar) { // only reco as D0bar + if (candidate1.isSelD0bar() >= selectionFlagD0bar) { // only reco as D0bar if (candidate1.flagMcMatchRec() == -(1 << aod::hf_cand_2prong::DecayType::D0ToPiK)) { // also matched as D0bar registry.fill(HIST("hMassD0barRecSig"), hfHelper.invMassD0barToKPi(candidate1), candidate1.pt(), efficiencyWeight); } else if (candidate1.flagMcMatchRec() == 1 << aod::hf_cand_2prong::DecayType::D0ToPiK) { From 8948e0c0d461fef9ca6b4fdd4bcecc98a1a656f2 Mon Sep 17 00:00:00 2001 From: ALICE Action Bot Date: Sat, 15 Jun 2024 19:01:18 +0000 Subject: [PATCH 05/11] Please consider the following formatting changes --- .../HFC/TableProducer/correlatorD0Hadrons.cxx | 917 ------------------ 1 file changed, 917 deletions(-) delete mode 100644 PWGHF/HFC/TableProducer/correlatorD0Hadrons.cxx diff --git a/PWGHF/HFC/TableProducer/correlatorD0Hadrons.cxx b/PWGHF/HFC/TableProducer/correlatorD0Hadrons.cxx deleted file mode 100644 index 9ec7e5fc7c2..00000000000 --- a/PWGHF/HFC/TableProducer/correlatorD0Hadrons.cxx +++ /dev/null @@ -1,917 +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 correlatorD0Hadrons.cxx -/// \brief D0-Hadron correlator task - data-like, MC-reco and MC-kine analyses. -/// -/// \author Samrangy Sadhu , INFN Bari -/// \author Swapnesh Santosh Khade , IIT Indore - -#include - -#include "CommonConstants/PhysicsConstants.h" -#include "Framework/AnalysisTask.h" -#include "Framework/HistogramRegistry.h" -#include "Framework/runDataProcessing.h" - -#include "Common/Core/TrackSelection.h" -#include "Common/DataModel/Centrality.h" -#include "Common/DataModel/EventSelection.h" -#include "Common/DataModel/Multiplicity.h" -#include "Common/DataModel/TrackSelectionTables.h" - -#include "PWGHF/Core/HfHelper.h" -#include "PWGHF/DataModel/CandidateReconstructionTables.h" -#include "PWGHF/DataModel/CandidateSelectionTables.h" -#include "PWGHF/HFC/DataModel/CorrelationTables.h" -#include "PWGHF/HFC/Utils/utilsCorrelations.h" - -using namespace o2; -using namespace o2::analysis; -using namespace o2::constants::physics; -using namespace o2::framework; -using namespace o2::framework::expressions; -using namespace o2::analysis::hf_correlations; - -/// -/// Returns deltaPhi value in range [-pi/2., 3.*pi/2], typically used for correlation studies -/// -double getDeltaPhi(double phiHadron, double phiD) -{ - return RecoDecay::constrainAngle(phiHadron - phiD, -o2::constants::math::PIHalf); -} - -const int nPtBinsMassAndEfficiency = o2::analysis::hf_cuts_d0_to_pi_k::nBinsPt; -const double efficiencyDmesonDefault[nPtBinsMassAndEfficiency] = {}; -auto vecEfficiencyDmeson = std::vector{efficiencyDmesonDefault, efficiencyDmesonDefault + nPtBinsMassAndEfficiency}; - -// histogram binning definition -const int massAxisNBins = 200; -const double massAxisMin = 1.3848; -const double massAxisMax = 2.3848; -const int phiAxisNBins = 32; -const double phiAxisMin = 0.; -const double phiAxisMax = o2::constants::math::TwoPI; -const int yAxisNBins = 100; -const double yAxisMin = -5.; -const double yAxisMax = 5.; -const int ptDAxisNBins = 180; -const double ptDAxisMin = 0.; -const double ptDAxisMax = 36.; - -// Types -using BinningType = ColumnBinningPolicy>; -using SelectedCollisions = soa::Filtered>; -using SelectedTracks = soa::Filtered>; -using SelectedCandidatesData = soa::Filtered>; -using SelectedTracksMcRec = soa::Filtered>; -using SelectedCandidatesMcRec = soa::Filtered>; -using SelectedCollisionsMcGen = soa::Filtered>; -using SelectedParticlesMcGen = soa::Join; - -// Code to select collisions with at least one D0 -struct HfCorrelatorD0HadronsSelection { - SliceCache cache; - - Produces d0Sel; - - Configurable selectionFlagD0{"selectionFlagD0", 1, "Selection Flag for D0"}; - Configurable selectionFlagD0bar{"selectionFlagD0bar", 1, "Selection Flag for D0bar"}; - Configurable yCandMax{"yCandMax", 4.0, "max. cand. rapidity"}; - Configurable ptCandMin{"ptCandMin", -1., "min. cand. pT"}; - - HfHelper hfHelper; - - Preslice perCol = aod::hf_cand::collisionId; - Partition> selectedD0Candidates = aod::hf_sel_candidate_d0::isSelD0 >= selectionFlagD0 || aod::hf_sel_candidate_d0::isSelD0bar >= selectionFlagD0bar; - Partition> selectedD0candidatesMc = aod::hf_sel_candidate_d0::isSelD0 >= selectionFlagD0 || aod::hf_sel_candidate_d0::isSelD0bar >= selectionFlagD0bar; - - void processD0SelectionData(aod::Collision const& collision, - soa::Join const&) - { - bool isD0Found = 0; - if (selectedD0Candidates.size() > 0) { - auto selectedD0CandidatesGrouped = selectedD0Candidates->sliceByCached(aod::hf_cand::collisionId, collision.globalIndex(), cache); - - for (const auto& candidate : selectedD0CandidatesGrouped) { - // check decay channel flag for candidate - if (!TESTBIT(candidate.hfflag(), aod::hf_cand_2prong::DecayType::D0ToPiK)) { - continue; - } - if (std::abs(hfHelper.yD0(candidate)) > yCandMax || candidate.pt() < ptCandMin) { - continue; - } - isD0Found = 1; - break; - } - } - d0Sel(isD0Found); - } - PROCESS_SWITCH(HfCorrelatorD0HadronsSelection, processD0SelectionData, "Process D0 Selection Data", false); - - void processD0SelectionMcRec(aod::Collision const& collision, - soa::Join const&) - { - bool isD0Found = 0; - if (selectedD0candidatesMc.size() > 0) { - auto selectedD0CandidatesGroupedMc = selectedD0candidatesMc->sliceByCached(aod::hf_cand::collisionId, collision.globalIndex(), cache); - for (const auto& candidate : selectedD0CandidatesGroupedMc) { - // check decay channel flag for candidate - if (!TESTBIT(candidate.hfflag(), aod::hf_cand_2prong::DecayType::D0ToPiK)) { - continue; - } - if (std::abs(hfHelper.yD0(candidate)) > yCandMax || candidate.pt() < ptCandMin) { - continue; - } - isD0Found = 1; - break; - } - } - d0Sel(isD0Found); - } - PROCESS_SWITCH(HfCorrelatorD0HadronsSelection, processD0SelectionMcRec, "Process D0 Selection MCRec", true); - - void processD0SelectionMcGen(aod::McCollision const&, - aod::McParticles const& mcParticles) - { - bool isD0Found = 0; - for (const auto& particle : mcParticles) { - if (std::abs(particle.pdgCode()) != Pdg::kD0) { - continue; - } - double yD = RecoDecay::y(particle.pVector(), MassD0); - if (std::abs(yD) > yCandMax || particle.pt() < ptCandMin) { - continue; - } - isD0Found = 1; - break; - } - d0Sel(isD0Found); - } - PROCESS_SWITCH(HfCorrelatorD0HadronsSelection, processD0SelectionMcGen, "Process D0 Selection MCGen", false); -}; - -/// D0-Hadron correlation pair builder - for real data and data-like analysis (i.e. reco-level w/o matching request via MC truth) -struct HfCorrelatorD0Hadrons { - SliceCache cache; - - Produces entryD0HadronPair; - Produces entryD0HadronRecoInfo; - - Configurable selectionFlagD0{"selectionFlagD0", 1, "Selection Flag for D0"}; - Configurable selectionFlagD0bar{"selectionFlagD0bar", 1, "Selection Flag for D0bar"}; - Configurable yCandMax{"yCandMax", 0.8, "max. cand. rapidity"}; - Configurable etaTrackMax{"etaTrackMax", 0.8, "max. eta of tracks"}; - Configurable dcaXYTrackMax{"dcaXYTrackMax", 1., "max. DCAxy of tracks"}; - Configurable dcaZTrackMax{"dcaZTrackMax", 1., "max. DCAz of tracks"}; - Configurable ptCandMin{"ptCandMin", 1., "min. cand. pT"}; - Configurable ptTrackMin{"ptTrackMin", 0.3, "min. track pT"}; - Configurable ptTrackMax{"ptTrackMax", 99., "max. track pT"}; - Configurable> bins{"ptBinsForMassAndEfficiency", std::vector{o2::analysis::hf_cuts_d0_to_pi_k::vecBinsPt}, "pT bin limits for candidate mass plots and efficiency"}; - Configurable> efficiencyDmeson{"efficiencyDmeson", std::vector{vecEfficiencyDmeson}, "Efficiency values for D0 meson"}; - Configurable applyEfficiency{"efficiencyFlagD", 1, "Flag for applying D-meson efficiency weights"}; - Configurable multMin{"multMin", 0., "minimum multiplicity accepted"}; - Configurable multMax{"multMax", 10000., "maximum multiplicity accepted"}; - Configurable ptSoftPionMax{"ptSoftPionMax", 3.f * 800.f * std::pow(10.f, -6.f), "max. pT cut for soft pion identification"}; - Configurable correlateD0WithLeadingParticle{"correlateD0WithLeadingParticle", false, "Switch for correlation of D0 mesons with leading particle only"}; - Configurable storeAutoCorrelationFlag{"storeAutoCorrelationFlag", false, "Store flag that indicates if the track is paired to its D-meson mother instead of skipping it"}; - Configurable numberEventsMixed{"numberEventsMixed", 5, "Number of events mixed in ME process"}; - ConfigurableAxis zPoolBins{"zPoolBins", {VARIABLE_WIDTH, -10.0f, -2.5f, 2.5f, 10.0f}, "z vertex position pools"}; - ConfigurableAxis multPoolBins{"multPoolBins", {VARIABLE_WIDTH, 0.0f, 2000.0f, 6000.0f, 10000.0f}, "event multiplicity pools (FT0M)"}; - ConfigurableAxis multPoolBinsMcGen{"multPoolBinsMcGen", {VARIABLE_WIDTH, 0.0f, 20.0f, 50.0f, 500.0f}, "Mixing bins - MC multiplicity"}; // In MCGen multiplicity is defined by counting tracks - - HfHelper hfHelper; - BinningType corrBinning{{zPoolBins, multPoolBins}, true}; - - int leadingIndex = 0; - double massD0{0.}; - double massPi{0.}; - double massK{0.}; - double softPiMass = 0.14543; // pion mass + Q-value of the D*->D0pi decay - - Preslice perCol = aod::hf_cand::collisionId; - - Partition> selectedD0Candidates = aod::hf_sel_candidate_d0::isSelD0 >= selectionFlagD0 || aod::hf_sel_candidate_d0::isSelD0bar >= selectionFlagD0bar; - Partition> selectedD0candidatesMc = aod::hf_sel_candidate_d0::isSelD0 >= selectionFlagD0 || aod::hf_sel_candidate_d0::isSelD0bar >= selectionFlagD0bar; - - Filter collisionFilter = aod::hf_selection_dmeson_collision::dmesonSel == true; - Filter trackFilter = requireGlobalTrackWoDCAInFilter() && (nabs(aod::track::eta) < etaTrackMax) && (aod::track::pt > ptTrackMin) && (aod::track::pt < ptTrackMax) && (nabs(aod::track::dcaXY) < dcaXYTrackMax) && (nabs(aod::track::dcaZ) < dcaZTrackMax); - Filter d0Filter = (aod::hf_sel_candidate_d0::isSelD0 >= 1) || (aod::hf_sel_candidate_d0::isSelD0bar >= 1); - Filter collisionFilterGen = aod::hf_selection_dmeson_collision::dmesonSel == true; - Filter particlesFilter = nabs(aod::mcparticle::pdgCode) == static_cast(Pdg::kD0) || ((aod::mcparticle::flags & (uint8_t)o2::aod::mcparticle::enums::PhysicalPrimary) == (uint8_t)o2::aod::mcparticle::enums::PhysicalPrimary); - - HistogramRegistry registry{ - "registry", - // NOTE: use hMassD0 for trigger normalisation (S*0.955), and hMass2DCorrelationPairs (in final task) for 2D-sideband-subtraction purposes - {{"hPtCand", "D0,D0bar candidates;candidate #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{ptDAxisNBins, ptDAxisMin, ptDAxisMax}}}}, - {"hPtProng0", "D0,D0bar candidates;prong 0 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{ptDAxisNBins, ptDAxisMin, ptDAxisMax}}}}, - {"hPtProng1", "D0,D0bar candidates;prong 1 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{ptDAxisNBins, ptDAxisMin, ptDAxisMax}}}}, - {"hSelectionStatus", "D0,D0bar candidates;selection status;entries", {HistType::kTH1F, {{4, -0.5, 3.5}}}}, - {"hEta", "D0,D0bar candidates;candidate #it{#eta};entries", {HistType::kTH1F, {{yAxisNBins, yAxisMin, yAxisMax}}}}, - {"hPhi", "D0,D0bar candidates;candidate #it{#varphi};entries", {HistType::kTH1F, {{phiAxisNBins, phiAxisMin, phiAxisMax}}}}, - {"hY", "D0,D0bar candidates;candidate #it{y};entries", {HistType::kTH1F, {{yAxisNBins, yAxisMin, yAxisMax}}}}, - {"hMultiplicityPreSelection", "multiplicity prior to selection;multiplicity;entries", {HistType::kTH1F, {{10000, 0., 10000.}}}}, - {"hMultiplicity", "multiplicity;multiplicity;entries", {HistType::kTH1F, {{10000, 0., 10000.}}}}, - {"hPtCandRec", "D0,D0bar candidates - MC reco;candidate #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{ptDAxisNBins, ptDAxisMin, ptDAxisMax}}}}, - {"hPtProng0Rec", "D0,D0bar candidates - MC reco;prong 0 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{ptDAxisNBins, ptDAxisMin, ptDAxisMax}}}}, - {"hPtProng1Rec", "D0,D0bar candidates - MC reco;prong 1 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{ptDAxisNBins, ptDAxisMin, ptDAxisMax}}}}, - {"hSelectionStatusRec", "D0,D0bar candidates - MC reco;selection status;entries", {HistType::kTH1F, {{4, -0.5, 3.5}}}}, - {"hSignalStatusMERec", "Signal Status - MC reco ME;candidate sidnalStatus;entries", {HistType::kTH1F, {{200, 0, 200}}}}, - {"hEtaRec", "D0,D0bar candidates - MC reco;candidate #it{#eta};entries", {HistType::kTH1F, {{yAxisNBins, yAxisMin, yAxisMax}}}}, - {"hPhiRec", "D0,D0bar candidates - MC reco;candidate #it{#varphi};entries", {HistType::kTH1F, {{phiAxisNBins, phiAxisMin, phiAxisMax}}}}, - {"hYRec", "D0,D0bar candidates - MC reco;candidate #it{y};entries", {HistType::kTH1F, {{yAxisNBins, yAxisMin, yAxisMax}}}}, - {"hEvtCountGen", "Event counter - MC gen;;entries", {HistType::kTH1F, {{1, -0.5, 0.5}}}}, - {"hPtCandGen", "D0,D0bar particles - MC gen;particle #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{ptDAxisNBins, ptDAxisMin, ptDAxisMax}}}}, - {"hEtaGen", "D0,D0bar particles - MC gen;particle #it{#eta};entries", {HistType::kTH1F, {{yAxisNBins, yAxisMin, yAxisMax}}}}, - {"hPhiGen", "D0,D0bar particles - MC gen;particle #it{#varphi};entries", {HistType::kTH1F, {{phiAxisNBins, phiAxisMin, phiAxisMax}}}}, - {"hYGen", "D0,D0bar candidates - MC gen;candidate #it{y};entries", {HistType::kTH1F, {{yAxisNBins, yAxisMin, yAxisMax}}}}, - {"hTrackCounter", "soft pion counter - Data", {HistType::kTH1F, {{5, 0., 5.}}}}, - {"hTrackCounterRec", "soft pion counter - MC rec", {HistType::kTH1F, {{5, 0., 5.}}}}, - {"hTrackCounterGen", "soft pion counter - MC gen", {HistType::kTH1F, {{5, 0., 5.}}}}, - {"hMultV0M", "multiplicity;multiplicity;entries", {HistType::kTH1F, {{10000, 0., 10000.}}}}, - {"hZvtx", "z vertex;z vertex;entries", {HistType::kTH1F, {{200, -20., 20.}}}}, - {"hMultFT0M", "Multiplicity FT0M", {HistType::kTH1F, {{10000, 0., 10000.}}}}, - {"hD0Bin", "D0 selected in pool Bin;pool Bin;entries", {HistType::kTH1F, {{9, 0., 9.}}}}, - {"hTracksBin", "Tracks selected in pool Bin;pool Bin;entries", {HistType::kTH1F, {{9, 0., 9.}}}}}}; - - void init(InitContext&) - { - massD0 = MassD0; - massPi = MassPiPlus; - massK = MassKPlus; - - auto vbins = (std::vector)bins; - registry.add("hMass", "D0,D0bar candidates;inv. mass (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{massAxisNBins, massAxisMin, massAxisMax}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); - registry.add("hMass1D", "D0,D0bar candidates;inv. mass (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH1F, {{massAxisNBins, massAxisMin, massAxisMax}}}); - registry.add("hMassD01D", "D0,D0bar candidates;inv. mass (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH1F, {{massAxisNBins, massAxisMin, massAxisMax}}}); - registry.add("hMassD0bar1D", "D0,D0bar candidates;inv. mass (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH1F, {{massAxisNBins, massAxisMin, massAxisMax}}}); - // mass histogram for D0 signal candidates - registry.add("hMassD0RecSig", "D0 signal candidates - MC reco;inv. mass (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{massAxisNBins, massAxisMin, massAxisMax}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); - // mass histogram for D0 Reflection candidates - registry.add("hMassD0RecRef", "D0 reflection candidates - MC reco;inv. mass (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{massAxisNBins, massAxisMin, massAxisMax}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); - // mass histogram for D0 background candidates - registry.add("hMassD0RecBg", "D0 background candidates - MC reco;inv. mass (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{massAxisNBins, massAxisMin, massAxisMax}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); - // mass histogram for D0bar signal candidates - registry.add("hMassD0barRecSig", "D0bar signal candidates - MC reco;inv. mass D0bar only (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{massAxisNBins, massAxisMin, massAxisMax}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); - // mass histogram for D0bar Reflection candidates - registry.add("hMassD0barRecRef", "D0bar reflection candidates - MC reco;inv. mass D0bar only (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{massAxisNBins, massAxisMin, massAxisMax}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); - // mass histogram for D0bar background candidates - registry.add("hMassD0barRecBg", "D0bar background candidates - MC reco;inv. mass D0bar only (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{massAxisNBins, massAxisMin, massAxisMax}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); - registry.add("hCountD0TriggersGen", "D0 trigger particles - MC gen;;N of trigger D0", {HistType::kTH2F, {{1, -0.5, 0.5}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); - } - - // ======= Process starts for Data, Same event ============ - - /// D0-h correlation pair builder - for real data and data-like analysis (i.e. reco-level w/o matching request via MC truth) - void processData(SelectedCollisions::iterator const& collision, - SelectedTracks const& tracks, - soa::Join const&) - { - // protection against empty tables to be sliced - if (selectedD0Candidates.size() == 0) { - return; - } - // find leading particle - if (correlateD0WithLeadingParticle) { - leadingIndex = findLeadingParticle(tracks, dcaXYTrackMax.value, dcaZTrackMax.value, etaTrackMax.value); - } - - int poolBin = corrBinning.getBin(std::make_tuple(collision.posZ(), collision.multFT0M())); - int nTracks = 0; - if (collision.numContrib() > 1) { - for (const auto& track : tracks) { - if (std::abs(track.eta()) > etaTrackMax) { - continue; - } - if (std::abs(track.dcaXY()) > dcaXYTrackMax || std::abs(track.dcaZ()) > dcaZTrackMax) { - continue; - } - nTracks++; - } - } - registry.fill(HIST("hMultiplicityPreSelection"), nTracks); - registry.fill(HIST("hZvtx"), collision.posZ()); - registry.fill(HIST("hMultFT0M"), collision.multFT0M()); - if (nTracks < multMin || nTracks > multMax) { - return; - } - registry.fill(HIST("hMultiplicity"), nTracks); - - auto selectedD0CandidatesGrouped = selectedD0Candidates->sliceByCached(aod::hf_cand::collisionId, collision.globalIndex(), cache); - - for (const auto& candidate1 : selectedD0CandidatesGrouped) { - if (std::abs(hfHelper.yD0(candidate1)) >= yCandMax || candidate1.pt() <= ptCandMin || candidate1.pt() >= ptTrackMax) { - continue; - } - // check decay channel flag for candidate1 - if (!TESTBIT(candidate1.hfflag(), aod::hf_cand_2prong::DecayType::D0ToPiK)) { - continue; - } - - // ========================== Define parameters for soft pion removal ================================ - auto ePiK = RecoDecay::e(candidate1.pVectorProng0(), massPi) + RecoDecay::e(candidate1.pVectorProng1(), massK); - auto eKPi = RecoDecay::e(candidate1.pVectorProng0(), massK) + RecoDecay::e(candidate1.pVectorProng1(), massPi); - - // ========================== trigger efficiency ================================ - double efficiencyWeight = 1.; - if (applyEfficiency) { - efficiencyWeight = 1. / efficiencyDmeson->at(o2::analysis::findBin(bins, candidate1.pt())); - } - // ========================== Fill mass histo ================================ - if (candidate1.isSelD0() >= selectionFlagD0) { - registry.fill(HIST("hMass"), hfHelper.invMassD0ToPiK(candidate1), candidate1.pt(), efficiencyWeight); - registry.fill(HIST("hMass1D"), hfHelper.invMassD0ToPiK(candidate1), efficiencyWeight); - registry.fill(HIST("hMassD01D"), hfHelper.invMassD0ToPiK(candidate1), efficiencyWeight); - } - if (candidate1.isSelD0bar() >= selectionFlagD0bar) { - registry.fill(HIST("hMass"), hfHelper.invMassD0barToKPi(candidate1), candidate1.pt(), efficiencyWeight); - registry.fill(HIST("hMass1D"), hfHelper.invMassD0barToKPi(candidate1), efficiencyWeight); - registry.fill(HIST("hMassD0bar1D"), hfHelper.invMassD0barToKPi(candidate1), efficiencyWeight); - } - // ========================== Fill general histos ================================ - registry.fill(HIST("hPtCand"), candidate1.pt()); - registry.fill(HIST("hPtProng0"), candidate1.ptProng0()); - registry.fill(HIST("hPtProng1"), candidate1.ptProng1()); - registry.fill(HIST("hEta"), candidate1.eta()); - registry.fill(HIST("hPhi"), candidate1.phi()); - registry.fill(HIST("hY"), hfHelper.yD0(candidate1)); - registry.fill(HIST("hSelectionStatus"), candidate1.isSelD0bar() + (candidate1.isSelD0() * 2)); - registry.fill(HIST("hD0Bin"), poolBin); - - // ============ D-h correlation dedicated section ================================== - - // ========================== track loop starts here ================================ - for (const auto& track : tracks) { - registry.fill(HIST("hTrackCounter"), 1); // fill total no. of tracks - // Remove D0 daughters by checking track indices - bool correlationStatus = false; - if ((candidate1.prong0Id() == track.globalIndex()) || (candidate1.prong1Id() == track.globalIndex())) { - if (!storeAutoCorrelationFlag) { - continue; - } - correlationStatus = true; - } - - registry.fill(HIST("hTrackCounter"), 2); // fill no. of tracks before soft pion removal - - // ========== soft pion removal =================================================== - double invMassDstar1 = 0., invMassDstar2 = 0.; - bool isSoftPiD0 = false, isSoftPiD0bar = false; - auto pSum2 = RecoDecay::p2(candidate1.pVector(), track.pVector()); - auto ePion = track.energy(massPi); - invMassDstar1 = std::sqrt((ePiK + ePion) * (ePiK + ePion) - pSum2); - invMassDstar2 = std::sqrt((eKPi + ePion) * (eKPi + ePion) - pSum2); - - if (candidate1.isSelD0() >= selectionFlagD0) { - if ((std::abs(invMassDstar1 - hfHelper.invMassD0ToPiK(candidate1)) - softPiMass) < ptSoftPionMax) { - isSoftPiD0 = true; - continue; - } - } - - if (candidate1.isSelD0bar() >= selectionFlagD0bar) { - if ((std::abs(invMassDstar2 - hfHelper.invMassD0barToKPi(candidate1)) - softPiMass) < ptSoftPionMax) { - isSoftPiD0bar = true; - continue; - } - } - registry.fill(HIST("hTrackCounter"), 3); // fill no. of tracks after soft pion removal - - int signalStatus = 0; - if ((candidate1.isSelD0() >= selectionFlagD0) && !isSoftPiD0) { - signalStatus += aod::hf_correlation_d0_hadron::ParticleTypeData::D0Only; - } - if ((candidate1.isSelD0bar() >= selectionFlagD0bar) && !isSoftPiD0bar) { - signalStatus += aod::hf_correlation_d0_hadron::ParticleTypeData::D0barOnly; - } - - if (correlateD0WithLeadingParticle) { - if (track.globalIndex() != leadingIndex) { - continue; - } - registry.fill(HIST("hTrackCounter"), 4); // fill no. of tracks have leading particle - } - entryD0HadronPair(getDeltaPhi(track.phi(), candidate1.phi()), - track.eta() - candidate1.eta(), - candidate1.pt(), - track.pt(), - poolBin, - correlationStatus); - entryD0HadronRecoInfo(hfHelper.invMassD0ToPiK(candidate1), hfHelper.invMassD0barToKPi(candidate1), signalStatus); - - } // end inner loop (tracks) - - } // end outer loop - } - PROCESS_SWITCH(HfCorrelatorD0Hadrons, processData, "Process data", false); - - // ================ Process starts for MCRec, same event ======================== - - void processMcRec(SelectedCollisions::iterator const& collision, - SelectedTracksMcRec const& tracks, - soa::Join const&) - { - // protection against empty tables to be sliced - if (selectedD0candidatesMc.size() == 0) { - return; - } - // find leading particle - if (correlateD0WithLeadingParticle) { - leadingIndex = findLeadingParticle(tracks, dcaXYTrackMax.value, dcaZTrackMax.value, etaTrackMax.value); - } - int poolBin = corrBinning.getBin(std::make_tuple(collision.posZ(), collision.multFT0M())); - int nTracks = 0; - if (collision.numContrib() > 1) { - for (const auto& track : tracks) { - if (std::abs(track.eta()) > etaTrackMax) { - continue; - } - if (std::abs(track.dcaXY()) > dcaXYTrackMax || std::abs(track.dcaZ()) > dcaZTrackMax) { - continue; - } - nTracks++; - } - } - registry.fill(HIST("hMultiplicityPreSelection"), nTracks); - if (nTracks < multMin || nTracks > multMax) { - return; - } - registry.fill(HIST("hMultiplicity"), nTracks); - - auto selectedD0CandidatesGroupedMc = selectedD0candidatesMc->sliceByCached(aod::hf_cand::collisionId, collision.globalIndex(), cache); - // MC reco level - bool flagD0 = false; - bool flagD0bar = false; - - for (const auto& candidate1 : selectedD0CandidatesGroupedMc) { - // check decay channel flag for candidate1 - if (!TESTBIT(candidate1.hfflag(), aod::hf_cand_2prong::DecayType::D0ToPiK)) { - continue; - } - if (std::abs(hfHelper.yD0(candidate1)) >= yCandMax || candidate1.pt() <= ptCandMin || candidate1.pt() >= ptTrackMax) { - continue; - } - - double efficiencyWeight = 1.; - if (applyEfficiency) { - efficiencyWeight = 1. / efficiencyDmeson->at(o2::analysis::findBin(bins, candidate1.pt())); - } - - if (std::abs(candidate1.flagMcMatchRec()) == 1 << aod::hf_cand_2prong::DecayType::D0ToPiK) { - // fill per-candidate distributions from D0/D0bar true candidates - registry.fill(HIST("hPtCandRec"), candidate1.pt()); - registry.fill(HIST("hPtProng0Rec"), candidate1.ptProng0()); - registry.fill(HIST("hPtProng1Rec"), candidate1.ptProng1()); - registry.fill(HIST("hEtaRec"), candidate1.eta()); - registry.fill(HIST("hPhiRec"), candidate1.phi()); - registry.fill(HIST("hYRec"), hfHelper.yD0(candidate1)); - registry.fill(HIST("hSelectionStatusRec"), candidate1.isSelD0bar() + (candidate1.isSelD0() * 2)); - } - // fill invariant mass plots from D0/D0bar signal and background candidates - if (candidate1.isSelD0() >= selectionFlagD0) { // only reco as D0 - if (candidate1.flagMcMatchRec() == 1 << aod::hf_cand_2prong::DecayType::D0ToPiK) { // also matched as D0 - registry.fill(HIST("hMassD0RecSig"), hfHelper.invMassD0ToPiK(candidate1), candidate1.pt(), efficiencyWeight); - } else if (candidate1.flagMcMatchRec() == -(1 << aod::hf_cand_2prong::DecayType::D0ToPiK)) { - registry.fill(HIST("hMassD0RecRef"), hfHelper.invMassD0ToPiK(candidate1), candidate1.pt(), efficiencyWeight); - } else { - registry.fill(HIST("hMassD0RecBg"), hfHelper.invMassD0ToPiK(candidate1), candidate1.pt(), efficiencyWeight); - } - } - if (candidate1.isSelD0bar() >= selectionFlagD0bar) { // only reco as D0bar - if (candidate1.flagMcMatchRec() == -(1 << aod::hf_cand_2prong::DecayType::D0ToPiK)) { // also matched as D0bar - registry.fill(HIST("hMassD0barRecSig"), hfHelper.invMassD0barToKPi(candidate1), candidate1.pt(), efficiencyWeight); - } else if (candidate1.flagMcMatchRec() == 1 << aod::hf_cand_2prong::DecayType::D0ToPiK) { - registry.fill(HIST("hMassD0barRecRef"), hfHelper.invMassD0barToKPi(candidate1), candidate1.pt(), efficiencyWeight); - } else { - registry.fill(HIST("hMassD0barRecBg"), hfHelper.invMassD0barToKPi(candidate1), candidate1.pt(), efficiencyWeight); - } - } - - // ===================== Define parameters for soft pion removal ======================== - auto ePiK = RecoDecay::e(candidate1.pVectorProng0(), massPi) + RecoDecay::e(candidate1.pVectorProng1(), massK); - auto eKPi = RecoDecay::e(candidate1.pVectorProng0(), massK) + RecoDecay::e(candidate1.pVectorProng1(), massPi); - - // ============== D-h correlation dedicated section ==================================== - - flagD0 = candidate1.flagMcMatchRec() == (1 << aod::hf_cand_2prong::DecayType::D0ToPiK); // flagD0Signal 'true' if candidate1 matched to D0 (particle) - flagD0bar = candidate1.flagMcMatchRec() == -(1 << aod::hf_cand_2prong::DecayType::D0ToPiK); // flagD0Reflection 'true' if candidate1, selected as D0 (particle), is matched to D0bar (antiparticle) - - // ========== track loop starts here ======================== - - for (const auto& track : tracks) { - registry.fill(HIST("hTrackCounterRec"), 1); // fill total no. of tracks - - // Removing D0 daughters by checking track indices - bool correlationStatus = false; - if ((candidate1.prong0Id() == track.globalIndex()) || (candidate1.prong1Id() == track.globalIndex())) { - if (!storeAutoCorrelationFlag) { - continue; - } - correlationStatus = true; - } - registry.fill(HIST("hTrackCounterRec"), 2); // fill no. of tracks before soft pion removal - - // ===== soft pion removal =================================================== - double invMassDstar1 = 0, invMassDstar2 = 0; - bool isSoftPiD0 = false, isSoftPiD0bar = false; - auto pSum2 = RecoDecay::p2(candidate1.pVector(), track.pVector()); - auto ePion = track.energy(massPi); - invMassDstar1 = std::sqrt((ePiK + ePion) * (ePiK + ePion) - pSum2); - invMassDstar2 = std::sqrt((eKPi + ePion) * (eKPi + ePion) - pSum2); - - if (candidate1.isSelD0() >= selectionFlagD0) { - if ((std::abs(invMassDstar1 - hfHelper.invMassD0ToPiK(candidate1)) - softPiMass) < ptSoftPionMax) { - isSoftPiD0 = true; - continue; - } - } - - if (candidate1.isSelD0bar() >= selectionFlagD0bar) { - if ((std::abs(invMassDstar2 - hfHelper.invMassD0barToKPi(candidate1)) - softPiMass) < ptSoftPionMax) { - isSoftPiD0bar = true; - continue; - } - } - - registry.fill(HIST("hTrackCounterRec"), 3); // fill no. of tracks after soft pion removal - - if (correlateD0WithLeadingParticle) { - if (track.globalIndex() != leadingIndex) { - continue; - } - registry.fill(HIST("hTrackCounterRec"), 4); // fill no. of tracks have leading particle - } - - int signalStatus = 0; - if (flagD0 && (candidate1.isSelD0() >= selectionFlagD0) && !isSoftPiD0) { - SETBIT(signalStatus, aod::hf_correlation_d0_hadron::ParticleTypeMcRec::D0Sig); - } // signal case D0 - if (flagD0bar && (candidate1.isSelD0() >= selectionFlagD0) && !isSoftPiD0) { - SETBIT(signalStatus, aod::hf_correlation_d0_hadron::ParticleTypeMcRec::D0Ref); - } // reflection case D0 - if (!flagD0 && !flagD0bar && (candidate1.isSelD0() >= selectionFlagD0) && !isSoftPiD0) { - SETBIT(signalStatus, aod::hf_correlation_d0_hadron::ParticleTypeMcRec::D0Bg); - } // background case D0 - - if (flagD0bar && (candidate1.isSelD0bar() >= selectionFlagD0bar) && !isSoftPiD0bar) { - SETBIT(signalStatus, aod::hf_correlation_d0_hadron::ParticleTypeMcRec::D0barSig); - } // signal case D0bar - if (flagD0 && (candidate1.isSelD0bar() >= selectionFlagD0bar) && !isSoftPiD0bar) { - SETBIT(signalStatus, aod::hf_correlation_d0_hadron::ParticleTypeMcRec::D0barRef); - } // reflection case D0bar - if (!flagD0 && !flagD0bar && (candidate1.isSelD0bar() >= selectionFlagD0bar) && !isSoftPiD0bar) { - SETBIT(signalStatus, aod::hf_correlation_d0_hadron::ParticleTypeMcRec::D0barBg); - } // background case D0bar - - entryD0HadronPair(getDeltaPhi(track.phi(), candidate1.phi()), - track.eta() - candidate1.eta(), - candidate1.pt(), - track.pt(), - poolBin, - correlationStatus); - entryD0HadronRecoInfo(hfHelper.invMassD0ToPiK(candidate1), hfHelper.invMassD0barToKPi(candidate1), signalStatus); - } // end inner loop (Tracks) - } // end of outer loop (D0) - registry.fill(HIST("hZvtx"), collision.posZ()); - registry.fill(HIST("hMultV0M"), collision.multFT0M()); - } - - PROCESS_SWITCH(HfCorrelatorD0Hadrons, processMcRec, "Process MC Reco mode", true); - - // ================= Process starts for MCGen, same event =================== - - void processMcGen(aod::McCollision const& mcCollision, - soa::Join const& mcParticles) - { - registry.fill(HIST("hEvtCountGen"), 0); - // MC gen level - // find leading particle - if (correlateD0WithLeadingParticle) { - leadingIndex = findLeadingParticleMcGen(mcParticles, etaTrackMax.value, ptTrackMin.value); - } - for (const auto& particle1 : mcParticles) { - // check if the particle is D0 or D0bar (for general plot filling and selection, so both cases are fine) - NOTE: decay channel is not probed! - if (std::abs(particle1.pdgCode()) != Pdg::kD0) { - continue; - } - double yD = RecoDecay::y(particle1.pVector(), MassD0); - if (yCandMax >= 0. && std::abs(yD) > yCandMax) { - continue; - } - if (ptCandMin >= 0. && particle1.pt() < ptCandMin) { - continue; - } - registry.fill(HIST("hPtCandGen"), particle1.pt()); - registry.fill(HIST("hEtaGen"), particle1.eta()); - registry.fill(HIST("hPhiGen"), particle1.phi()); - registry.fill(HIST("hYGen"), yD); - - // =============== D-h correlation dedicated section ===================== - - if (std::abs(particle1.pdgCode()) != Pdg::kD0) { // just checking the particle PDG, not the decay channel (differently from Reco: you have a BR factor btw such levels!) - continue; - } - registry.fill(HIST("hCountD0TriggersGen"), 0, particle1.pt()); // to count trigger D0 (for normalisation) - - for (const auto& particle2 : mcParticles) { - registry.fill(HIST("hTrackCounterGen"), 1); // total no. of tracks - if (std::abs(particle2.eta()) > etaTrackMax) { - continue; - } - if (particle2.pt() < ptTrackMin) { - continue; - } - if ((std::abs(particle2.pdgCode()) != kElectron) && (std::abs(particle2.pdgCode()) != kMuonMinus) && (std::abs(particle2.pdgCode()) != kPiPlus) && (std::abs(particle2.pdgCode()) != kKPlus) && (std::abs(particle2.pdgCode()) != kProton)) { - continue; - } - - // ==============================soft pion removal================================ - registry.fill(HIST("hTrackCounterGen"), 2); // fill before soft pi removal - // method used: indexMother = -1 by default if the mother doesn't match with given PID of the mother. We find mother of pion if it is D* and mother of D0 if it is D*. If they are both positive and they both match each other, then it is detected as a soft pion - - auto indexMotherPi = RecoDecay::getMother(mcParticles, particle2, Pdg::kDStar, true, nullptr, 1); // last arguement 1 is written to consider immediate decay mother only - auto indexMotherD0 = RecoDecay::getMother(mcParticles, particle1, Pdg::kDStar, true, nullptr, 1); - bool correlationStatus = false; - if (std::abs(particle2.pdgCode()) == kPiPlus && indexMotherPi >= 0 && indexMotherD0 >= 0 && indexMotherPi == indexMotherD0) { - if (!storeAutoCorrelationFlag) { - continue; - } - correlationStatus = true; - } - - registry.fill(HIST("hTrackCounterGen"), 3); // fill after soft pion removal - - if (correlateD0WithLeadingParticle) { - if (particle2.globalIndex() != leadingIndex) { - continue; - } - registry.fill(HIST("hTrackCounterGen"), 4); // fill no. of tracks have leading particle - } - - auto getTracksSize = [&mcParticles](aod::McCollision const& /*collision*/) { - int nTracks = 0; - for (const auto& track : mcParticles) { - if (track.isPhysicalPrimary() && std::abs(track.eta()) < 1.0) { - nTracks++; - } - } - return nTracks; - }; - using BinningTypeMcGen = FlexibleBinningPolicy, aod::mccollision::PosZ, decltype(getTracksSize)>; - BinningTypeMcGen corrBinningMcGen{{getTracksSize}, {zPoolBins, multPoolBinsMcGen}, true}; - int poolBin = corrBinningMcGen.getBin(std::make_tuple(mcCollision.posZ(), getTracksSize(mcCollision))); - entryD0HadronPair(getDeltaPhi(particle2.phi(), particle1.phi()), - particle2.eta() - particle1.eta(), - particle1.pt(), - particle2.pt(), - poolBin, - correlationStatus); - entryD0HadronRecoInfo(massD0, massD0, 0); // dummy info - } // end inner loop (Tracks) - } // end outer loop (D0) - } - - PROCESS_SWITCH(HfCorrelatorD0Hadrons, processMcGen, "Process MC Gen mode", false); - - // ====================== Implement Event mixing on Data =================================== - - void processDataMixedEvent(SelectedCollisions const& collisions, - SelectedCandidatesData const& candidates, - SelectedTracks const& tracks) - { - for (const auto& collision : collisions) { - registry.fill(HIST("hMultFT0M"), collision.multFT0M()); - registry.fill(HIST("hZvtx"), collision.posZ()); - } - - auto tracksTuple = std::make_tuple(candidates, tracks); - Pair pairData{corrBinning, numberEventsMixed, -1, collisions, tracksTuple, &cache}; - - for (const auto& [c1, tracks1, c2, tracks2] : pairData) { - // LOGF(info, "Mixed event collisions: Index = (%d, %d), tracks Size: (%d, %d), Z Vertex: (%f, %f), Pool Bin: (%d, %d)", c1.globalIndex(), c2.globalIndex(), tracks1.size(), tracks2.size(), c1.posZ(), c2.posZ(), corrBinning.getBin(std::make_tuple(c1.posZ(), c1.multFT0M())),corrBinning.getBin(std::make_tuple(c2.posZ(), c2.multFT0M()))); // For debug - int poolBin = corrBinning.getBin(std::make_tuple(c2.posZ(), c2.multFT0M())); - for (const auto& [t1, t2] : o2::soa::combinations(o2::soa::CombinationsFullIndexPolicy(tracks1, tracks2))) { - - if (std::abs(hfHelper.yD0(t1)) >= yCandMax) { - continue; - } - - // soft pion removal, signal status 1,3 for D0 and 2,3 for D0bar (SoftPi removed), signal status 11,13 for D0 and 12,13 for D0bar (only SoftPi) - auto ePiK = RecoDecay::e(t1.pVectorProng0(), massPi) + RecoDecay::e(t1.pVectorProng1(), massK); - auto eKPi = RecoDecay::e(t1.pVectorProng0(), massK) + RecoDecay::e(t1.pVectorProng1(), massPi); - double invMassDstar1 = 0., invMassDstar2 = 0.; - bool isSoftPiD0 = false, isSoftPiD0bar = false; - auto pSum2 = RecoDecay::p2(t1.pVector(), t2.pVector()); - auto ePion = t2.energy(massPi); - invMassDstar1 = std::sqrt((ePiK + ePion) * (ePiK + ePion) - pSum2); - invMassDstar2 = std::sqrt((eKPi + ePion) * (eKPi + ePion) - pSum2); - - if (t1.isSelD0() >= selectionFlagD0) { - if ((std::abs(invMassDstar1 - hfHelper.invMassD0ToPiK(t1)) - softPiMass) < ptSoftPionMax) { - isSoftPiD0 = true; - } - } - if (t1.isSelD0bar() >= selectionFlagD0bar) { - if ((std::abs(invMassDstar2 - hfHelper.invMassD0barToKPi(t1)) - softPiMass) < ptSoftPionMax) { - isSoftPiD0bar = true; - } - } - - int signalStatus = 0; - if (t1.isSelD0() >= selectionFlagD0) { - if (!isSoftPiD0) { - signalStatus += aod::hf_correlation_d0_hadron::ParticleTypeData::D0Only; - } else { - signalStatus += aod::hf_correlation_d0_hadron::ParticleTypeData::D0OnlySoftPi; - } - } - if (t1.isSelD0bar() >= selectionFlagD0bar) { - if (!isSoftPiD0bar) { - signalStatus += aod::hf_correlation_d0_hadron::ParticleTypeData::D0barOnly; - } else { - signalStatus += aod::hf_correlation_d0_hadron::ParticleTypeData::D0barOnlySoftPi; - } - } - bool correlationStatus = false; - entryD0HadronPair(getDeltaPhi(t1.phi(), t2.phi()), t1.eta() - t2.eta(), t1.pt(), t2.pt(), poolBin, correlationStatus); - entryD0HadronRecoInfo(hfHelper.invMassD0ToPiK(t1), hfHelper.invMassD0barToKPi(t1), signalStatus); - } - } - } - PROCESS_SWITCH(HfCorrelatorD0Hadrons, processDataMixedEvent, "Process data mixed event", false); - - // ====================== Implement Event mixing on McRec =================================== - - void processMcRecMixedEvent(SelectedCollisions const& collisions, - SelectedCandidatesMcRec const& candidates, - SelectedTracksMcRec const& tracks) - { - auto tracksTuple = std::make_tuple(candidates, tracks); - Pair pairMcRec{corrBinning, numberEventsMixed, -1, collisions, tracksTuple, &cache}; - bool flagD0 = false; - bool flagD0bar = false; - for (const auto& [c1, tracks1, c2, tracks2] : pairMcRec) { - int poolBin = corrBinning.getBin(std::make_tuple(c2.posZ(), c2.multFT0M())); - - for (const auto& [t1, t2] : o2::soa::combinations(o2::soa::CombinationsFullIndexPolicy(tracks1, tracks2))) { - - if (std::abs(hfHelper.yD0(t1)) >= yCandMax) { - continue; - } - - // soft pion removal - auto ePiK = RecoDecay::e(t1.pVectorProng0(), massPi) + RecoDecay::e(t1.pVectorProng1(), massK); - auto eKPi = RecoDecay::e(t1.pVectorProng0(), massK) + RecoDecay::e(t1.pVectorProng1(), massPi); - double invMassDstar1 = 0., invMassDstar2 = 0.; - bool isSoftPiD0 = false, isSoftPiD0bar = false; - auto pSum2 = RecoDecay::p2(t1.pVector(), t2.pVector()); - auto ePion = t2.energy(massPi); - invMassDstar1 = std::sqrt((ePiK + ePion) * (ePiK + ePion) - pSum2); - invMassDstar2 = std::sqrt((eKPi + ePion) * (eKPi + ePion) - pSum2); - - if (t1.isSelD0() >= selectionFlagD0) { - if ((std::abs(invMassDstar1 - hfHelper.invMassD0ToPiK(t1)) - softPiMass) < ptSoftPionMax) { - isSoftPiD0 = true; - } - } - - if (t1.isSelD0bar() >= selectionFlagD0bar) { - if ((std::abs(invMassDstar2 - hfHelper.invMassD0barToKPi(t1)) - softPiMass) < ptSoftPionMax) { - isSoftPiD0bar = true; - } - } - - flagD0 = t1.flagMcMatchRec() == (1 << aod::hf_cand_2prong::DecayType::D0ToPiK); // flagD0Signal 'true' if candidate1 matched to D0 (particle) - flagD0bar = t1.flagMcMatchRec() == -(1 << aod::hf_cand_2prong::DecayType::D0ToPiK); // flagD0Reflection 'true' if candidate1, selected as D0 (particle), is matched to D0bar (antiparticle) - int signalStatus = 0; - - if (flagD0 && (t1.isSelD0() >= selectionFlagD0)) { - if (!isSoftPiD0) { - SETBIT(signalStatus, aod::hf_correlation_d0_hadron::ParticleTypeMcRec::D0Sig); // signalStatus += 1; - } else { - SETBIT(signalStatus, aod::hf_correlation_d0_hadron::ParticleTypeMcRec::SoftPi); // signalStatus += 64; - } - } // signal case D0 - - if (flagD0bar && (t1.isSelD0() >= selectionFlagD0)) { - if (!isSoftPiD0) { - SETBIT(signalStatus, aod::hf_correlation_d0_hadron::ParticleTypeMcRec::D0Ref); // signalStatus += 2; - } else { - SETBIT(signalStatus, aod::hf_correlation_d0_hadron::ParticleTypeMcRec::SoftPi); // signalStatus += 64; - } - } // reflection case D0 - - if (!flagD0 && !flagD0bar && (t1.isSelD0() >= selectionFlagD0)) { - if (!isSoftPiD0) { - SETBIT(signalStatus, aod::hf_correlation_d0_hadron::ParticleTypeMcRec::D0Bg); // signalStatus += 4; - } else { - SETBIT(signalStatus, aod::hf_correlation_d0_hadron::ParticleTypeMcRec::SoftPi); - } - } // background case D0 - - if (flagD0bar && (t1.isSelD0bar() >= selectionFlagD0bar)) { - if (!isSoftPiD0bar) { - SETBIT(signalStatus, aod::hf_correlation_d0_hadron::ParticleTypeMcRec::D0barSig); // signalStatus += 8; - } else { - SETBIT(signalStatus, aod::hf_correlation_d0_hadron::ParticleTypeMcRec::SoftPi); - } - } // signal case D0bar - - if (flagD0 && (t1.isSelD0bar() >= selectionFlagD0bar)) { - if (!isSoftPiD0bar) { - SETBIT(signalStatus, aod::hf_correlation_d0_hadron::ParticleTypeMcRec::D0barRef); // signalStatus += 16; - } else { - SETBIT(signalStatus, aod::hf_correlation_d0_hadron::ParticleTypeMcRec::SoftPi); - } - } // reflection case D0bar - - if (!flagD0 && !flagD0bar && (t1.isSelD0bar() >= selectionFlagD0bar)) { - if (!isSoftPiD0bar) { - SETBIT(signalStatus, aod::hf_correlation_d0_hadron::ParticleTypeMcRec::D0barBg); // signalStatus += 32; - } else { - SETBIT(signalStatus, aod::hf_correlation_d0_hadron::ParticleTypeMcRec::SoftPi); - } - } // background case D0bar - registry.fill(HIST("hSignalStatusMERec"), signalStatus); - bool correlationStatus = false; - entryD0HadronPair(getDeltaPhi(t1.phi(), t2.phi()), t1.eta() - t2.eta(), t1.pt(), t2.pt(), poolBin, correlationStatus); - entryD0HadronRecoInfo(hfHelper.invMassD0ToPiK(t1), hfHelper.invMassD0barToKPi(t1), signalStatus); - } - } - } - PROCESS_SWITCH(HfCorrelatorD0Hadrons, processMcRecMixedEvent, "Process Mixed Event MCRec", false); - - // ====================== Implement Event mixing on McGen =================================== - - void processMcGenMixedEvent(SelectedCollisionsMcGen const& collisions, - SelectedParticlesMcGen const& mcParticles) - { - - auto getTracksSize = [&mcParticles, this](SelectedCollisionsMcGen::iterator const& collision) { - int nTracks = 0; - auto associatedTracks = mcParticles.sliceByCached(o2::aod::mcparticle::mcCollisionId, collision.globalIndex(), this->cache); - for (const auto& track : associatedTracks) { - if (track.isPhysicalPrimary() && std::abs(track.eta()) < 1.0) { - nTracks++; - } - } - return nTracks; - }; - - using BinningTypeMcGen = FlexibleBinningPolicy, aod::mccollision::PosZ, decltype(getTracksSize)>; - BinningTypeMcGen corrBinningMcGen{{getTracksSize}, {zPoolBins, multPoolBinsMcGen}, true}; - - auto tracksTuple = std::make_tuple(mcParticles, mcParticles); - Pair pairMcGen{corrBinningMcGen, numberEventsMixed, -1, collisions, tracksTuple, &cache}; - - for (const auto& [c1, tracks1, c2, tracks2] : pairMcGen) { - for (const auto& [t1, t2] : o2::soa::combinations(o2::soa::CombinationsFullIndexPolicy(tracks1, tracks2))) { - - // Check track t1 is D0 - if (std::abs(t1.pdgCode()) != Pdg::kD0) { - continue; - } - - double yD = RecoDecay::y(t1.pVector(), MassD0); - if (std::abs(yD) >= yCandMax || t1.pt() <= ptCandMin || std::abs(t2.eta()) >= etaTrackMax || t2.pt() <= ptTrackMin) { - continue; - } - if ((std::abs(t2.pdgCode()) != kElectron) && (std::abs(t2.pdgCode()) != kMuonMinus) && (std::abs(t2.pdgCode()) != kPiPlus) && (std::abs(t2.pdgCode()) != kKPlus) && (std::abs(t2.pdgCode()) != kProton)) { - continue; - } - if (!t2.isPhysicalPrimary()) { - continue; - } - - // ==============================soft pion removal================================ - // method used: indexMother = -1 by default if the mother doesn't match with given PID of the mother. We find mother of pion if it is D* and mother of D0 if it is D*. If they are both positive and they both match each other, then it is detected as a soft pion - - auto indexMotherPi = RecoDecay::getMother(mcParticles, t2, Pdg::kDStar, true, nullptr, 1); // last arguement 1 is written to consider immediate decay mother only - auto indexMotherD0 = RecoDecay::getMother(mcParticles, t1, Pdg::kDStar, true, nullptr, 1); - if (std::abs(t2.pdgCode()) == kPiPlus && indexMotherPi >= 0 && indexMotherD0 >= 0 && indexMotherPi == indexMotherD0) { - continue; - } - int poolBin = corrBinningMcGen.getBin(std::make_tuple(c2.posZ(), getTracksSize(c2))); - bool correlationStatus = false; - entryD0HadronPair(getDeltaPhi(t2.phi(), t1.phi()), t2.eta() - t1.eta(), t1.pt(), t2.pt(), poolBin, correlationStatus); - entryD0HadronRecoInfo(massD0, massD0, 0); // dummy info - } - } - } - PROCESS_SWITCH(HfCorrelatorD0Hadrons, processMcGenMixedEvent, "Process Mixed Event MCGen", false); -}; - -WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) -{ - return WorkflowSpec{ - adaptAnalysisTask(cfgc), - adaptAnalysisTask(cfgc)}; -} From 6cc727800e7ea4c94e11fef36ec77b4480ae99ce Mon Sep 17 00:00:00 2001 From: samrangy Date: Fri, 29 Nov 2024 18:32:41 +0100 Subject: [PATCH 06/11] Add requireGlobalTrackWoDCAInFilter() to trackFilter --- .../HFC/TableProducer/correlatorD0Hadrons.cxx | 915 ++++++++++++++++++ 1 file changed, 915 insertions(+) create mode 100644 PWGHF/HFC/TableProducer/correlatorD0Hadrons.cxx diff --git a/PWGHF/HFC/TableProducer/correlatorD0Hadrons.cxx b/PWGHF/HFC/TableProducer/correlatorD0Hadrons.cxx new file mode 100644 index 00000000000..6ad93e81d42 --- /dev/null +++ b/PWGHF/HFC/TableProducer/correlatorD0Hadrons.cxx @@ -0,0 +1,915 @@ +// 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 correlatorD0Hadrons.cxx +/// \brief D0-Hadron correlator task - data-like, MC-reco and MC-kine analyses. +/// +/// \author Samrangy Sadhu , INFN Bari +/// \author Swapnesh Santosh Khade , IIT Indore + +#include "CommonConstants/PhysicsConstants.h" +#include "Framework/AnalysisTask.h" +#include "Framework/HistogramRegistry.h" +#include "Framework/runDataProcessing.h" + +#include "Common/Core/TrackSelection.h" +#include "Common/DataModel/Centrality.h" +#include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/Multiplicity.h" +#include "Common/DataModel/TrackSelectionTables.h" + +#include "PWGHF/Core/HfHelper.h" +#include "PWGHF/DataModel/CandidateReconstructionTables.h" +#include "PWGHF/DataModel/CandidateSelectionTables.h" +#include "PWGHF/HFC/DataModel/CorrelationTables.h" + +using namespace o2; +using namespace o2::analysis; +using namespace o2::constants::physics; +using namespace o2::framework; +using namespace o2::framework::expressions; + +/// +/// Returns deltaPhi value in range [-pi/2., 3.*pi/2], typically used for correlation studies +/// +double getDeltaPhi(double phiHadron, double phiD) +{ + return RecoDecay::constrainAngle(phiHadron - phiD, -o2::constants::math::PIHalf); +} + +const int nPtBinsMassAndEfficiency = o2::analysis::hf_cuts_d0_to_pi_k::nBinsPt; +const double efficiencyDmesonDefault[nPtBinsMassAndEfficiency] = {}; +auto vecEfficiencyDmeson = std::vector{efficiencyDmesonDefault, efficiencyDmesonDefault + nPtBinsMassAndEfficiency}; + +// histogram binning definition +const int massAxisNBins = 200; +const double massAxisMin = 1.3848; +const double massAxisMax = 2.3848; +const int phiAxisNBins = 32; +const double phiAxisMin = 0.; +const double phiAxisMax = o2::constants::math::TwoPI; +const int yAxisNBins = 100; +const double yAxisMin = -5.; +const double yAxisMax = 5.; +const int ptDAxisNBins = 180; +const double ptDAxisMin = 0.; +const double ptDAxisMax = 36.; + +// Types +using BinningType = ColumnBinningPolicy>; +using SelectedCollisions = soa::Filtered>; +using SelectedTracks = soa::Filtered>; +using SelectedCandidatesData = soa::Filtered>; +using SelectedTracksMcRec = soa::Filtered>; +using SelectedCandidatesMcRec = soa::Filtered>; +using SelectedCollisionsMcGen = soa::Filtered>; +using SelectedParticlesMcGen = soa::Join; + +// Code to select collisions with at least one D0 +struct HfCorrelatorD0HadronsSelection { + SliceCache cache; + + Produces d0Sel; + + Configurable selectionFlagD0{"selectionFlagD0", 1, "Selection Flag for D0"}; + Configurable selectionFlagD0bar{"selectionFlagD0bar", 1, "Selection Flag for D0bar"}; + Configurable yCandMax{"yCandMax", 4.0, "max. cand. rapidity"}; + Configurable ptCandMin{"ptCandMin", -1., "min. cand. pT"}; + + HfHelper hfHelper; + + Preslice perCol = aod::hf_cand::collisionId; + Partition> selectedD0Candidates = aod::hf_sel_candidate_d0::isSelD0 >= selectionFlagD0 || aod::hf_sel_candidate_d0::isSelD0bar >= selectionFlagD0bar; + Partition> selectedD0candidatesMc = aod::hf_sel_candidate_d0::isSelD0 >= selectionFlagD0 || aod::hf_sel_candidate_d0::isSelD0bar >= selectionFlagD0bar; + + void processD0SelectionData(aod::Collision const& collision, + soa::Join const&) + { + bool isD0Found = 0; + if (selectedD0Candidates.size() > 0) { + auto selectedD0CandidatesGrouped = selectedD0Candidates->sliceByCached(aod::hf_cand::collisionId, collision.globalIndex(), cache); + + for (const auto& candidate : selectedD0CandidatesGrouped) { + // check decay channel flag for candidate + if (!TESTBIT(candidate.hfflag(), aod::hf_cand_2prong::DecayType::D0ToPiK)) { + continue; + } + if (std::abs(hfHelper.yD0(candidate)) > yCandMax || candidate.pt() < ptCandMin) { + continue; + } + isD0Found = 1; + break; + } + } + d0Sel(isD0Found); + } + PROCESS_SWITCH(HfCorrelatorD0HadronsSelection, processD0SelectionData, "Process D0 Selection Data", false); + + void processD0SelectionMcRec(aod::Collision const& collision, + soa::Join const&) + { + bool isD0Found = 0; + if (selectedD0candidatesMc.size() > 0) { + auto selectedD0CandidatesGroupedMc = selectedD0candidatesMc->sliceByCached(aod::hf_cand::collisionId, collision.globalIndex(), cache); + for (const auto& candidate : selectedD0CandidatesGroupedMc) { + // check decay channel flag for candidate + if (!TESTBIT(candidate.hfflag(), aod::hf_cand_2prong::DecayType::D0ToPiK)) { + continue; + } + if (std::abs(hfHelper.yD0(candidate)) > yCandMax || candidate.pt() < ptCandMin) { + continue; + } + isD0Found = 1; + break; + } + } + d0Sel(isD0Found); + } + PROCESS_SWITCH(HfCorrelatorD0HadronsSelection, processD0SelectionMcRec, "Process D0 Selection MCRec", true); + + void processD0SelectionMcGen(aod::McCollision const&, + aod::McParticles const& mcParticles) + { + bool isD0Found = 0; + for (const auto& particle : mcParticles) { + if (std::abs(particle.pdgCode()) != Pdg::kD0) { + continue; + } + double yD = RecoDecay::y(particle.pVector(), MassD0); + if (std::abs(yD) > yCandMax || particle.pt() < ptCandMin) { + continue; + } + isD0Found = 1; + break; + } + d0Sel(isD0Found); + } + PROCESS_SWITCH(HfCorrelatorD0HadronsSelection, processD0SelectionMcGen, "Process D0 Selection MCGen", false); +}; + +/// D0-Hadron correlation pair builder - for real data and data-like analysis (i.e. reco-level w/o matching request via MC truth) +struct HfCorrelatorD0Hadrons { + SliceCache cache; + + Produces entryD0HadronPair; + Produces entryD0HadronRecoInfo; + + Configurable selectionFlagD0{"selectionFlagD0", 1, "Selection Flag for D0"}; + Configurable selectionFlagD0bar{"selectionFlagD0bar", 1, "Selection Flag for D0bar"}; + Configurable yCandMax{"yCandMax", 0.8, "max. cand. rapidity"}; + Configurable etaTrackMax{"etaTrackMax", 0.8, "max. eta of tracks"}; + Configurable dcaXYTrackMax{"dcaXYTrackMax", 1., "max. DCAxy of tracks"}; + Configurable dcaZTrackMax{"dcaZTrackMax", 1., "max. DCAz of tracks"}; + Configurable ptCandMin{"ptCandMin", 1., "min. cand. pT"}; + Configurable ptTrackMin{"ptTrackMin", 0.3, "min. track pT"}; + Configurable ptTrackMax{"ptTrackMax", 99., "max. track pT"}; + Configurable> bins{"ptBinsForMassAndEfficiency", std::vector{o2::analysis::hf_cuts_d0_to_pi_k::vecBinsPt}, "pT bin limits for candidate mass plots and efficiency"}; + Configurable> efficiencyDmeson{"efficiencyDmeson", std::vector{vecEfficiencyDmeson}, "Efficiency values for D0 meson"}; + Configurable applyEfficiency{"efficiencyFlagD", 1, "Flag for applying D-meson efficiency weights"}; + Configurable multMin{"multMin", 0., "minimum multiplicity accepted"}; + Configurable multMax{"multMax", 10000., "maximum multiplicity accepted"}; + Configurable ptSoftPionMax{"ptSoftPionMax", 3 * 800. * pow(10., -6.), "max. pT cut for soft pion identification"}; + Configurable correlateD0WithLeadingParticle{"correlateD0WithLeadingParticle", false, "Switch for correlation of D0 mesons with leading particle only"}; + Configurable storeAutoCorrelationFlag{"storeAutoCorrelationFlag", false, "Store flag that indicates if the track is paired to its D-meson mother instead of skipping it"}; + Configurable numberEventsMixed{"numberEventsMixed", 5, "Number of events mixed in ME process"}; + ConfigurableAxis zPoolBins{"zPoolBins", {VARIABLE_WIDTH, -10.0f, -2.5f, 2.5f, 10.0f}, "z vertex position pools"}; + ConfigurableAxis multPoolBins{"multPoolBins", {VARIABLE_WIDTH, 0.0f, 2000.0f, 6000.0f, 10000.0f}, "event multiplicity pools (FT0M)"}; + ConfigurableAxis multPoolBinsMcGen{"multPoolBinsMcGen", {VARIABLE_WIDTH, 0.0f, 20.0f, 50.0f, 500.0f}, "Mixing bins - MC multiplicity"}; // In MCGen multiplicity is defined by counting tracks + + HfHelper hfHelper; + BinningType corrBinning{{zPoolBins, multPoolBins}, true}; + + int leadingIndex = 0; + double massD0{0.}; + double massPi{0.}; + double massK{0.}; + double softPiMass = 0.14543; // pion mass + Q-value of the D*->D0pi decay + + Preslice perCol = aod::hf_cand::collisionId; + + Partition> selectedD0Candidates = aod::hf_sel_candidate_d0::isSelD0 >= selectionFlagD0 || aod::hf_sel_candidate_d0::isSelD0bar >= selectionFlagD0bar; + Partition> selectedD0candidatesMc = aod::hf_sel_candidate_d0::isSelD0 >= selectionFlagD0 || aod::hf_sel_candidate_d0::isSelD0bar >= selectionFlagD0bar; + + Filter collisionFilter = aod::hf_selection_dmeson_collision::dmesonSel == true; + Filter trackFilter = requireGlobalTrackWoDCAInFilter() && (nabs(aod::track::eta) < etaTrackMax) && (aod::track::pt > ptTrackMin) && (aod::track::pt < ptTrackMax) && (nabs(aod::track::dcaXY) < dcaXYTrackMax) && (nabs(aod::track::dcaZ) < dcaZTrackMax); + Filter d0Filter = (aod::hf_sel_candidate_d0::isSelD0 >= 1) || (aod::hf_sel_candidate_d0::isSelD0bar >= 1); + Filter collisionFilterGen = aod::hf_selection_dmeson_collision::dmesonSel == true; + Filter particlesFilter = nabs(aod::mcparticle::pdgCode) == static_cast(Pdg::kD0) || ((aod::mcparticle::flags & (uint8_t)o2::aod::mcparticle::enums::PhysicalPrimary) == (uint8_t)o2::aod::mcparticle::enums::PhysicalPrimary); + + HistogramRegistry registry{ + "registry", + // NOTE: use hMassD0 for trigger normalisation (S*0.955), and hMass2DCorrelationPairs (in final task) for 2D-sideband-subtraction purposes + {{"hPtCand", "D0,D0bar candidates;candidate #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{ptDAxisNBins, ptDAxisMin, ptDAxisMax}}}}, + {"hPtProng0", "D0,D0bar candidates;prong 0 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{ptDAxisNBins, ptDAxisMin, ptDAxisMax}}}}, + {"hPtProng1", "D0,D0bar candidates;prong 1 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{ptDAxisNBins, ptDAxisMin, ptDAxisMax}}}}, + {"hSelectionStatus", "D0,D0bar candidates;selection status;entries", {HistType::kTH1F, {{4, -0.5, 3.5}}}}, + {"hEta", "D0,D0bar candidates;candidate #it{#eta};entries", {HistType::kTH1F, {{yAxisNBins, yAxisMin, yAxisMax}}}}, + {"hPhi", "D0,D0bar candidates;candidate #it{#varphi};entries", {HistType::kTH1F, {{phiAxisNBins, phiAxisMin, phiAxisMax}}}}, + {"hY", "D0,D0bar candidates;candidate #it{y};entries", {HistType::kTH1F, {{yAxisNBins, yAxisMin, yAxisMax}}}}, + {"hMultiplicityPreSelection", "multiplicity prior to selection;multiplicity;entries", {HistType::kTH1F, {{10000, 0., 10000.}}}}, + {"hMultiplicity", "multiplicity;multiplicity;entries", {HistType::kTH1F, {{10000, 0., 10000.}}}}, + {"hPtCandRec", "D0,D0bar candidates - MC reco;candidate #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{ptDAxisNBins, ptDAxisMin, ptDAxisMax}}}}, + {"hPtProng0Rec", "D0,D0bar candidates - MC reco;prong 0 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{ptDAxisNBins, ptDAxisMin, ptDAxisMax}}}}, + {"hPtProng1Rec", "D0,D0bar candidates - MC reco;prong 1 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{ptDAxisNBins, ptDAxisMin, ptDAxisMax}}}}, + {"hSelectionStatusRec", "D0,D0bar candidates - MC reco;selection status;entries", {HistType::kTH1F, {{4, -0.5, 3.5}}}}, + {"hSignalStatusMERec", "Signal Status - MC reco ME;candidate sidnalStatus;entries", {HistType::kTH1F, {{200, 0, 200}}}}, + {"hEtaRec", "D0,D0bar candidates - MC reco;candidate #it{#eta};entries", {HistType::kTH1F, {{yAxisNBins, yAxisMin, yAxisMax}}}}, + {"hPhiRec", "D0,D0bar candidates - MC reco;candidate #it{#varphi};entries", {HistType::kTH1F, {{phiAxisNBins, phiAxisMin, phiAxisMax}}}}, + {"hYRec", "D0,D0bar candidates - MC reco;candidate #it{y};entries", {HistType::kTH1F, {{yAxisNBins, yAxisMin, yAxisMax}}}}, + {"hEvtCountGen", "Event counter - MC gen;;entries", {HistType::kTH1F, {{1, -0.5, 0.5}}}}, + {"hPtCandGen", "D0,D0bar particles - MC gen;particle #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{ptDAxisNBins, ptDAxisMin, ptDAxisMax}}}}, + {"hEtaGen", "D0,D0bar particles - MC gen;particle #it{#eta};entries", {HistType::kTH1F, {{yAxisNBins, yAxisMin, yAxisMax}}}}, + {"hPhiGen", "D0,D0bar particles - MC gen;particle #it{#varphi};entries", {HistType::kTH1F, {{phiAxisNBins, phiAxisMin, phiAxisMax}}}}, + {"hYGen", "D0,D0bar candidates - MC gen;candidate #it{y};entries", {HistType::kTH1F, {{yAxisNBins, yAxisMin, yAxisMax}}}}, + {"hTrackCounter", "soft pion counter - Data", {HistType::kTH1F, {{5, 0., 5.}}}}, + {"hTrackCounterRec", "soft pion counter - MC rec", {HistType::kTH1F, {{5, 0., 5.}}}}, + {"hTrackCounterGen", "soft pion counter - MC gen", {HistType::kTH1F, {{5, 0., 5.}}}}, + {"hMultV0M", "multiplicity;multiplicity;entries", {HistType::kTH1F, {{10000, 0., 10000.}}}}, + {"hZvtx", "z vertex;z vertex;entries", {HistType::kTH1F, {{200, -20., 20.}}}}, + {"hMultFT0M", "Multiplicity FT0M", {HistType::kTH1F, {{10000, 0., 10000.}}}}, + {"hD0Bin", "D0 selected in pool Bin;pool Bin;entries", {HistType::kTH1F, {{9, 0., 9.}}}}, + {"hTracksBin", "Tracks selected in pool Bin;pool Bin;entries", {HistType::kTH1F, {{9, 0., 9.}}}}}}; + + void init(InitContext&) + { + massD0 = MassD0; + massPi = MassPiPlus; + massK = MassKPlus; + + auto vbins = (std::vector)bins; + registry.add("hMass", "D0,D0bar candidates;inv. mass (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{massAxisNBins, massAxisMin, massAxisMax}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hMass1D", "D0,D0bar candidates;inv. mass (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH1F, {{massAxisNBins, massAxisMin, massAxisMax}}}); + registry.add("hMassD01D", "D0,D0bar candidates;inv. mass (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH1F, {{massAxisNBins, massAxisMin, massAxisMax}}}); + registry.add("hMassD0bar1D", "D0,D0bar candidates;inv. mass (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH1F, {{massAxisNBins, massAxisMin, massAxisMax}}}); + // mass histogram for D0 signal candidates + registry.add("hMassD0RecSig", "D0 signal candidates - MC reco;inv. mass (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{massAxisNBins, massAxisMin, massAxisMax}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + // mass histogram for D0 Reflection candidates + registry.add("hMassD0RecRef", "D0 reflection candidates - MC reco;inv. mass (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{massAxisNBins, massAxisMin, massAxisMax}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + // mass histogram for D0 background candidates + registry.add("hMassD0RecBg", "D0 background candidates - MC reco;inv. mass (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{massAxisNBins, massAxisMin, massAxisMax}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + // mass histogram for D0bar signal candidates + registry.add("hMassD0barRecSig", "D0bar signal candidates - MC reco;inv. mass D0bar only (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{massAxisNBins, massAxisMin, massAxisMax}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + // mass histogram for D0bar Reflection candidates + registry.add("hMassD0barRecRef", "D0bar reflection candidates - MC reco;inv. mass D0bar only (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{massAxisNBins, massAxisMin, massAxisMax}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + // mass histogram for D0bar background candidates + registry.add("hMassD0barRecBg", "D0bar background candidates - MC reco;inv. mass D0bar only (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{massAxisNBins, massAxisMin, massAxisMax}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hCountD0TriggersGen", "D0 trigger particles - MC gen;;N of trigger D0", {HistType::kTH2F, {{1, -0.5, 0.5}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + } + + // Find Leading Particle + template + int findLeadingParticle(TTracks const& tracks) + { + auto leadingParticle = tracks.begin(); + for (auto const& track : tracks) { + if (std::abs(track.dcaXY()) >= 1. || std::abs(track.dcaZ()) >= 1.) { + continue; + } + if (track.pt() > leadingParticle.pt()) { + leadingParticle = track; + } + } + int leadingIndex = leadingParticle.globalIndex(); + return leadingIndex; + } + // ======= Process starts for Data, Same event ============ + + /// D0-h correlation pair builder - for real data and data-like analysis (i.e. reco-level w/o matching request via MC truth) + void processData(SelectedCollisions::iterator const& collision, + SelectedTracks const& tracks, + soa::Join const&) + { + // protection against empty tables to be sliced + if (selectedD0Candidates.size() == 0) { + return; + } + // find leading particle + if (correlateD0WithLeadingParticle) { + leadingIndex = findLeadingParticle(tracks); + } + + int poolBin = corrBinning.getBin(std::make_tuple(collision.posZ(), collision.multFT0M())); + int nTracks = 0; + if (collision.numContrib() > 1) { + for (const auto& track : tracks) { + if (std::abs(track.eta()) > etaTrackMax) { + continue; + } + if (std::abs(track.dcaXY()) > dcaXYTrackMax || std::abs(track.dcaZ()) > dcaZTrackMax) { + continue; + } + nTracks++; + } + } + registry.fill(HIST("hMultiplicityPreSelection"), nTracks); + registry.fill(HIST("hZvtx"), collision.posZ()); + registry.fill(HIST("hMultFT0M"), collision.multFT0M()); + if (nTracks < multMin || nTracks > multMax) { + return; + } + registry.fill(HIST("hMultiplicity"), nTracks); + + auto selectedD0CandidatesGrouped = selectedD0Candidates->sliceByCached(aod::hf_cand::collisionId, collision.globalIndex(), cache); + + for (const auto& candidate1 : selectedD0CandidatesGrouped) { + if (std::abs(hfHelper.yD0(candidate1)) >= yCandMax || candidate1.pt() <= ptCandMin || candidate1.pt() >= ptTrackMax) { + continue; + } + // check decay channel flag for candidate1 + if (!TESTBIT(candidate1.hfflag(), aod::hf_cand_2prong::DecayType::D0ToPiK)) { + continue; + } + + // ========================== Define parameters for soft pion removal ================================ + auto ePiK = RecoDecay::e(candidate1.pVectorProng0(), massPi) + RecoDecay::e(candidate1.pVectorProng1(), massK); + auto eKPi = RecoDecay::e(candidate1.pVectorProng0(), massK) + RecoDecay::e(candidate1.pVectorProng1(), massPi); + + // ========================== trigger efficiency ================================ + double efficiencyWeight = 1.; + if (applyEfficiency) { + efficiencyWeight = 1. / efficiencyDmeson->at(o2::analysis::findBin(bins, candidate1.pt())); + } + // ========================== Fill mass histo ================================ + if (candidate1.isSelD0() >= selectionFlagD0) { + registry.fill(HIST("hMass"), hfHelper.invMassD0ToPiK(candidate1), candidate1.pt(), efficiencyWeight); + registry.fill(HIST("hMass1D"), hfHelper.invMassD0ToPiK(candidate1), efficiencyWeight); + registry.fill(HIST("hMassD01D"), hfHelper.invMassD0ToPiK(candidate1), efficiencyWeight); + } + if (candidate1.isSelD0bar() >= selectionFlagD0bar) { + registry.fill(HIST("hMass"), hfHelper.invMassD0barToKPi(candidate1), candidate1.pt(), efficiencyWeight); + registry.fill(HIST("hMass1D"), hfHelper.invMassD0barToKPi(candidate1), efficiencyWeight); + registry.fill(HIST("hMassD0bar1D"), hfHelper.invMassD0barToKPi(candidate1), efficiencyWeight); + } + // ========================== Fill general histos ================================ + registry.fill(HIST("hPtCand"), candidate1.pt()); + registry.fill(HIST("hPtProng0"), candidate1.ptProng0()); + registry.fill(HIST("hPtProng1"), candidate1.ptProng1()); + registry.fill(HIST("hEta"), candidate1.eta()); + registry.fill(HIST("hPhi"), candidate1.phi()); + registry.fill(HIST("hY"), hfHelper.yD0(candidate1)); + registry.fill(HIST("hSelectionStatus"), candidate1.isSelD0bar() + (candidate1.isSelD0() * 2)); + registry.fill(HIST("hD0Bin"), poolBin); + + // ============ D-h correlation dedicated section ================================== + + // ========================== track loop starts here ================================ + for (const auto& track : tracks) { + registry.fill(HIST("hTrackCounter"), 1); // fill total no. of tracks + // Remove D0 daughters by checking track indices + bool correlationStatus = false; + if ((candidate1.prong0Id() == track.globalIndex()) || (candidate1.prong1Id() == track.globalIndex())) { + if (!storeAutoCorrelationFlag) { + continue; + } + correlationStatus = true; + } + + registry.fill(HIST("hTrackCounter"), 2); // fill no. of tracks before soft pion removal + + // ========== soft pion removal =================================================== + double invMassDstar1 = 0., invMassDstar2 = 0.; + bool isSoftPiD0 = false, isSoftPiD0bar = false; + auto pSum2 = RecoDecay::p2(candidate1.pVector(), track.pVector()); + auto ePion = track.energy(massPi); + invMassDstar1 = std::sqrt((ePiK + ePion) * (ePiK + ePion) - pSum2); + invMassDstar2 = std::sqrt((eKPi + ePion) * (eKPi + ePion) - pSum2); + + if (candidate1.isSelD0() >= selectionFlagD0) { + if ((std::abs(invMassDstar1 - hfHelper.invMassD0ToPiK(candidate1)) - softPiMass) < ptSoftPionMax) { + isSoftPiD0 = true; + continue; + } + } + + if (candidate1.isSelD0bar() >= selectionFlagD0bar) { + if ((std::abs(invMassDstar2 - hfHelper.invMassD0barToKPi(candidate1)) - softPiMass) < ptSoftPionMax) { + isSoftPiD0bar = true; + continue; + } + } + registry.fill(HIST("hTrackCounter"), 3); // fill no. of tracks after soft pion removal + + int signalStatus = 0; + if ((candidate1.isSelD0() >= selectionFlagD0) && !isSoftPiD0) { + signalStatus += aod::hf_correlation_d0_hadron::ParticleTypeData::D0Only; + } + if ((candidate1.isSelD0bar() >= selectionFlagD0bar) && !isSoftPiD0bar) { + signalStatus += aod::hf_correlation_d0_hadron::ParticleTypeData::D0barOnly; + } + + if (correlateD0WithLeadingParticle) { + if (track.globalIndex() != leadingIndex) { + continue; + } + registry.fill(HIST("hTrackCounter"), 4); // fill no. of tracks have leading particle + } + entryD0HadronPair(getDeltaPhi(track.phi(), candidate1.phi()), + track.eta() - candidate1.eta(), + candidate1.pt(), + track.pt(), + poolBin, + correlationStatus); + entryD0HadronRecoInfo(hfHelper.invMassD0ToPiK(candidate1), hfHelper.invMassD0barToKPi(candidate1), signalStatus); + + } // end inner loop (tracks) + + } // end outer loop + } + PROCESS_SWITCH(HfCorrelatorD0Hadrons, processData, "Process data", false); + + // ================ Process starts for MCRec, same event ======================== + + void processMcRec(SelectedCollisions::iterator const& collision, + SelectedTracksMcRec const& tracks, + soa::Join const&) + { + // protection against empty tables to be sliced + if (selectedD0candidatesMc.size() == 0) { + return; + } + // find leading particle + if (correlateD0WithLeadingParticle) { + leadingIndex = findLeadingParticle(tracks); + } + int poolBin = corrBinning.getBin(std::make_tuple(collision.posZ(), collision.multFT0M())); + int nTracks = 0; + if (collision.numContrib() > 1) { + for (const auto& track : tracks) { + if (std::abs(track.eta()) > etaTrackMax) { + continue; + } + if (std::abs(track.dcaXY()) > dcaXYTrackMax || std::abs(track.dcaZ()) > dcaZTrackMax) { + continue; + } + nTracks++; + } + } + registry.fill(HIST("hMultiplicityPreSelection"), nTracks); + if (nTracks < multMin || nTracks > multMax) { + return; + } + registry.fill(HIST("hMultiplicity"), nTracks); + + auto selectedD0CandidatesGroupedMc = selectedD0candidatesMc->sliceByCached(aod::hf_cand::collisionId, collision.globalIndex(), cache); + // MC reco level + bool flagD0 = false; + bool flagD0bar = false; + + for (const auto& candidate1 : selectedD0CandidatesGroupedMc) { + // check decay channel flag for candidate1 + if (!TESTBIT(candidate1.hfflag(), aod::hf_cand_2prong::DecayType::D0ToPiK)) { + continue; + } + if (std::abs(hfHelper.yD0(candidate1)) >= yCandMax || candidate1.pt() <= ptCandMin || candidate1.pt() >= ptTrackMax) { + continue; + } + + double efficiencyWeight = 1.; + if (applyEfficiency) { + efficiencyWeight = 1. / efficiencyDmeson->at(o2::analysis::findBin(bins, candidate1.pt())); + } + + if (std::abs(candidate1.flagMcMatchRec()) == 1 << aod::hf_cand_2prong::DecayType::D0ToPiK) { + // fill per-candidate distributions from D0/D0bar true candidates + registry.fill(HIST("hPtCandRec"), candidate1.pt()); + registry.fill(HIST("hPtProng0Rec"), candidate1.ptProng0()); + registry.fill(HIST("hPtProng1Rec"), candidate1.ptProng1()); + registry.fill(HIST("hEtaRec"), candidate1.eta()); + registry.fill(HIST("hPhiRec"), candidate1.phi()); + registry.fill(HIST("hYRec"), hfHelper.yD0(candidate1)); + registry.fill(HIST("hSelectionStatusRec"), candidate1.isSelD0bar() + (candidate1.isSelD0() * 2)); + } + // fill invariant mass plots from D0/D0bar signal and background candidates + if (candidate1.isSelD0() >= selectionFlagD0) { // only reco as D0 + if (candidate1.flagMcMatchRec() == 1 << aod::hf_cand_2prong::DecayType::D0ToPiK) { // also matched as D0 + registry.fill(HIST("hMassD0RecSig"), hfHelper.invMassD0ToPiK(candidate1), candidate1.pt(), efficiencyWeight); + } else if (candidate1.flagMcMatchRec() == -(1 << aod::hf_cand_2prong::DecayType::D0ToPiK)) { + registry.fill(HIST("hMassD0RecRef"), hfHelper.invMassD0ToPiK(candidate1), candidate1.pt(), efficiencyWeight); + } else { + registry.fill(HIST("hMassD0RecBg"), hfHelper.invMassD0ToPiK(candidate1), candidate1.pt(), efficiencyWeight); + } + } + if (candidate1.isSelD0bar() >= selectionFlagD0bar) { // only reco as D0bar + if (candidate1.flagMcMatchRec() == -(1 << aod::hf_cand_2prong::DecayType::D0ToPiK)) { // also matched as D0bar + registry.fill(HIST("hMassD0barRecSig"), hfHelper.invMassD0barToKPi(candidate1), candidate1.pt(), efficiencyWeight); + } else if (candidate1.flagMcMatchRec() == 1 << aod::hf_cand_2prong::DecayType::D0ToPiK) { + registry.fill(HIST("hMassD0barRecRef"), hfHelper.invMassD0barToKPi(candidate1), candidate1.pt(), efficiencyWeight); + } else { + registry.fill(HIST("hMassD0barRecBg"), hfHelper.invMassD0barToKPi(candidate1), candidate1.pt(), efficiencyWeight); + } + } + + // ===================== Define parameters for soft pion removal ======================== + auto ePiK = RecoDecay::e(candidate1.pVectorProng0(), massPi) + RecoDecay::e(candidate1.pVectorProng1(), massK); + auto eKPi = RecoDecay::e(candidate1.pVectorProng0(), massK) + RecoDecay::e(candidate1.pVectorProng1(), massPi); + + // ============== D-h correlation dedicated section ==================================== + + flagD0 = candidate1.flagMcMatchRec() == (1 << aod::hf_cand_2prong::DecayType::D0ToPiK); // flagD0Signal 'true' if candidate1 matched to D0 (particle) + flagD0bar = candidate1.flagMcMatchRec() == -(1 << aod::hf_cand_2prong::DecayType::D0ToPiK); // flagD0Reflection 'true' if candidate1, selected as D0 (particle), is matched to D0bar (antiparticle) + + // ========== track loop starts here ======================== + + for (const auto& track : tracks) { + registry.fill(HIST("hTrackCounterRec"), 1); // fill total no. of tracks + + // Removing D0 daughters by checking track indices + bool correlationStatus = false; + if ((candidate1.prong0Id() == track.globalIndex()) || (candidate1.prong1Id() == track.globalIndex())) { + if (!storeAutoCorrelationFlag) { + continue; + } + correlationStatus = true; + } + registry.fill(HIST("hTrackCounterRec"), 2); // fill no. of tracks before soft pion removal + + // ===== soft pion removal =================================================== + double invMassDstar1 = 0, invMassDstar2 = 0; + bool isSoftPiD0 = false, isSoftPiD0bar = false; + auto pSum2 = RecoDecay::p2(candidate1.pVector(), track.pVector()); + auto ePion = track.energy(massPi); + invMassDstar1 = std::sqrt((ePiK + ePion) * (ePiK + ePion) - pSum2); + invMassDstar2 = std::sqrt((eKPi + ePion) * (eKPi + ePion) - pSum2); + + if (candidate1.isSelD0() >= selectionFlagD0) { + if ((std::abs(invMassDstar1 - hfHelper.invMassD0ToPiK(candidate1)) - softPiMass) < ptSoftPionMax) { + isSoftPiD0 = true; + continue; + } + } + + if (candidate1.isSelD0bar() >= selectionFlagD0bar) { + if ((std::abs(invMassDstar2 - hfHelper.invMassD0barToKPi(candidate1)) - softPiMass) < ptSoftPionMax) { + isSoftPiD0bar = true; + continue; + } + } + + registry.fill(HIST("hTrackCounterRec"), 3); // fill no. of tracks after soft pion removal + + if (correlateD0WithLeadingParticle) { + if (track.globalIndex() != leadingIndex) { + continue; + } + registry.fill(HIST("hTrackCounter"), 4); // fill no. of tracks have leading particle + } + + int signalStatus = 0; + if (flagD0 && (candidate1.isSelD0() >= selectionFlagD0) && !isSoftPiD0) { + SETBIT(signalStatus, aod::hf_correlation_d0_hadron::ParticleTypeMcRec::D0Sig); + } // signal case D0 + if (flagD0bar && (candidate1.isSelD0() >= selectionFlagD0) && !isSoftPiD0) { + SETBIT(signalStatus, aod::hf_correlation_d0_hadron::ParticleTypeMcRec::D0Ref); + } // reflection case D0 + if (!flagD0 && !flagD0bar && (candidate1.isSelD0() >= selectionFlagD0) && !isSoftPiD0) { + SETBIT(signalStatus, aod::hf_correlation_d0_hadron::ParticleTypeMcRec::D0Bg); + } // background case D0 + + if (flagD0bar && (candidate1.isSelD0bar() >= selectionFlagD0bar) && !isSoftPiD0bar) { + SETBIT(signalStatus, aod::hf_correlation_d0_hadron::ParticleTypeMcRec::D0barSig); + } // signal case D0bar + if (flagD0 && (candidate1.isSelD0bar() >= selectionFlagD0bar) && !isSoftPiD0bar) { + SETBIT(signalStatus, aod::hf_correlation_d0_hadron::ParticleTypeMcRec::D0barRef); + } // reflection case D0bar + if (!flagD0 && !flagD0bar && (candidate1.isSelD0bar() >= selectionFlagD0bar) && !isSoftPiD0bar) { + SETBIT(signalStatus, aod::hf_correlation_d0_hadron::ParticleTypeMcRec::D0barBg); + } // background case D0bar + + entryD0HadronPair(getDeltaPhi(track.phi(), candidate1.phi()), + track.eta() - candidate1.eta(), + candidate1.pt(), + track.pt(), + poolBin, + correlationStatus); + entryD0HadronRecoInfo(hfHelper.invMassD0ToPiK(candidate1), hfHelper.invMassD0barToKPi(candidate1), signalStatus); + } // end inner loop (Tracks) + } // end of outer loop (D0) + registry.fill(HIST("hZvtx"), collision.posZ()); + registry.fill(HIST("hMultV0M"), collision.multFT0M()); + } + + PROCESS_SWITCH(HfCorrelatorD0Hadrons, processMcRec, "Process MC Reco mode", true); + + // ================= Process starts for MCGen, same event =================== + + void processMcGen(aod::McCollision const& mcCollision, + soa::Join const& mcParticles) + { + registry.fill(HIST("hEvtCountGen"), 0); + // MC gen level + for (const auto& particle1 : mcParticles) { + // check if the particle is D0 or D0bar (for general plot filling and selection, so both cases are fine) - NOTE: decay channel is not probed! + if (std::abs(particle1.pdgCode()) != Pdg::kD0) { + continue; + } + double yD = RecoDecay::y(particle1.pVector(), MassD0); + if (yCandMax >= 0. && std::abs(yD) > yCandMax) { + continue; + } + if (ptCandMin >= 0. && particle1.pt() < ptCandMin) { + continue; + } + registry.fill(HIST("hPtCandGen"), particle1.pt()); + registry.fill(HIST("hEtaGen"), particle1.eta()); + registry.fill(HIST("hPhiGen"), particle1.phi()); + registry.fill(HIST("hYGen"), yD); + + // =============== D-h correlation dedicated section ===================== + + if (std::abs(particle1.pdgCode()) != Pdg::kD0) { // just checking the particle PDG, not the decay channel (differently from Reco: you have a BR factor btw such levels!) + continue; + } + registry.fill(HIST("hCountD0TriggersGen"), 0, particle1.pt()); // to count trigger D0 (for normalisation) + + for (const auto& particle2 : mcParticles) { + registry.fill(HIST("hTrackCounterGen"), 1); // total no. of tracks + if (std::abs(particle2.eta()) > etaTrackMax) { + continue; + } + if (particle2.pt() < ptTrackMin) { + continue; + } + if ((std::abs(particle2.pdgCode()) != kElectron) && (std::abs(particle2.pdgCode()) != kMuonMinus) && (std::abs(particle2.pdgCode()) != kPiPlus) && (std::abs(particle2.pdgCode()) != kKPlus) && (std::abs(particle2.pdgCode()) != kProton)) { + continue; + } + + // ==============================soft pion removal================================ + registry.fill(HIST("hTrackCounterGen"), 2); // fill before soft pi removal + // method used: indexMother = -1 by default if the mother doesn't match with given PID of the mother. We find mother of pion if it is D* and mother of D0 if it is D*. If they are both positive and they both match each other, then it is detected as a soft pion + + auto indexMotherPi = RecoDecay::getMother(mcParticles, particle2, Pdg::kDStar, true, nullptr, 1); // last arguement 1 is written to consider immediate decay mother only + auto indexMotherD0 = RecoDecay::getMother(mcParticles, particle1, Pdg::kDStar, true, nullptr, 1); + if (std::abs(particle2.pdgCode()) == kPiPlus && indexMotherPi >= 0 && indexMotherD0 >= 0 && indexMotherPi == indexMotherD0) + continue; + + registry.fill(HIST("hTrackCounterGen"), 3); // fill after soft pion removal + + auto getTracksSize = [&mcParticles](aod::McCollision const& /*collision*/) { + int nTracks = 0; + for (const auto& track : mcParticles) { + if (track.isPhysicalPrimary() && std::abs(track.eta()) < 1.0) { + nTracks++; + } + } + return nTracks; + }; + using BinningTypeMcGen = FlexibleBinningPolicy, aod::mccollision::PosZ, decltype(getTracksSize)>; + BinningTypeMcGen corrBinningMcGen{{getTracksSize}, {zPoolBins, multPoolBinsMcGen}, true}; + int poolBin = corrBinningMcGen.getBin(std::make_tuple(mcCollision.posZ(), getTracksSize(mcCollision))); + + bool correlationStatus = false; + entryD0HadronPair(getDeltaPhi(particle2.phi(), particle1.phi()), + particle2.eta() - particle1.eta(), + particle1.pt(), + particle2.pt(), + poolBin, + correlationStatus); + entryD0HadronRecoInfo(massD0, massD0, 0); // dummy info + } // end inner loop (Tracks) + } // end outer loop (D0) + } + + PROCESS_SWITCH(HfCorrelatorD0Hadrons, processMcGen, "Process MC Gen mode", false); + + // ====================== Implement Event mixing on Data =================================== + + void processDataMixedEvent(SelectedCollisions const& collisions, + SelectedCandidatesData const& candidates, + SelectedTracks const& tracks) + { + for (const auto& collision : collisions) { + registry.fill(HIST("hMultFT0M"), collision.multFT0M()); + registry.fill(HIST("hZvtx"), collision.posZ()); + } + + auto tracksTuple = std::make_tuple(candidates, tracks); + Pair pairData{corrBinning, numberEventsMixed, -1, collisions, tracksTuple, &cache}; + + for (const auto& [c1, tracks1, c2, tracks2] : pairData) { + // LOGF(info, "Mixed event collisions: Index = (%d, %d), tracks Size: (%d, %d), Z Vertex: (%f, %f), Pool Bin: (%d, %d)", c1.globalIndex(), c2.globalIndex(), tracks1.size(), tracks2.size(), c1.posZ(), c2.posZ(), corrBinning.getBin(std::make_tuple(c1.posZ(), c1.multFT0M())),corrBinning.getBin(std::make_tuple(c2.posZ(), c2.multFT0M()))); // For debug + int poolBin = corrBinning.getBin(std::make_tuple(c2.posZ(), c2.multFT0M())); + for (const auto& [t1, t2] : o2::soa::combinations(o2::soa::CombinationsFullIndexPolicy(tracks1, tracks2))) { + + if (std::abs(hfHelper.yD0(t1)) >= yCandMax) { + continue; + } + + // soft pion removal, signal status 1,3 for D0 and 2,3 for D0bar (SoftPi removed), signal status 11,13 for D0 and 12,13 for D0bar (only SoftPi) + auto ePiK = RecoDecay::e(t1.pVectorProng0(), massPi) + RecoDecay::e(t1.pVectorProng1(), massK); + auto eKPi = RecoDecay::e(t1.pVectorProng0(), massK) + RecoDecay::e(t1.pVectorProng1(), massPi); + double invMassDstar1 = 0., invMassDstar2 = 0.; + bool isSoftPiD0 = false, isSoftPiD0bar = false; + auto pSum2 = RecoDecay::p2(t1.pVector(), t2.pVector()); + auto ePion = t2.energy(massPi); + invMassDstar1 = std::sqrt((ePiK + ePion) * (ePiK + ePion) - pSum2); + invMassDstar2 = std::sqrt((eKPi + ePion) * (eKPi + ePion) - pSum2); + + if (t1.isSelD0() >= selectionFlagD0) { + if ((std::abs(invMassDstar1 - hfHelper.invMassD0ToPiK(t1)) - softPiMass) < ptSoftPionMax) { + isSoftPiD0 = true; + } + } + if (t1.isSelD0bar() >= selectionFlagD0bar) { + if ((std::abs(invMassDstar2 - hfHelper.invMassD0barToKPi(t1)) - softPiMass) < ptSoftPionMax) { + isSoftPiD0bar = true; + } + } + + int signalStatus = 0; + if (t1.isSelD0() >= selectionFlagD0) { + if (!isSoftPiD0) { + signalStatus += aod::hf_correlation_d0_hadron::ParticleTypeData::D0Only; + } else { + signalStatus += aod::hf_correlation_d0_hadron::ParticleTypeData::D0OnlySoftPi; + } + } + if (t1.isSelD0bar() >= selectionFlagD0bar) { + if (!isSoftPiD0bar) { + signalStatus += aod::hf_correlation_d0_hadron::ParticleTypeData::D0barOnly; + } else { + signalStatus += aod::hf_correlation_d0_hadron::ParticleTypeData::D0barOnlySoftPi; + } + } + bool correlationStatus = false; + entryD0HadronPair(getDeltaPhi(t1.phi(), t2.phi()), t1.eta() - t2.eta(), t1.pt(), t2.pt(), poolBin, correlationStatus); + entryD0HadronRecoInfo(hfHelper.invMassD0ToPiK(t1), hfHelper.invMassD0barToKPi(t1), signalStatus); + } + } + } + PROCESS_SWITCH(HfCorrelatorD0Hadrons, processDataMixedEvent, "Process data mixed event", false); + + // ====================== Implement Event mixing on McRec =================================== + + void processMcRecMixedEvent(SelectedCollisions const& collisions, + SelectedCandidatesMcRec const& candidates, + SelectedTracksMcRec const& tracks) + { + auto tracksTuple = std::make_tuple(candidates, tracks); + Pair pairMcRec{corrBinning, numberEventsMixed, -1, collisions, tracksTuple, &cache}; + bool flagD0 = false; + bool flagD0bar = false; + for (const auto& [c1, tracks1, c2, tracks2] : pairMcRec) { + int poolBin = corrBinning.getBin(std::make_tuple(c2.posZ(), c2.multFT0M())); + + for (const auto& [t1, t2] : o2::soa::combinations(o2::soa::CombinationsFullIndexPolicy(tracks1, tracks2))) { + + if (std::abs(hfHelper.yD0(t1)) >= yCandMax) { + continue; + } + + // soft pion removal + auto ePiK = RecoDecay::e(t1.pVectorProng0(), massPi) + RecoDecay::e(t1.pVectorProng1(), massK); + auto eKPi = RecoDecay::e(t1.pVectorProng0(), massK) + RecoDecay::e(t1.pVectorProng1(), massPi); + double invMassDstar1 = 0., invMassDstar2 = 0.; + bool isSoftPiD0 = false, isSoftPiD0bar = false; + auto pSum2 = RecoDecay::p2(t1.pVector(), t2.pVector()); + auto ePion = t2.energy(massPi); + invMassDstar1 = std::sqrt((ePiK + ePion) * (ePiK + ePion) - pSum2); + invMassDstar2 = std::sqrt((eKPi + ePion) * (eKPi + ePion) - pSum2); + + if (t1.isSelD0() >= selectionFlagD0) { + if ((std::abs(invMassDstar1 - hfHelper.invMassD0ToPiK(t1)) - softPiMass) < ptSoftPionMax) { + isSoftPiD0 = true; + } + } + + if (t1.isSelD0bar() >= selectionFlagD0bar) { + if ((std::abs(invMassDstar2 - hfHelper.invMassD0barToKPi(t1)) - softPiMass) < ptSoftPionMax) { + isSoftPiD0bar = true; + } + } + + flagD0 = t1.flagMcMatchRec() == (1 << aod::hf_cand_2prong::DecayType::D0ToPiK); // flagD0Signal 'true' if candidate1 matched to D0 (particle) + flagD0bar = t1.flagMcMatchRec() == -(1 << aod::hf_cand_2prong::DecayType::D0ToPiK); // flagD0Reflection 'true' if candidate1, selected as D0 (particle), is matched to D0bar (antiparticle) + int signalStatus = 0; + + if (flagD0 && (t1.isSelD0() >= selectionFlagD0)) { + if (!isSoftPiD0) { + SETBIT(signalStatus, aod::hf_correlation_d0_hadron::ParticleTypeMcRec::D0Sig); // signalStatus += 1; + } else { + SETBIT(signalStatus, aod::hf_correlation_d0_hadron::ParticleTypeMcRec::SoftPi); // signalStatus += 64; + } + } // signal case D0 + + if (flagD0bar && (t1.isSelD0() >= selectionFlagD0)) { + if (!isSoftPiD0) { + SETBIT(signalStatus, aod::hf_correlation_d0_hadron::ParticleTypeMcRec::D0Ref); // signalStatus += 2; + } else { + SETBIT(signalStatus, aod::hf_correlation_d0_hadron::ParticleTypeMcRec::SoftPi); // signalStatus += 64; + } + } // reflection case D0 + + if (!flagD0 && !flagD0bar && (t1.isSelD0() >= selectionFlagD0)) { + if (!isSoftPiD0) { + SETBIT(signalStatus, aod::hf_correlation_d0_hadron::ParticleTypeMcRec::D0Bg); // signalStatus += 4; + } else { + SETBIT(signalStatus, aod::hf_correlation_d0_hadron::ParticleTypeMcRec::SoftPi); + } + } // background case D0 + + if (flagD0bar && (t1.isSelD0bar() >= selectionFlagD0bar)) { + if (!isSoftPiD0bar) { + SETBIT(signalStatus, aod::hf_correlation_d0_hadron::ParticleTypeMcRec::D0barSig); // signalStatus += 8; + } else { + SETBIT(signalStatus, aod::hf_correlation_d0_hadron::ParticleTypeMcRec::SoftPi); + } + } // signal case D0bar + + if (flagD0 && (t1.isSelD0bar() >= selectionFlagD0bar)) { + if (!isSoftPiD0bar) { + SETBIT(signalStatus, aod::hf_correlation_d0_hadron::ParticleTypeMcRec::D0barRef); // signalStatus += 16; + } else { + SETBIT(signalStatus, aod::hf_correlation_d0_hadron::ParticleTypeMcRec::SoftPi); + } + } // reflection case D0bar + + if (!flagD0 && !flagD0bar && (t1.isSelD0bar() >= selectionFlagD0bar)) { + if (!isSoftPiD0bar) { + SETBIT(signalStatus, aod::hf_correlation_d0_hadron::ParticleTypeMcRec::D0barBg); // signalStatus += 32; + } else { + SETBIT(signalStatus, aod::hf_correlation_d0_hadron::ParticleTypeMcRec::SoftPi); + } + } // background case D0bar + registry.fill(HIST("hSignalStatusMERec"), signalStatus); + bool correlationStatus = false; + entryD0HadronPair(getDeltaPhi(t1.phi(), t2.phi()), t1.eta() - t2.eta(), t1.pt(), t2.pt(), poolBin, correlationStatus); + entryD0HadronRecoInfo(hfHelper.invMassD0ToPiK(t1), hfHelper.invMassD0barToKPi(t1), signalStatus); + } + } + } + PROCESS_SWITCH(HfCorrelatorD0Hadrons, processMcRecMixedEvent, "Process Mixed Event MCRec", false); + + // ====================== Implement Event mixing on McGen =================================== + + void processMcGenMixedEvent(SelectedCollisionsMcGen const& collisions, + SelectedParticlesMcGen const& mcParticles) + { + + auto getTracksSize = [&mcParticles, this](SelectedCollisionsMcGen::iterator const& collision) { + int nTracks = 0; + auto associatedTracks = mcParticles.sliceByCached(o2::aod::mcparticle::mcCollisionId, collision.globalIndex(), this->cache); + for (const auto& track : associatedTracks) { + if (track.isPhysicalPrimary() && std::abs(track.eta()) < 1.0) { + nTracks++; + } + } + return nTracks; + }; + + using BinningTypeMcGen = FlexibleBinningPolicy, aod::mccollision::PosZ, decltype(getTracksSize)>; + BinningTypeMcGen corrBinningMcGen{{getTracksSize}, {zPoolBins, multPoolBinsMcGen}, true}; + + auto tracksTuple = std::make_tuple(mcParticles, mcParticles); + Pair pairMcGen{corrBinningMcGen, numberEventsMixed, -1, collisions, tracksTuple, &cache}; + + for (const auto& [c1, tracks1, c2, tracks2] : pairMcGen) { + for (const auto& [t1, t2] : o2::soa::combinations(o2::soa::CombinationsFullIndexPolicy(tracks1, tracks2))) { + + // Check track t1 is D0 + if (std::abs(t1.pdgCode()) != Pdg::kD0) { + continue; + } + + double yD = RecoDecay::y(t1.pVector(), MassD0); + if (std::abs(yD) >= yCandMax || t1.pt() <= ptCandMin || std::abs(t2.eta()) >= etaTrackMax || t2.pt() <= ptTrackMin) { + continue; + } + if ((std::abs(t2.pdgCode()) != kElectron) && (std::abs(t2.pdgCode()) != kMuonMinus) && (std::abs(t2.pdgCode()) != kPiPlus) && (std::abs(t2.pdgCode()) != kKPlus) && (std::abs(t2.pdgCode()) != kProton)) { + continue; + } + if (!t2.isPhysicalPrimary()) { + continue; + } + + // ==============================soft pion removal================================ + // method used: indexMother = -1 by default if the mother doesn't match with given PID of the mother. We find mother of pion if it is D* and mother of D0 if it is D*. If they are both positive and they both match each other, then it is detected as a soft pion + + auto indexMotherPi = RecoDecay::getMother(mcParticles, t2, Pdg::kDStar, true, nullptr, 1); // last arguement 1 is written to consider immediate decay mother only + auto indexMotherD0 = RecoDecay::getMother(mcParticles, t1, Pdg::kDStar, true, nullptr, 1); + if (std::abs(t2.pdgCode()) == kPiPlus && indexMotherPi >= 0 && indexMotherD0 >= 0 && indexMotherPi == indexMotherD0) { + continue; + } + int poolBin = corrBinningMcGen.getBin(std::make_tuple(c2.posZ(), getTracksSize(c2))); + bool correlationStatus = false; + entryD0HadronPair(getDeltaPhi(t2.phi(), t1.phi()), t2.eta() - t1.eta(), t1.pt(), t2.pt(), poolBin, correlationStatus); + entryD0HadronRecoInfo(massD0, massD0, 0); // dummy info + } + } + } + PROCESS_SWITCH(HfCorrelatorD0Hadrons, processMcGenMixedEvent, "Process Mixed Event MCGen", false); +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{ + adaptAnalysisTask(cfgc), + adaptAnalysisTask(cfgc)}; +} From bb03560bbdf480ac47c19ebacb262bed95002c0c Mon Sep 17 00:00:00 2001 From: samrangy Date: Fri, 29 Nov 2024 18:45:19 +0100 Subject: [PATCH 07/11] Implement ML models in data --- PWGHF/HFC/DataModel/CorrelationTables.h | 27 +- .../HFC/TableProducer/correlatorD0Hadrons.cxx | 296 +++++++++--------- PWGHF/HFC/Tasks/taskCorrelationD0Hadrons.cxx | 47 ++- 3 files changed, 221 insertions(+), 149 deletions(-) diff --git a/PWGHF/HFC/DataModel/CorrelationTables.h b/PWGHF/HFC/DataModel/CorrelationTables.h index 3c5b8cf86fd..cc121170f3e 100644 --- a/PWGHF/HFC/DataModel/CorrelationTables.h +++ b/PWGHF/HFC/DataModel/CorrelationTables.h @@ -57,6 +57,12 @@ DECLARE_SOA_COLUMN(PtD, ptD, float); //! Transverse mom DECLARE_SOA_COLUMN(PtHadron, ptHadron, float); //! Transverse momentum of Hadron DECLARE_SOA_COLUMN(MD, mD, float); //! Invariant mass of D0 DECLARE_SOA_COLUMN(MDbar, mDbar, float); //! Invariant mass of D0bar +DECLARE_SOA_COLUMN(MlScoreBkgD0, mlScoreBkgD0, float); //! ML background score for D0 selection +DECLARE_SOA_COLUMN(MlScoreNonPromptD0, mlScoreNonPromptD0, float); //! ML prompt score for D0 selection +DECLARE_SOA_COLUMN(MlScorePromptD0, mlScorePromptD0, float); //! ML prompt score for D0 selection +DECLARE_SOA_COLUMN(MlScoreBkgD0bar, mlScoreBkgD0bar, float); //! ML background score for D0 selection +DECLARE_SOA_COLUMN(MlScoreNonPromptD0bar, mlScoreNonPromptD0bar, float); //! ML prompt score for D0 selection +DECLARE_SOA_COLUMN(MlScorePromptD0bar, mlScorePromptD0bar, float); //! ML prompt score for D0 selection DECLARE_SOA_COLUMN(SignalStatus, signalStatus, int); //! Tag for D0,D0bar DECLARE_SOA_COLUMN(PoolBin, poolBin, int); //! Pool Bin for the MixedEvent DECLARE_SOA_COLUMN(IsAutoCorrelated, isAutoCorrelated, bool); //! Correlation Status @@ -81,7 +87,7 @@ enum ParticleTypeMcRec { }; } // namespace hf_correlation_d0_hadron -DECLARE_SOA_TABLE(DHadronPair, "AOD", "DHADRONPAIR", //! D0-Hadrons pairs Informations +DECLARE_SOA_TABLE(D0HadronPair, "AOD", "D0HPAIR", //! D0-Hadrons pairs Informations aod::hf_correlation_d0_hadron::DeltaPhi, aod::hf_correlation_d0_hadron::DeltaEta, aod::hf_correlation_d0_hadron::PtD, @@ -89,11 +95,28 @@ DECLARE_SOA_TABLE(DHadronPair, "AOD", "DHADRONPAIR", //! D0-Hadrons pairs Inform aod::hf_correlation_d0_hadron::PoolBin, aod::hf_correlation_d0_hadron::IsAutoCorrelated); -DECLARE_SOA_TABLE(DHadronRecoInfo, "AOD", "DHADRONRECOINFO", //! D0-Hadrons pairs Reconstructed Informations +DECLARE_SOA_TABLE(D0HadronRecoInfo, "AOD", "D0HRECOINFO", //! D0-Hadrons pairs Reconstructed Informations aod::hf_correlation_d0_hadron::MD, aod::hf_correlation_d0_hadron::MDbar, aod::hf_correlation_d0_hadron::SignalStatus); +DECLARE_SOA_TABLE(D0HadronMlInfo, "AOD", "D0HMLINFO", //! D0-Hadrons pairs Machine Learning Information + aod::hf_correlation_d0_hadron::MlScoreBkgD0, + aod::hf_correlation_d0_hadron::MlScoreNonPromptD0, + aod::hf_correlation_d0_hadron::MlScorePromptD0, + aod::hf_correlation_d0_hadron::MlScoreBkgD0bar, + aod::hf_correlation_d0_hadron::MlScoreNonPromptD0bar, + aod::hf_correlation_d0_hadron::MlScorePromptD0bar); + +DECLARE_SOA_TABLE(D0CandRecoInfo, "AOD", "D0CANDRECOINFO", //! Ds candidates Reconstructed Information + aod::hf_correlation_d0_hadron::MD, + aod::hf_correlation_d0_hadron::MDbar, + aod::hf_correlation_d0_hadron::PtD, + aod::hf_correlation_d0_hadron::MlScoreBkgD0, + aod::hf_correlation_d0_hadron::MlScorePromptD0, + aod::hf_correlation_d0_hadron::MlScoreBkgD0bar, + aod::hf_correlation_d0_hadron::MlScorePromptD0bar); + // Note: definition of columns and tables for Lc-Hadron correlation pairs namespace hf_correlation_lc_hadron { diff --git a/PWGHF/HFC/TableProducer/correlatorD0Hadrons.cxx b/PWGHF/HFC/TableProducer/correlatorD0Hadrons.cxx index 6ad93e81d42..206bf849388 100644 --- a/PWGHF/HFC/TableProducer/correlatorD0Hadrons.cxx +++ b/PWGHF/HFC/TableProducer/correlatorD0Hadrons.cxx @@ -15,6 +15,8 @@ /// \author Samrangy Sadhu , INFN Bari /// \author Swapnesh Santosh Khade , IIT Indore +#include + #include "CommonConstants/PhysicsConstants.h" #include "Framework/AnalysisTask.h" #include "Framework/HistogramRegistry.h" @@ -30,12 +32,14 @@ #include "PWGHF/DataModel/CandidateReconstructionTables.h" #include "PWGHF/DataModel/CandidateSelectionTables.h" #include "PWGHF/HFC/DataModel/CorrelationTables.h" +#include "PWGHF/HFC/Utils/utilsCorrelations.h" using namespace o2; using namespace o2::analysis; using namespace o2::constants::physics; using namespace o2::framework; using namespace o2::framework::expressions; +using namespace o2::analysis::hf_correlations; /// /// Returns deltaPhi value in range [-pi/2., 3.*pi/2], typically used for correlation studies @@ -68,6 +72,7 @@ using BinningType = ColumnBinningPolicy>; using SelectedTracks = soa::Filtered>; using SelectedCandidatesData = soa::Filtered>; +using SelectedCandidatesDataMl = soa::Filtered>; using SelectedTracksMcRec = soa::Filtered>; using SelectedCandidatesMcRec = soa::Filtered>; using SelectedCollisionsMcGen = soa::Filtered>; @@ -145,8 +150,8 @@ struct HfCorrelatorD0HadronsSelection { } double yD = RecoDecay::y(particle.pVector(), MassD0); if (std::abs(yD) > yCandMax || particle.pt() < ptCandMin) { - continue; - } + continue; + } isD0Found = 1; break; } @@ -159,8 +164,10 @@ struct HfCorrelatorD0HadronsSelection { struct HfCorrelatorD0Hadrons { SliceCache cache; - Produces entryD0HadronPair; - Produces entryD0HadronRecoInfo; + Produces entryD0HadronPair; + Produces entryD0HadronRecoInfo; + Produces entryD0HadronMlInfo; + Produces entryD0CandRecoInfo; Configurable selectionFlagD0{"selectionFlagD0", 1, "Selection Flag for D0"}; Configurable selectionFlagD0bar{"selectionFlagD0bar", 1, "Selection Flag for D0bar"}; @@ -176,7 +183,8 @@ struct HfCorrelatorD0Hadrons { Configurable applyEfficiency{"efficiencyFlagD", 1, "Flag for applying D-meson efficiency weights"}; Configurable multMin{"multMin", 0., "minimum multiplicity accepted"}; Configurable multMax{"multMax", 10000., "maximum multiplicity accepted"}; - Configurable ptSoftPionMax{"ptSoftPionMax", 3 * 800. * pow(10., -6.), "max. pT cut for soft pion identification"}; + Configurable ptSoftPionMax{"ptSoftPionMax", 3.f * 800.f * std::pow(10.f, -6.f), "max. pT cut for soft pion identification"}; + Configurable> classMl{"classMl", {0, 1, 2}, "Indexes of ML scores to be stored. Three indexes max."}; Configurable correlateD0WithLeadingParticle{"correlateD0WithLeadingParticle", false, "Switch for correlation of D0 mesons with leading particle only"}; Configurable storeAutoCorrelationFlag{"storeAutoCorrelationFlag", false, "Store flag that indicates if the track is paired to its D-meson mother instead of skipping it"}; Configurable numberEventsMixed{"numberEventsMixed", 5, "Number of events mixed in ME process"}; @@ -195,9 +203,6 @@ struct HfCorrelatorD0Hadrons { Preslice perCol = aod::hf_cand::collisionId; - Partition> selectedD0Candidates = aod::hf_sel_candidate_d0::isSelD0 >= selectionFlagD0 || aod::hf_sel_candidate_d0::isSelD0bar >= selectionFlagD0bar; - Partition> selectedD0candidatesMc = aod::hf_sel_candidate_d0::isSelD0 >= selectionFlagD0 || aod::hf_sel_candidate_d0::isSelD0bar >= selectionFlagD0bar; - Filter collisionFilter = aod::hf_selection_dmeson_collision::dmesonSel == true; Filter trackFilter = requireGlobalTrackWoDCAInFilter() && (nabs(aod::track::eta) < etaTrackMax) && (aod::track::pt > ptTrackMin) && (aod::track::pt < ptTrackMax) && (nabs(aod::track::dcaXY) < dcaXYTrackMax) && (nabs(aod::track::dcaZ) < dcaZTrackMax); Filter d0Filter = (aod::hf_sel_candidate_d0::isSelD0 >= 1) || (aod::hf_sel_candidate_d0::isSelD0bar >= 1); @@ -264,36 +269,16 @@ struct HfCorrelatorD0Hadrons { registry.add("hCountD0TriggersGen", "D0 trigger particles - MC gen;;N of trigger D0", {HistType::kTH2F, {{1, -0.5, 0.5}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); } - // Find Leading Particle - template - int findLeadingParticle(TTracks const& tracks) - { - auto leadingParticle = tracks.begin(); - for (auto const& track : tracks) { - if (std::abs(track.dcaXY()) >= 1. || std::abs(track.dcaZ()) >= 1.) { - continue; - } - if (track.pt() > leadingParticle.pt()) { - leadingParticle = track; - } - } - int leadingIndex = leadingParticle.globalIndex(); - return leadingIndex; - } // ======= Process starts for Data, Same event ============ /// D0-h correlation pair builder - for real data and data-like analysis (i.e. reco-level w/o matching request via MC truth) void processData(SelectedCollisions::iterator const& collision, SelectedTracks const& tracks, - soa::Join const&) + SelectedCandidatesDataMl const& candidates) { - // protection against empty tables to be sliced - if (selectedD0Candidates.size() == 0) { - return; - } // find leading particle if (correlateD0WithLeadingParticle) { - leadingIndex = findLeadingParticle(tracks); + leadingIndex = findLeadingParticle(tracks, dcaXYTrackMax.value, dcaZTrackMax.value, etaTrackMax.value); } int poolBin = corrBinning.getBin(std::make_tuple(collision.posZ(), collision.multFT0M())); @@ -316,46 +301,54 @@ struct HfCorrelatorD0Hadrons { return; } registry.fill(HIST("hMultiplicity"), nTracks); - - auto selectedD0CandidatesGrouped = selectedD0Candidates->sliceByCached(aod::hf_cand::collisionId, collision.globalIndex(), cache); - - for (const auto& candidate1 : selectedD0CandidatesGrouped) { - if (std::abs(hfHelper.yD0(candidate1)) >= yCandMax || candidate1.pt() <= ptCandMin || candidate1.pt() >= ptTrackMax) { - continue; - } - // check decay channel flag for candidate1 - if (!TESTBIT(candidate1.hfflag(), aod::hf_cand_2prong::DecayType::D0ToPiK)) { + std::vector outputMlD0 = {-1., -1., -1.}; + std::vector outputMlD0bar = {-1., -1., -1.}; + + for (const auto& candidate : candidates) { + if (std::abs(hfHelper.yD0(candidate)) >= yCandMax || candidate.pt() <= ptCandMin || candidate.pt() >= ptTrackMax) { + continue; + } + // check decay channel flag for candidate + if (!TESTBIT(candidate.hfflag(), aod::hf_cand_2prong::DecayType::D0ToPiK)) { continue; } // ========================== Define parameters for soft pion removal ================================ - auto ePiK = RecoDecay::e(candidate1.pVectorProng0(), massPi) + RecoDecay::e(candidate1.pVectorProng1(), massK); - auto eKPi = RecoDecay::e(candidate1.pVectorProng0(), massK) + RecoDecay::e(candidate1.pVectorProng1(), massPi); + auto ePiK = RecoDecay::e(candidate.pVectorProng0(), massPi) + RecoDecay::e(candidate.pVectorProng1(), massK); + auto eKPi = RecoDecay::e(candidate.pVectorProng0(), massK) + RecoDecay::e(candidate.pVectorProng1(), massPi); // ========================== trigger efficiency ================================ double efficiencyWeight = 1.; if (applyEfficiency) { - efficiencyWeight = 1. / efficiencyDmeson->at(o2::analysis::findBin(bins, candidate1.pt())); + efficiencyWeight = 1. / efficiencyDmeson->at(o2::analysis::findBin(bins, candidate.pt())); } // ========================== Fill mass histo ================================ - if (candidate1.isSelD0() >= selectionFlagD0) { - registry.fill(HIST("hMass"), hfHelper.invMassD0ToPiK(candidate1), candidate1.pt(), efficiencyWeight); - registry.fill(HIST("hMass1D"), hfHelper.invMassD0ToPiK(candidate1), efficiencyWeight); - registry.fill(HIST("hMassD01D"), hfHelper.invMassD0ToPiK(candidate1), efficiencyWeight); + if (candidate.isSelD0() >= selectionFlagD0) { + registry.fill(HIST("hMass"), hfHelper.invMassD0ToPiK(candidate), candidate.pt(), efficiencyWeight); + registry.fill(HIST("hMass1D"), hfHelper.invMassD0ToPiK(candidate), efficiencyWeight); + registry.fill(HIST("hMassD01D"), hfHelper.invMassD0ToPiK(candidate), efficiencyWeight); + for (unsigned int iclass = 0; iclass < classMl->size(); iclass++) { + outputMlD0[iclass] = candidate.mlProbD0()[classMl->at(iclass)]; + } } - if (candidate1.isSelD0bar() >= selectionFlagD0bar) { - registry.fill(HIST("hMass"), hfHelper.invMassD0barToKPi(candidate1), candidate1.pt(), efficiencyWeight); - registry.fill(HIST("hMass1D"), hfHelper.invMassD0barToKPi(candidate1), efficiencyWeight); - registry.fill(HIST("hMassD0bar1D"), hfHelper.invMassD0barToKPi(candidate1), efficiencyWeight); + if (candidate.isSelD0bar() >= selectionFlagD0bar) { + registry.fill(HIST("hMass"), hfHelper.invMassD0barToKPi(candidate), candidate.pt(), efficiencyWeight); + registry.fill(HIST("hMass1D"), hfHelper.invMassD0barToKPi(candidate), efficiencyWeight); + registry.fill(HIST("hMassD0bar1D"), hfHelper.invMassD0barToKPi(candidate), efficiencyWeight); + for (unsigned int iclass = 0; iclass < classMl->size(); iclass++) { + outputMlD0bar[iclass] = candidate.mlProbD0bar()[classMl->at(iclass)]; + } } + entryD0CandRecoInfo(hfHelper.invMassD0ToPiK(candidate), hfHelper.invMassD0barToKPi(candidate), candidate.pt(), outputMlD0[0], outputMlD0[2], outputMlD0bar[0], outputMlD0bar[2]); + // ========================== Fill general histos ================================ - registry.fill(HIST("hPtCand"), candidate1.pt()); - registry.fill(HIST("hPtProng0"), candidate1.ptProng0()); - registry.fill(HIST("hPtProng1"), candidate1.ptProng1()); - registry.fill(HIST("hEta"), candidate1.eta()); - registry.fill(HIST("hPhi"), candidate1.phi()); - registry.fill(HIST("hY"), hfHelper.yD0(candidate1)); - registry.fill(HIST("hSelectionStatus"), candidate1.isSelD0bar() + (candidate1.isSelD0() * 2)); + registry.fill(HIST("hPtCand"), candidate.pt()); + registry.fill(HIST("hPtProng0"), candidate.ptProng0()); + registry.fill(HIST("hPtProng1"), candidate.ptProng1()); + registry.fill(HIST("hEta"), candidate.eta()); + registry.fill(HIST("hPhi"), candidate.phi()); + registry.fill(HIST("hY"), hfHelper.yD0(candidate)); + registry.fill(HIST("hSelectionStatus"), candidate.isSelD0bar() + (candidate.isSelD0() * 2)); registry.fill(HIST("hD0Bin"), poolBin); // ============ D-h correlation dedicated section ================================== @@ -365,7 +358,7 @@ struct HfCorrelatorD0Hadrons { registry.fill(HIST("hTrackCounter"), 1); // fill total no. of tracks // Remove D0 daughters by checking track indices bool correlationStatus = false; - if ((candidate1.prong0Id() == track.globalIndex()) || (candidate1.prong1Id() == track.globalIndex())) { + if ((candidate.prong0Id() == track.globalIndex()) || (candidate.prong1Id() == track.globalIndex())) { if (!storeAutoCorrelationFlag) { continue; } @@ -377,20 +370,20 @@ struct HfCorrelatorD0Hadrons { // ========== soft pion removal =================================================== double invMassDstar1 = 0., invMassDstar2 = 0.; bool isSoftPiD0 = false, isSoftPiD0bar = false; - auto pSum2 = RecoDecay::p2(candidate1.pVector(), track.pVector()); + auto pSum2 = RecoDecay::p2(candidate.pVector(), track.pVector()); auto ePion = track.energy(massPi); invMassDstar1 = std::sqrt((ePiK + ePion) * (ePiK + ePion) - pSum2); invMassDstar2 = std::sqrt((eKPi + ePion) * (eKPi + ePion) - pSum2); - if (candidate1.isSelD0() >= selectionFlagD0) { - if ((std::abs(invMassDstar1 - hfHelper.invMassD0ToPiK(candidate1)) - softPiMass) < ptSoftPionMax) { + if (candidate.isSelD0() >= selectionFlagD0) { + if ((std::abs(invMassDstar1 - hfHelper.invMassD0ToPiK(candidate)) - softPiMass) < ptSoftPionMax) { isSoftPiD0 = true; continue; } } - if (candidate1.isSelD0bar() >= selectionFlagD0bar) { - if ((std::abs(invMassDstar2 - hfHelper.invMassD0barToKPi(candidate1)) - softPiMass) < ptSoftPionMax) { + if (candidate.isSelD0bar() >= selectionFlagD0bar) { + if ((std::abs(invMassDstar2 - hfHelper.invMassD0barToKPi(candidate)) - softPiMass) < ptSoftPionMax) { isSoftPiD0bar = true; continue; } @@ -398,10 +391,10 @@ struct HfCorrelatorD0Hadrons { registry.fill(HIST("hTrackCounter"), 3); // fill no. of tracks after soft pion removal int signalStatus = 0; - if ((candidate1.isSelD0() >= selectionFlagD0) && !isSoftPiD0) { + if ((candidate.isSelD0() >= selectionFlagD0) && !isSoftPiD0) { signalStatus += aod::hf_correlation_d0_hadron::ParticleTypeData::D0Only; } - if ((candidate1.isSelD0bar() >= selectionFlagD0bar) && !isSoftPiD0bar) { + if ((candidate.isSelD0bar() >= selectionFlagD0bar) && !isSoftPiD0bar) { signalStatus += aod::hf_correlation_d0_hadron::ParticleTypeData::D0barOnly; } @@ -411,13 +404,14 @@ struct HfCorrelatorD0Hadrons { } registry.fill(HIST("hTrackCounter"), 4); // fill no. of tracks have leading particle } - entryD0HadronPair(getDeltaPhi(track.phi(), candidate1.phi()), - track.eta() - candidate1.eta(), - candidate1.pt(), + entryD0HadronPair(getDeltaPhi(track.phi(), candidate.phi()), + track.eta() - candidate.eta(), + candidate.pt(), track.pt(), poolBin, correlationStatus); - entryD0HadronRecoInfo(hfHelper.invMassD0ToPiK(candidate1), hfHelper.invMassD0barToKPi(candidate1), signalStatus); + entryD0HadronRecoInfo(hfHelper.invMassD0ToPiK(candidate), hfHelper.invMassD0barToKPi(candidate), signalStatus); + entryD0HadronMlInfo(outputMlD0[0],outputMlD0[1],outputMlD0[2],outputMlD0bar[0],outputMlD0bar[1],outputMlD0bar[2]); } // end inner loop (tracks) @@ -429,15 +423,11 @@ struct HfCorrelatorD0Hadrons { void processMcRec(SelectedCollisions::iterator const& collision, SelectedTracksMcRec const& tracks, - soa::Join const&) + SelectedCandidatesMcRec const& candidates) { - // protection against empty tables to be sliced - if (selectedD0candidatesMc.size() == 0) { - return; - } // find leading particle if (correlateD0WithLeadingParticle) { - leadingIndex = findLeadingParticle(tracks); + leadingIndex = findLeadingParticle(tracks, dcaXYTrackMax.value, dcaZTrackMax.value, etaTrackMax.value); } int poolBin = corrBinning.getBin(std::make_tuple(collision.posZ(), collision.multFT0M())); int nTracks = 0; @@ -458,72 +448,71 @@ struct HfCorrelatorD0Hadrons { } registry.fill(HIST("hMultiplicity"), nTracks); - auto selectedD0CandidatesGroupedMc = selectedD0candidatesMc->sliceByCached(aod::hf_cand::collisionId, collision.globalIndex(), cache); // MC reco level bool flagD0 = false; bool flagD0bar = false; - for (const auto& candidate1 : selectedD0CandidatesGroupedMc) { - // check decay channel flag for candidate1 - if (!TESTBIT(candidate1.hfflag(), aod::hf_cand_2prong::DecayType::D0ToPiK)) { + for (const auto& candidate : candidates) { + // check decay channel flag for candidate + if (!TESTBIT(candidate.hfflag(), aod::hf_cand_2prong::DecayType::D0ToPiK)) { + continue; + } + if (std::abs(hfHelper.yD0(candidate)) >= yCandMax || candidate.pt() <= ptCandMin || candidate.pt() >= ptTrackMax) { continue; } - if (std::abs(hfHelper.yD0(candidate1)) >= yCandMax || candidate1.pt() <= ptCandMin || candidate1.pt() >= ptTrackMax) { - continue; - } double efficiencyWeight = 1.; if (applyEfficiency) { - efficiencyWeight = 1. / efficiencyDmeson->at(o2::analysis::findBin(bins, candidate1.pt())); + efficiencyWeight = 1. / efficiencyDmeson->at(o2::analysis::findBin(bins, candidate.pt())); } - if (std::abs(candidate1.flagMcMatchRec()) == 1 << aod::hf_cand_2prong::DecayType::D0ToPiK) { + if (std::abs(candidate.flagMcMatchRec()) == 1 << aod::hf_cand_2prong::DecayType::D0ToPiK) { // fill per-candidate distributions from D0/D0bar true candidates - registry.fill(HIST("hPtCandRec"), candidate1.pt()); - registry.fill(HIST("hPtProng0Rec"), candidate1.ptProng0()); - registry.fill(HIST("hPtProng1Rec"), candidate1.ptProng1()); - registry.fill(HIST("hEtaRec"), candidate1.eta()); - registry.fill(HIST("hPhiRec"), candidate1.phi()); - registry.fill(HIST("hYRec"), hfHelper.yD0(candidate1)); - registry.fill(HIST("hSelectionStatusRec"), candidate1.isSelD0bar() + (candidate1.isSelD0() * 2)); + registry.fill(HIST("hPtCandRec"), candidate.pt()); + registry.fill(HIST("hPtProng0Rec"), candidate.ptProng0()); + registry.fill(HIST("hPtProng1Rec"), candidate.ptProng1()); + registry.fill(HIST("hEtaRec"), candidate.eta()); + registry.fill(HIST("hPhiRec"), candidate.phi()); + registry.fill(HIST("hYRec"), hfHelper.yD0(candidate)); + registry.fill(HIST("hSelectionStatusRec"), candidate.isSelD0bar() + (candidate.isSelD0() * 2)); } // fill invariant mass plots from D0/D0bar signal and background candidates - if (candidate1.isSelD0() >= selectionFlagD0) { // only reco as D0 - if (candidate1.flagMcMatchRec() == 1 << aod::hf_cand_2prong::DecayType::D0ToPiK) { // also matched as D0 - registry.fill(HIST("hMassD0RecSig"), hfHelper.invMassD0ToPiK(candidate1), candidate1.pt(), efficiencyWeight); - } else if (candidate1.flagMcMatchRec() == -(1 << aod::hf_cand_2prong::DecayType::D0ToPiK)) { - registry.fill(HIST("hMassD0RecRef"), hfHelper.invMassD0ToPiK(candidate1), candidate1.pt(), efficiencyWeight); + if (candidate.isSelD0() >= selectionFlagD0) { // only reco as D0 + if (candidate.flagMcMatchRec() == 1 << aod::hf_cand_2prong::DecayType::D0ToPiK) { // also matched as D0 + registry.fill(HIST("hMassD0RecSig"), hfHelper.invMassD0ToPiK(candidate), candidate.pt(), efficiencyWeight); + } else if (candidate.flagMcMatchRec() == -(1 << aod::hf_cand_2prong::DecayType::D0ToPiK)) { + registry.fill(HIST("hMassD0RecRef"), hfHelper.invMassD0ToPiK(candidate), candidate.pt(), efficiencyWeight); } else { - registry.fill(HIST("hMassD0RecBg"), hfHelper.invMassD0ToPiK(candidate1), candidate1.pt(), efficiencyWeight); + registry.fill(HIST("hMassD0RecBg"), hfHelper.invMassD0ToPiK(candidate), candidate.pt(), efficiencyWeight); } } - if (candidate1.isSelD0bar() >= selectionFlagD0bar) { // only reco as D0bar - if (candidate1.flagMcMatchRec() == -(1 << aod::hf_cand_2prong::DecayType::D0ToPiK)) { // also matched as D0bar - registry.fill(HIST("hMassD0barRecSig"), hfHelper.invMassD0barToKPi(candidate1), candidate1.pt(), efficiencyWeight); - } else if (candidate1.flagMcMatchRec() == 1 << aod::hf_cand_2prong::DecayType::D0ToPiK) { - registry.fill(HIST("hMassD0barRecRef"), hfHelper.invMassD0barToKPi(candidate1), candidate1.pt(), efficiencyWeight); + if (candidate.isSelD0bar() >= selectionFlagD0bar) { // only reco as D0bar + if (candidate.flagMcMatchRec() == -(1 << aod::hf_cand_2prong::DecayType::D0ToPiK)) { // also matched as D0bar + registry.fill(HIST("hMassD0barRecSig"), hfHelper.invMassD0barToKPi(candidate), candidate.pt(), efficiencyWeight); + } else if (candidate.flagMcMatchRec() == 1 << aod::hf_cand_2prong::DecayType::D0ToPiK) { + registry.fill(HIST("hMassD0barRecRef"), hfHelper.invMassD0barToKPi(candidate), candidate.pt(), efficiencyWeight); } else { - registry.fill(HIST("hMassD0barRecBg"), hfHelper.invMassD0barToKPi(candidate1), candidate1.pt(), efficiencyWeight); + registry.fill(HIST("hMassD0barRecBg"), hfHelper.invMassD0barToKPi(candidate), candidate.pt(), efficiencyWeight); } } // ===================== Define parameters for soft pion removal ======================== - auto ePiK = RecoDecay::e(candidate1.pVectorProng0(), massPi) + RecoDecay::e(candidate1.pVectorProng1(), massK); - auto eKPi = RecoDecay::e(candidate1.pVectorProng0(), massK) + RecoDecay::e(candidate1.pVectorProng1(), massPi); + auto ePiK = RecoDecay::e(candidate.pVectorProng0(), massPi) + RecoDecay::e(candidate.pVectorProng1(), massK); + auto eKPi = RecoDecay::e(candidate.pVectorProng0(), massK) + RecoDecay::e(candidate.pVectorProng1(), massPi); // ============== D-h correlation dedicated section ==================================== - flagD0 = candidate1.flagMcMatchRec() == (1 << aod::hf_cand_2prong::DecayType::D0ToPiK); // flagD0Signal 'true' if candidate1 matched to D0 (particle) - flagD0bar = candidate1.flagMcMatchRec() == -(1 << aod::hf_cand_2prong::DecayType::D0ToPiK); // flagD0Reflection 'true' if candidate1, selected as D0 (particle), is matched to D0bar (antiparticle) + flagD0 = candidate.flagMcMatchRec() == (1 << aod::hf_cand_2prong::DecayType::D0ToPiK); // flagD0Signal 'true' if candidate matched to D0 (particle) + flagD0bar = candidate.flagMcMatchRec() == -(1 << aod::hf_cand_2prong::DecayType::D0ToPiK); // flagD0Reflection 'true' if candidate, selected as D0 (particle), is matched to D0bar (antiparticle) // ========== track loop starts here ======================== for (const auto& track : tracks) { registry.fill(HIST("hTrackCounterRec"), 1); // fill total no. of tracks - + // Removing D0 daughters by checking track indices bool correlationStatus = false; - if ((candidate1.prong0Id() == track.globalIndex()) || (candidate1.prong1Id() == track.globalIndex())) { + if ((candidate.prong0Id() == track.globalIndex()) || (candidate.prong1Id() == track.globalIndex())) { if (!storeAutoCorrelationFlag) { continue; } @@ -534,20 +523,20 @@ struct HfCorrelatorD0Hadrons { // ===== soft pion removal =================================================== double invMassDstar1 = 0, invMassDstar2 = 0; bool isSoftPiD0 = false, isSoftPiD0bar = false; - auto pSum2 = RecoDecay::p2(candidate1.pVector(), track.pVector()); + auto pSum2 = RecoDecay::p2(candidate.pVector(), track.pVector()); auto ePion = track.energy(massPi); invMassDstar1 = std::sqrt((ePiK + ePion) * (ePiK + ePion) - pSum2); invMassDstar2 = std::sqrt((eKPi + ePion) * (eKPi + ePion) - pSum2); - if (candidate1.isSelD0() >= selectionFlagD0) { - if ((std::abs(invMassDstar1 - hfHelper.invMassD0ToPiK(candidate1)) - softPiMass) < ptSoftPionMax) { + if (candidate.isSelD0() >= selectionFlagD0) { + if ((std::abs(invMassDstar1 - hfHelper.invMassD0ToPiK(candidate)) - softPiMass) < ptSoftPionMax) { isSoftPiD0 = true; continue; } } - if (candidate1.isSelD0bar() >= selectionFlagD0bar) { - if ((std::abs(invMassDstar2 - hfHelper.invMassD0barToKPi(candidate1)) - softPiMass) < ptSoftPionMax) { + if (candidate.isSelD0bar() >= selectionFlagD0bar) { + if ((std::abs(invMassDstar2 - hfHelper.invMassD0barToKPi(candidate)) - softPiMass) < ptSoftPionMax) { isSoftPiD0bar = true; continue; } @@ -559,37 +548,37 @@ struct HfCorrelatorD0Hadrons { if (track.globalIndex() != leadingIndex) { continue; } - registry.fill(HIST("hTrackCounter"), 4); // fill no. of tracks have leading particle + registry.fill(HIST("hTrackCounterRec"), 4); // fill no. of tracks have leading particle } int signalStatus = 0; - if (flagD0 && (candidate1.isSelD0() >= selectionFlagD0) && !isSoftPiD0) { + if (flagD0 && (candidate.isSelD0() >= selectionFlagD0) && !isSoftPiD0) { SETBIT(signalStatus, aod::hf_correlation_d0_hadron::ParticleTypeMcRec::D0Sig); } // signal case D0 - if (flagD0bar && (candidate1.isSelD0() >= selectionFlagD0) && !isSoftPiD0) { + if (flagD0bar && (candidate.isSelD0() >= selectionFlagD0) && !isSoftPiD0) { SETBIT(signalStatus, aod::hf_correlation_d0_hadron::ParticleTypeMcRec::D0Ref); } // reflection case D0 - if (!flagD0 && !flagD0bar && (candidate1.isSelD0() >= selectionFlagD0) && !isSoftPiD0) { + if (!flagD0 && !flagD0bar && (candidate.isSelD0() >= selectionFlagD0) && !isSoftPiD0) { SETBIT(signalStatus, aod::hf_correlation_d0_hadron::ParticleTypeMcRec::D0Bg); } // background case D0 - if (flagD0bar && (candidate1.isSelD0bar() >= selectionFlagD0bar) && !isSoftPiD0bar) { + if (flagD0bar && (candidate.isSelD0bar() >= selectionFlagD0bar) && !isSoftPiD0bar) { SETBIT(signalStatus, aod::hf_correlation_d0_hadron::ParticleTypeMcRec::D0barSig); } // signal case D0bar - if (flagD0 && (candidate1.isSelD0bar() >= selectionFlagD0bar) && !isSoftPiD0bar) { + if (flagD0 && (candidate.isSelD0bar() >= selectionFlagD0bar) && !isSoftPiD0bar) { SETBIT(signalStatus, aod::hf_correlation_d0_hadron::ParticleTypeMcRec::D0barRef); } // reflection case D0bar - if (!flagD0 && !flagD0bar && (candidate1.isSelD0bar() >= selectionFlagD0bar) && !isSoftPiD0bar) { + if (!flagD0 && !flagD0bar && (candidate.isSelD0bar() >= selectionFlagD0bar) && !isSoftPiD0bar) { SETBIT(signalStatus, aod::hf_correlation_d0_hadron::ParticleTypeMcRec::D0barBg); } // background case D0bar - entryD0HadronPair(getDeltaPhi(track.phi(), candidate1.phi()), - track.eta() - candidate1.eta(), - candidate1.pt(), + entryD0HadronPair(getDeltaPhi(track.phi(), candidate.phi()), + track.eta() - candidate.eta(), + candidate.pt(), track.pt(), poolBin, correlationStatus); - entryD0HadronRecoInfo(hfHelper.invMassD0ToPiK(candidate1), hfHelper.invMassD0barToKPi(candidate1), signalStatus); + entryD0HadronRecoInfo(hfHelper.invMassD0ToPiK(candidate), hfHelper.invMassD0barToKPi(candidate), signalStatus); } // end inner loop (Tracks) } // end of outer loop (D0) registry.fill(HIST("hZvtx"), collision.posZ()); @@ -605,6 +594,10 @@ struct HfCorrelatorD0Hadrons { { registry.fill(HIST("hEvtCountGen"), 0); // MC gen level + // find leading particle + if (correlateD0WithLeadingParticle) { + leadingIndex = findLeadingParticleMcGen(mcParticles, etaTrackMax.value, ptTrackMin.value); + } for (const auto& particle1 : mcParticles) { // check if the particle is D0 or D0bar (for general plot filling and selection, so both cases are fine) - NOTE: decay channel is not probed! if (std::abs(particle1.pdgCode()) != Pdg::kD0) { @@ -647,11 +640,23 @@ struct HfCorrelatorD0Hadrons { auto indexMotherPi = RecoDecay::getMother(mcParticles, particle2, Pdg::kDStar, true, nullptr, 1); // last arguement 1 is written to consider immediate decay mother only auto indexMotherD0 = RecoDecay::getMother(mcParticles, particle1, Pdg::kDStar, true, nullptr, 1); - if (std::abs(particle2.pdgCode()) == kPiPlus && indexMotherPi >= 0 && indexMotherD0 >= 0 && indexMotherPi == indexMotherD0) - continue; + bool correlationStatus = false; + if (std::abs(particle2.pdgCode()) == kPiPlus && indexMotherPi >= 0 && indexMotherD0 >= 0 && indexMotherPi == indexMotherD0) { + if (!storeAutoCorrelationFlag) { + continue; + } + correlationStatus = true; + } registry.fill(HIST("hTrackCounterGen"), 3); // fill after soft pion removal + if (correlateD0WithLeadingParticle) { + if (particle2.globalIndex() != leadingIndex) { + continue; + } + registry.fill(HIST("hTrackCounterGen"), 4); // fill no. of tracks have leading particle + } + auto getTracksSize = [&mcParticles](aod::McCollision const& /*collision*/) { int nTracks = 0; for (const auto& track : mcParticles) { @@ -664,8 +669,6 @@ struct HfCorrelatorD0Hadrons { using BinningTypeMcGen = FlexibleBinningPolicy, aod::mccollision::PosZ, decltype(getTracksSize)>; BinningTypeMcGen corrBinningMcGen{{getTracksSize}, {zPoolBins, multPoolBinsMcGen}, true}; int poolBin = corrBinningMcGen.getBin(std::make_tuple(mcCollision.posZ(), getTracksSize(mcCollision))); - - bool correlationStatus = false; entryD0HadronPair(getDeltaPhi(particle2.phi(), particle1.phi()), particle2.eta() - particle1.eta(), particle1.pt(), @@ -682,23 +685,23 @@ struct HfCorrelatorD0Hadrons { // ====================== Implement Event mixing on Data =================================== void processDataMixedEvent(SelectedCollisions const& collisions, - SelectedCandidatesData const& candidates, + SelectedCandidatesDataMl const& candidates, SelectedTracks const& tracks) { for (const auto& collision : collisions) { registry.fill(HIST("hMultFT0M"), collision.multFT0M()); registry.fill(HIST("hZvtx"), collision.posZ()); } - + auto tracksTuple = std::make_tuple(candidates, tracks); - Pair pairData{corrBinning, numberEventsMixed, -1, collisions, tracksTuple, &cache}; + Pair pairData{corrBinning, numberEventsMixed, -1, collisions, tracksTuple, &cache}; for (const auto& [c1, tracks1, c2, tracks2] : pairData) { // LOGF(info, "Mixed event collisions: Index = (%d, %d), tracks Size: (%d, %d), Z Vertex: (%f, %f), Pool Bin: (%d, %d)", c1.globalIndex(), c2.globalIndex(), tracks1.size(), tracks2.size(), c1.posZ(), c2.posZ(), corrBinning.getBin(std::make_tuple(c1.posZ(), c1.multFT0M())),corrBinning.getBin(std::make_tuple(c2.posZ(), c2.multFT0M()))); // For debug int poolBin = corrBinning.getBin(std::make_tuple(c2.posZ(), c2.multFT0M())); for (const auto& [t1, t2] : o2::soa::combinations(o2::soa::CombinationsFullIndexPolicy(tracks1, tracks2))) { - if (std::abs(hfHelper.yD0(t1)) >= yCandMax) { + if (std::abs(hfHelper.yD0(t1)) >= yCandMax) { continue; } @@ -711,16 +714,24 @@ struct HfCorrelatorD0Hadrons { auto ePion = t2.energy(massPi); invMassDstar1 = std::sqrt((ePiK + ePion) * (ePiK + ePion) - pSum2); invMassDstar2 = std::sqrt((eKPi + ePion) * (eKPi + ePion) - pSum2); - + std::vector outputMlD0 = {-1., -1., -1.}; + std::vector outputMlD0bar = {-1., -1., -1.}; + if (t1.isSelD0() >= selectionFlagD0) { if ((std::abs(invMassDstar1 - hfHelper.invMassD0ToPiK(t1)) - softPiMass) < ptSoftPionMax) { isSoftPiD0 = true; } + for (unsigned int iclass = 0; iclass < classMl->size(); iclass++) { + outputMlD0[iclass] = t1.mlProbD0()[classMl->at(iclass)]; + } } if (t1.isSelD0bar() >= selectionFlagD0bar) { if ((std::abs(invMassDstar2 - hfHelper.invMassD0barToKPi(t1)) - softPiMass) < ptSoftPionMax) { isSoftPiD0bar = true; } + for (unsigned int iclass = 0; iclass < classMl->size(); iclass++) { + outputMlD0bar[iclass] = t1.mlProbD0bar()[classMl->at(iclass)]; + } } int signalStatus = 0; @@ -741,6 +752,7 @@ struct HfCorrelatorD0Hadrons { bool correlationStatus = false; entryD0HadronPair(getDeltaPhi(t1.phi(), t2.phi()), t1.eta() - t2.eta(), t1.pt(), t2.pt(), poolBin, correlationStatus); entryD0HadronRecoInfo(hfHelper.invMassD0ToPiK(t1), hfHelper.invMassD0barToKPi(t1), signalStatus); + entryD0HadronMlInfo(outputMlD0[0],outputMlD0[1],outputMlD0[2],outputMlD0bar[0],outputMlD0bar[1],outputMlD0bar[2]); } } } @@ -787,8 +799,8 @@ struct HfCorrelatorD0Hadrons { } } - flagD0 = t1.flagMcMatchRec() == (1 << aod::hf_cand_2prong::DecayType::D0ToPiK); // flagD0Signal 'true' if candidate1 matched to D0 (particle) - flagD0bar = t1.flagMcMatchRec() == -(1 << aod::hf_cand_2prong::DecayType::D0ToPiK); // flagD0Reflection 'true' if candidate1, selected as D0 (particle), is matched to D0bar (antiparticle) + flagD0 = t1.flagMcMatchRec() == (1 << aod::hf_cand_2prong::DecayType::D0ToPiK); // flagD0Signal 'true' if candidate matched to D0 (particle) + flagD0bar = t1.flagMcMatchRec() == -(1 << aod::hf_cand_2prong::DecayType::D0ToPiK); // flagD0Reflection 'true' if candidate, selected as D0 (particle), is matched to D0bar (antiparticle) int signalStatus = 0; if (flagD0 && (t1.isSelD0() >= selectionFlagD0)) { @@ -883,12 +895,12 @@ struct HfCorrelatorD0Hadrons { continue; } if ((std::abs(t2.pdgCode()) != kElectron) && (std::abs(t2.pdgCode()) != kMuonMinus) && (std::abs(t2.pdgCode()) != kPiPlus) && (std::abs(t2.pdgCode()) != kKPlus) && (std::abs(t2.pdgCode()) != kProton)) { - continue; - } - if (!t2.isPhysicalPrimary()) { - continue; - } - + continue; + } + if (!t2.isPhysicalPrimary()) { + continue; + } + // ==============================soft pion removal================================ // method used: indexMother = -1 by default if the mother doesn't match with given PID of the mother. We find mother of pion if it is D* and mother of D0 if it is D*. If they are both positive and they both match each other, then it is detected as a soft pion diff --git a/PWGHF/HFC/Tasks/taskCorrelationD0Hadrons.cxx b/PWGHF/HFC/Tasks/taskCorrelationD0Hadrons.cxx index 416460e4689..bc99470cb54 100644 --- a/PWGHF/HFC/Tasks/taskCorrelationD0Hadrons.cxx +++ b/PWGHF/HFC/Tasks/taskCorrelationD0Hadrons.cxx @@ -36,7 +36,8 @@ using namespace o2::analysis::hf_correlations; namespace o2::aod { -using DHadronPairFull = soa::Join; +using D0HadronPairFull = soa::Join; +using D0HadronPairFullMl = soa::Join; } // namespace o2::aod // string definitions, used for histogram axis labels @@ -58,6 +59,8 @@ AxisSpec axisPtHadron = {11, 0., 11., ""}; AxisSpec axisPoolBin = {9, 0., 9., ""}; AxisSpec axisCorrelationState = {2, 0., 2., ""}; ConfigurableAxis axisMass{"axisMass", {250, 1.65f, 2.15f}, ""}; +ConfigurableAxis binsBdtScore{"binsBdtScore", {100, 0., 1.}, "Bdt output scores"}; +AxisSpec axisBdtScore = {binsBdtScore, "Bdt score"}; // definition of vectors for standard ptbin and invariant mass configurables const int nPtBinsCorrelations = 12; @@ -86,6 +89,11 @@ struct HfTaskCorrelationD0Hadrons { Configurable> binsCorrelations{"ptBinsForCorrelations", std::vector{vecPtBinsCorrelations}, "pT bin limits for correlation plots"}; // pT bins for effiencies: same as above Configurable> binsEfficiency{"ptBinsForEfficiency", std::vector{o2::analysis::hf_cuts_d0_to_pi_k::vecBinsPt}, "pT bin limits for efficiency"}; + Configurable> classMl{"classMl", {0, 1, 2}, "Indexes of ML scores to be stored. Three indexes max."}; + Configurable> mlOutputPromptD0{"mlScorePromptD0", {0.5, 0.5, 0.5, 0.5}, "Machine learning scores for prompt"}; + Configurable> mlOutputBkgD0{"mlScoreBkgD0", {0.5, 0.5, 0.5, 0.5}, "Machine learning scores for bkg"}; + Configurable> mlOutputPromptD0bar{"mlScorePromptD0bar", {0.5, 0.5, 0.5, 0.5}, "Machine learning scores for prompt"}; + Configurable> mlOutputBkgD0bar{"mlScoreBkgD0bar", {0.5, 0.5, 0.5, 0.5}, "Machine learning scores for bkg"}; // signal and sideband region edges, to be defined via json file (initialised to empty) Configurable> signalRegionLeft{"signalRegionLeft", std::vector{vecsignalRegionLeft}, "Inner values of signal region vs pT"}; Configurable> signalRegionRight{"signalRegionRight", std::vector{vecsignalRegionRight}, "Outer values of signal region vs pT"}; @@ -191,6 +199,10 @@ struct HfTaskCorrelationD0Hadrons { { int nBinsPtAxis = binsCorrelations->size() - 1; const double* valuesPtAxis = binsCorrelations->data(); + registry.add("hBdtScorePromptD0", "D0 BDT prompt score", {HistType::kTH1F, {axisBdtScore}}); + registry.add("hBdtScoreBkgD0", "D0 BDT bkg score", {HistType::kTH1F, {axisBdtScore}}); + registry.add("hBdtScorePromptD0bar", "D0bar BDT prompt score", {HistType::kTH1F, {axisBdtScore}}); + registry.add("hBdtScoreBkgD0bar", "D0bar BDT bkg score", {HistType::kTH1F, {axisBdtScore}}); registry.get(HIST("hCorrel2DVsPtSignalRegion"))->GetAxis(2)->Set(nBinsPtAxis, valuesPtAxis); registry.get(HIST("hCorrel2DVsPtSidebands"))->GetAxis(2)->Set(nBinsPtAxis, valuesPtAxis); registry.get(HIST("hCorrel2DVsPtSignalRegion"))->Sumw2(); @@ -233,8 +245,26 @@ struct HfTaskCorrelationD0Hadrons { /// D-h correlation pair filling task, from pair tables - for real data and data-like analysis (i.e. reco-level w/o matching request via MC truth) /// Works on both USL and LS analyses pair tables - void processData(aod::DHadronPairFull const& pairEntries) + void processData(aod::D0HadronPairFullMl const& pairEntries, + aod::D0CandRecoInfo const& candidates) { + for (const auto& candidate : candidates) { + float ptD = candidate.ptD(); + float bdtScorePromptD0 = candidate.mlScorePromptD0(); + float bdtScoreBkgD0 = candidate.mlScoreBkgD0(); + float bdtScorePromptD0bar = candidate.mlScorePromptD0bar(); + float bdtScoreBkgD0bar = candidate.mlScoreBkgD0bar(); + int effBinD = o2::analysis::findBin(binsEfficiency, ptD); + if (bdtScorePromptD0 < mlOutputPromptD0->at(effBinD) || bdtScoreBkgD0 > mlOutputBkgD0->at(effBinD) || + bdtScorePromptD0bar < mlOutputPromptD0bar->at(effBinD) || bdtScoreBkgD0bar > mlOutputBkgD0bar->at(effBinD)) { + continue; + } + registry.fill(HIST("hBdtScorePromptD0"), bdtScorePromptD0); + registry.fill(HIST("hBdtScoreBkgD0"), bdtScoreBkgD0); + registry.fill(HIST("hBdtScorePromptD0bar"), bdtScorePromptD0bar); + registry.fill(HIST("hBdtScoreBkgD0bar"), bdtScoreBkgD0bar); + } + for (const auto& pairEntry : pairEntries) { // define variables for widely used quantities double deltaPhi = pairEntry.deltaPhi(); @@ -248,6 +278,10 @@ struct HfTaskCorrelationD0Hadrons { int ptBinD = o2::analysis::findBin(binsCorrelations, ptD); int poolBin = pairEntry.poolBin(); bool isAutoCorrelated = pairEntry.isAutoCorrelated(); + float bdtScorePromptD0 = pairEntry.mlScorePromptD0(); + float bdtScoreBkgD0 = pairEntry.mlScoreBkgD0(); + float bdtScorePromptD0bar = pairEntry.mlScorePromptD0bar(); + float bdtScoreBkgD0bar = pairEntry.mlScoreBkgD0bar(); // reject entries outside pT ranges of interest if (ptBinD < 0 || effBinD < 0) { continue; @@ -255,7 +289,10 @@ struct HfTaskCorrelationD0Hadrons { if (ptHadron > ptHadronMax) { ptHadron = ptHadronMax + 0.5; } - + if (bdtScorePromptD0 < mlOutputPromptD0->at(ptBinD) || bdtScoreBkgD0 > mlOutputBkgD0->at(ptBinD) || + bdtScorePromptD0bar < mlOutputPromptD0bar->at(ptBinD) || bdtScoreBkgD0bar > mlOutputBkgD0bar->at(ptBinD)) { + continue; + } double efficiencyWeight = 1.; if (applyEfficiency) { efficiencyWeight = 1. / (efficiencyDmeson->at(o2::analysis::findBin(binsEfficiency, ptD))); // ***** track efficiency to be implemented ***** @@ -362,7 +399,7 @@ struct HfTaskCorrelationD0Hadrons { } PROCESS_SWITCH(HfTaskCorrelationD0Hadrons, processData, "Process data", false); - void processMcRec(aod::DHadronPairFull const& pairEntries) + void processMcRec(aod::D0HadronPairFull const& pairEntries) { for (const auto& pairEntry : pairEntries) { // define variables for widely used quantities @@ -580,7 +617,7 @@ struct HfTaskCorrelationD0Hadrons { PROCESS_SWITCH(HfTaskCorrelationD0Hadrons, processMcRec, "Process MC Reco mode", true); /// D-Hadron correlation pair filling task, from pair tables - for MC gen-level analysis (no filter/selection, only true signal) - void processMcGen(aod::DHadronPairFull const& pairEntries) + void processMcGen(aod::D0HadronPairFull const& pairEntries) { for (const auto& pairEntry : pairEntries) { // define variables for widely used quantities From b042d2593c56322e150a6546241d977b3cfb5247 Mon Sep 17 00:00:00 2001 From: ALICE Action Bot Date: Fri, 6 Dec 2024 20:36:06 +0000 Subject: [PATCH 08/11] Please consider the following formatting changes --- PWGHF/HFC/DataModel/CorrelationTables.h | 12 ++++++------ .../HFC/TableProducer/correlatorD0Hadrons.cxx | 18 +++++++++--------- PWGHF/HFC/Tasks/taskCorrelationD0Hadrons.cxx | 4 ++-- 3 files changed, 17 insertions(+), 17 deletions(-) diff --git a/PWGHF/HFC/DataModel/CorrelationTables.h b/PWGHF/HFC/DataModel/CorrelationTables.h index cc121170f3e..3b8e68cdbd8 100644 --- a/PWGHF/HFC/DataModel/CorrelationTables.h +++ b/PWGHF/HFC/DataModel/CorrelationTables.h @@ -57,12 +57,12 @@ DECLARE_SOA_COLUMN(PtD, ptD, float); //! Transverse mom DECLARE_SOA_COLUMN(PtHadron, ptHadron, float); //! Transverse momentum of Hadron DECLARE_SOA_COLUMN(MD, mD, float); //! Invariant mass of D0 DECLARE_SOA_COLUMN(MDbar, mDbar, float); //! Invariant mass of D0bar -DECLARE_SOA_COLUMN(MlScoreBkgD0, mlScoreBkgD0, float); //! ML background score for D0 selection -DECLARE_SOA_COLUMN(MlScoreNonPromptD0, mlScoreNonPromptD0, float); //! ML prompt score for D0 selection -DECLARE_SOA_COLUMN(MlScorePromptD0, mlScorePromptD0, float); //! ML prompt score for D0 selection -DECLARE_SOA_COLUMN(MlScoreBkgD0bar, mlScoreBkgD0bar, float); //! ML background score for D0 selection -DECLARE_SOA_COLUMN(MlScoreNonPromptD0bar, mlScoreNonPromptD0bar, float); //! ML prompt score for D0 selection -DECLARE_SOA_COLUMN(MlScorePromptD0bar, mlScorePromptD0bar, float); //! ML prompt score for D0 selection +DECLARE_SOA_COLUMN(MlScoreBkgD0, mlScoreBkgD0, float); //! ML background score for D0 selection +DECLARE_SOA_COLUMN(MlScoreNonPromptD0, mlScoreNonPromptD0, float); //! ML prompt score for D0 selection +DECLARE_SOA_COLUMN(MlScorePromptD0, mlScorePromptD0, float); //! ML prompt score for D0 selection +DECLARE_SOA_COLUMN(MlScoreBkgD0bar, mlScoreBkgD0bar, float); //! ML background score for D0 selection +DECLARE_SOA_COLUMN(MlScoreNonPromptD0bar, mlScoreNonPromptD0bar, float); //! ML prompt score for D0 selection +DECLARE_SOA_COLUMN(MlScorePromptD0bar, mlScorePromptD0bar, float); //! ML prompt score for D0 selection DECLARE_SOA_COLUMN(SignalStatus, signalStatus, int); //! Tag for D0,D0bar DECLARE_SOA_COLUMN(PoolBin, poolBin, int); //! Pool Bin for the MixedEvent DECLARE_SOA_COLUMN(IsAutoCorrelated, isAutoCorrelated, bool); //! Correlation Status diff --git a/PWGHF/HFC/TableProducer/correlatorD0Hadrons.cxx b/PWGHF/HFC/TableProducer/correlatorD0Hadrons.cxx index 206bf849388..c4043c8ed53 100644 --- a/PWGHF/HFC/TableProducer/correlatorD0Hadrons.cxx +++ b/PWGHF/HFC/TableProducer/correlatorD0Hadrons.cxx @@ -303,7 +303,7 @@ struct HfCorrelatorD0Hadrons { registry.fill(HIST("hMultiplicity"), nTracks); std::vector outputMlD0 = {-1., -1., -1.}; std::vector outputMlD0bar = {-1., -1., -1.}; - + for (const auto& candidate : candidates) { if (std::abs(hfHelper.yD0(candidate)) >= yCandMax || candidate.pt() <= ptCandMin || candidate.pt() >= ptTrackMax) { continue; @@ -339,8 +339,8 @@ struct HfCorrelatorD0Hadrons { outputMlD0bar[iclass] = candidate.mlProbD0bar()[classMl->at(iclass)]; } } - entryD0CandRecoInfo(hfHelper.invMassD0ToPiK(candidate), hfHelper.invMassD0barToKPi(candidate), candidate.pt(), outputMlD0[0], outputMlD0[2], outputMlD0bar[0], outputMlD0bar[2]); - + entryD0CandRecoInfo(hfHelper.invMassD0ToPiK(candidate), hfHelper.invMassD0barToKPi(candidate), candidate.pt(), outputMlD0[0], outputMlD0[2], outputMlD0bar[0], outputMlD0bar[2]); + // ========================== Fill general histos ================================ registry.fill(HIST("hPtCand"), candidate.pt()); registry.fill(HIST("hPtProng0"), candidate.ptProng0()); @@ -411,7 +411,7 @@ struct HfCorrelatorD0Hadrons { poolBin, correlationStatus); entryD0HadronRecoInfo(hfHelper.invMassD0ToPiK(candidate), hfHelper.invMassD0barToKPi(candidate), signalStatus); - entryD0HadronMlInfo(outputMlD0[0],outputMlD0[1],outputMlD0[2],outputMlD0bar[0],outputMlD0bar[1],outputMlD0bar[2]); + entryD0HadronMlInfo(outputMlD0[0], outputMlD0[1], outputMlD0[2], outputMlD0bar[0], outputMlD0bar[1], outputMlD0bar[2]); } // end inner loop (tracks) @@ -580,7 +580,7 @@ struct HfCorrelatorD0Hadrons { correlationStatus); entryD0HadronRecoInfo(hfHelper.invMassD0ToPiK(candidate), hfHelper.invMassD0barToKPi(candidate), signalStatus); } // end inner loop (Tracks) - } // end of outer loop (D0) + } // end of outer loop (D0) registry.fill(HIST("hZvtx"), collision.posZ()); registry.fill(HIST("hMultV0M"), collision.multFT0M()); } @@ -716,13 +716,13 @@ struct HfCorrelatorD0Hadrons { invMassDstar2 = std::sqrt((eKPi + ePion) * (eKPi + ePion) - pSum2); std::vector outputMlD0 = {-1., -1., -1.}; std::vector outputMlD0bar = {-1., -1., -1.}; - + if (t1.isSelD0() >= selectionFlagD0) { if ((std::abs(invMassDstar1 - hfHelper.invMassD0ToPiK(t1)) - softPiMass) < ptSoftPionMax) { isSoftPiD0 = true; } for (unsigned int iclass = 0; iclass < classMl->size(); iclass++) { - outputMlD0[iclass] = t1.mlProbD0()[classMl->at(iclass)]; + outputMlD0[iclass] = t1.mlProbD0()[classMl->at(iclass)]; } } if (t1.isSelD0bar() >= selectionFlagD0bar) { @@ -730,7 +730,7 @@ struct HfCorrelatorD0Hadrons { isSoftPiD0bar = true; } for (unsigned int iclass = 0; iclass < classMl->size(); iclass++) { - outputMlD0bar[iclass] = t1.mlProbD0bar()[classMl->at(iclass)]; + outputMlD0bar[iclass] = t1.mlProbD0bar()[classMl->at(iclass)]; } } @@ -752,7 +752,7 @@ struct HfCorrelatorD0Hadrons { bool correlationStatus = false; entryD0HadronPair(getDeltaPhi(t1.phi(), t2.phi()), t1.eta() - t2.eta(), t1.pt(), t2.pt(), poolBin, correlationStatus); entryD0HadronRecoInfo(hfHelper.invMassD0ToPiK(t1), hfHelper.invMassD0barToKPi(t1), signalStatus); - entryD0HadronMlInfo(outputMlD0[0],outputMlD0[1],outputMlD0[2],outputMlD0bar[0],outputMlD0bar[1],outputMlD0bar[2]); + entryD0HadronMlInfo(outputMlD0[0], outputMlD0[1], outputMlD0[2], outputMlD0bar[0], outputMlD0bar[1], outputMlD0bar[2]); } } } diff --git a/PWGHF/HFC/Tasks/taskCorrelationD0Hadrons.cxx b/PWGHF/HFC/Tasks/taskCorrelationD0Hadrons.cxx index bc99470cb54..5c1ca9852bd 100644 --- a/PWGHF/HFC/Tasks/taskCorrelationD0Hadrons.cxx +++ b/PWGHF/HFC/Tasks/taskCorrelationD0Hadrons.cxx @@ -264,7 +264,7 @@ struct HfTaskCorrelationD0Hadrons { registry.fill(HIST("hBdtScorePromptD0bar"), bdtScorePromptD0bar); registry.fill(HIST("hBdtScoreBkgD0bar"), bdtScoreBkgD0bar); } - + for (const auto& pairEntry : pairEntries) { // define variables for widely used quantities double deltaPhi = pairEntry.deltaPhi(); @@ -289,7 +289,7 @@ struct HfTaskCorrelationD0Hadrons { if (ptHadron > ptHadronMax) { ptHadron = ptHadronMax + 0.5; } - if (bdtScorePromptD0 < mlOutputPromptD0->at(ptBinD) || bdtScoreBkgD0 > mlOutputBkgD0->at(ptBinD) || + if (bdtScorePromptD0 < mlOutputPromptD0->at(ptBinD) || bdtScoreBkgD0 > mlOutputBkgD0->at(ptBinD) || bdtScorePromptD0bar < mlOutputPromptD0bar->at(ptBinD) || bdtScoreBkgD0bar > mlOutputBkgD0bar->at(ptBinD)) { continue; } From 184ed136260d5f1730fcaaec34747098feaa270a Mon Sep 17 00:00:00 2001 From: samrangy Date: Tue, 10 Dec 2024 22:34:38 +0100 Subject: [PATCH 09/11] Fix O2 Linter warnings --- PWGHF/HFC/DataModel/CorrelationTables.h | 1 + .../HFC/TableProducer/correlatorD0Hadrons.cxx | 30 +++++++++---------- PWGHF/HFC/Tasks/taskCorrelationD0Hadrons.cxx | 15 +++++----- 3 files changed, 24 insertions(+), 22 deletions(-) diff --git a/PWGHF/HFC/DataModel/CorrelationTables.h b/PWGHF/HFC/DataModel/CorrelationTables.h index 3b8e68cdbd8..71e983cf194 100644 --- a/PWGHF/HFC/DataModel/CorrelationTables.h +++ b/PWGHF/HFC/DataModel/CorrelationTables.h @@ -423,3 +423,4 @@ DECLARE_SOA_TABLE(HfEHadronPair, "AOD", "HFEHADRONPAIR", //! Hfe-Hadrons pairs I } // namespace o2::aod #endif // PWGHF_HFC_DATAMODEL_CORRELATIONTABLES_H_ + diff --git a/PWGHF/HFC/TableProducer/correlatorD0Hadrons.cxx b/PWGHF/HFC/TableProducer/correlatorD0Hadrons.cxx index c4043c8ed53..923f543a841 100644 --- a/PWGHF/HFC/TableProducer/correlatorD0Hadrons.cxx +++ b/PWGHF/HFC/TableProducer/correlatorD0Hadrons.cxx @@ -80,8 +80,6 @@ using SelectedParticlesMcGen = soa::Join d0Sel; Configurable selectionFlagD0{"selectionFlagD0", 1, "Selection Flag for D0"}; @@ -90,8 +88,9 @@ struct HfCorrelatorD0HadronsSelection { Configurable ptCandMin{"ptCandMin", -1., "min. cand. pT"}; HfHelper hfHelper; - + SliceCache cache; Preslice perCol = aod::hf_cand::collisionId; + Partition> selectedD0Candidates = aod::hf_sel_candidate_d0::isSelD0 >= selectionFlagD0 || aod::hf_sel_candidate_d0::isSelD0bar >= selectionFlagD0bar; Partition> selectedD0candidatesMc = aod::hf_sel_candidate_d0::isSelD0 >= selectionFlagD0 || aod::hf_sel_candidate_d0::isSelD0bar >= selectionFlagD0bar; @@ -162,8 +161,6 @@ struct HfCorrelatorD0HadronsSelection { /// D0-Hadron correlation pair builder - for real data and data-like analysis (i.e. reco-level w/o matching request via MC truth) struct HfCorrelatorD0Hadrons { - SliceCache cache; - Produces entryD0HadronPair; Produces entryD0HadronRecoInfo; Produces entryD0HadronMlInfo; @@ -178,9 +175,9 @@ struct HfCorrelatorD0Hadrons { Configurable ptCandMin{"ptCandMin", 1., "min. cand. pT"}; Configurable ptTrackMin{"ptTrackMin", 0.3, "min. track pT"}; Configurable ptTrackMax{"ptTrackMax", 99., "max. track pT"}; - Configurable> bins{"ptBinsForMassAndEfficiency", std::vector{o2::analysis::hf_cuts_d0_to_pi_k::vecBinsPt}, "pT bin limits for candidate mass plots and efficiency"}; + Configurable> bins{"bins", std::vector{o2::analysis::hf_cuts_d0_to_pi_k::vecBinsPt}, "pT bin limits for candidate mass plots and efficiency"}; Configurable> efficiencyDmeson{"efficiencyDmeson", std::vector{vecEfficiencyDmeson}, "Efficiency values for D0 meson"}; - Configurable applyEfficiency{"efficiencyFlagD", 1, "Flag for applying D-meson efficiency weights"}; + Configurable applyEfficiency{"applyEfficiency", 1, "Flag for applying D-meson efficiency weights"}; Configurable multMin{"multMin", 0., "minimum multiplicity accepted"}; Configurable multMax{"multMax", 10000., "maximum multiplicity accepted"}; Configurable ptSoftPionMax{"ptSoftPionMax", 3.f * 800.f * std::pow(10.f, -6.f), "max. pT cut for soft pion identification"}; @@ -188,12 +185,6 @@ struct HfCorrelatorD0Hadrons { Configurable correlateD0WithLeadingParticle{"correlateD0WithLeadingParticle", false, "Switch for correlation of D0 mesons with leading particle only"}; Configurable storeAutoCorrelationFlag{"storeAutoCorrelationFlag", false, "Store flag that indicates if the track is paired to its D-meson mother instead of skipping it"}; Configurable numberEventsMixed{"numberEventsMixed", 5, "Number of events mixed in ME process"}; - ConfigurableAxis zPoolBins{"zPoolBins", {VARIABLE_WIDTH, -10.0f, -2.5f, 2.5f, 10.0f}, "z vertex position pools"}; - ConfigurableAxis multPoolBins{"multPoolBins", {VARIABLE_WIDTH, 0.0f, 2000.0f, 6000.0f, 10000.0f}, "event multiplicity pools (FT0M)"}; - ConfigurableAxis multPoolBinsMcGen{"multPoolBinsMcGen", {VARIABLE_WIDTH, 0.0f, 20.0f, 50.0f, 500.0f}, "Mixing bins - MC multiplicity"}; // In MCGen multiplicity is defined by counting tracks - - HfHelper hfHelper; - BinningType corrBinning{{zPoolBins, multPoolBins}, true}; int leadingIndex = 0; double massD0{0.}; @@ -201,13 +192,21 @@ struct HfCorrelatorD0Hadrons { double massK{0.}; double softPiMass = 0.14543; // pion mass + Q-value of the D*->D0pi decay - Preslice perCol = aod::hf_cand::collisionId; - Filter collisionFilter = aod::hf_selection_dmeson_collision::dmesonSel == true; Filter trackFilter = requireGlobalTrackWoDCAInFilter() && (nabs(aod::track::eta) < etaTrackMax) && (aod::track::pt > ptTrackMin) && (aod::track::pt < ptTrackMax) && (nabs(aod::track::dcaXY) < dcaXYTrackMax) && (nabs(aod::track::dcaZ) < dcaZTrackMax); Filter d0Filter = (aod::hf_sel_candidate_d0::isSelD0 >= 1) || (aod::hf_sel_candidate_d0::isSelD0bar >= 1); Filter collisionFilterGen = aod::hf_selection_dmeson_collision::dmesonSel == true; Filter particlesFilter = nabs(aod::mcparticle::pdgCode) == static_cast(Pdg::kD0) || ((aod::mcparticle::flags & (uint8_t)o2::aod::mcparticle::enums::PhysicalPrimary) == (uint8_t)o2::aod::mcparticle::enums::PhysicalPrimary); + + HfHelper hfHelper; + SliceCache cache; + Preslice perCol = aod::hf_cand::collisionId; + + ConfigurableAxis zPoolBins{"zPoolBins", {VARIABLE_WIDTH, -10.0f, -2.5f, 2.5f, 10.0f}, "z vertex position pools"}; + ConfigurableAxis multPoolBins{"multPoolBins", {VARIABLE_WIDTH, 0.0f, 2000.0f, 6000.0f, 10000.0f}, "event multiplicity pools (FT0M)"}; + ConfigurableAxis multPoolBinsMcGen{"multPoolBinsMcGen", {VARIABLE_WIDTH, 0.0f, 20.0f, 50.0f, 500.0f}, "Mixing bins - MC multiplicity"}; // In MCGen multiplicity is defined by counting tracks + + BinningType corrBinning{{zPoolBins, multPoolBins}, true}; HistogramRegistry registry{ "registry", @@ -925,3 +924,4 @@ WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) adaptAnalysisTask(cfgc), adaptAnalysisTask(cfgc)}; } + diff --git a/PWGHF/HFC/Tasks/taskCorrelationD0Hadrons.cxx b/PWGHF/HFC/Tasks/taskCorrelationD0Hadrons.cxx index 5c1ca9852bd..8829c6684a5 100644 --- a/PWGHF/HFC/Tasks/taskCorrelationD0Hadrons.cxx +++ b/PWGHF/HFC/Tasks/taskCorrelationD0Hadrons.cxx @@ -86,14 +86,14 @@ const double ptHadronMax = 10.0; struct HfTaskCorrelationD0Hadrons { // pT ranges for correlation plots: the default values are those embedded in hf_cuts_d0_to_pi_k (i.e. the mass pT bins), but can be redefined via json files - Configurable> binsCorrelations{"ptBinsForCorrelations", std::vector{vecPtBinsCorrelations}, "pT bin limits for correlation plots"}; + Configurable> binsCorrelations{"binsCorrelations", std::vector{vecPtBinsCorrelations}, "pT bin limits for correlation plots"}; // pT bins for effiencies: same as above - Configurable> binsEfficiency{"ptBinsForEfficiency", std::vector{o2::analysis::hf_cuts_d0_to_pi_k::vecBinsPt}, "pT bin limits for efficiency"}; + Configurable> binsEfficiency{"binsEfficiency", std::vector{o2::analysis::hf_cuts_d0_to_pi_k::vecBinsPt}, "pT bin limits for efficiency"}; Configurable> classMl{"classMl", {0, 1, 2}, "Indexes of ML scores to be stored. Three indexes max."}; - Configurable> mlOutputPromptD0{"mlScorePromptD0", {0.5, 0.5, 0.5, 0.5}, "Machine learning scores for prompt"}; - Configurable> mlOutputBkgD0{"mlScoreBkgD0", {0.5, 0.5, 0.5, 0.5}, "Machine learning scores for bkg"}; - Configurable> mlOutputPromptD0bar{"mlScorePromptD0bar", {0.5, 0.5, 0.5, 0.5}, "Machine learning scores for prompt"}; - Configurable> mlOutputBkgD0bar{"mlScoreBkgD0bar", {0.5, 0.5, 0.5, 0.5}, "Machine learning scores for bkg"}; + Configurable> mlOutputPromptD0{"mlOutputPromptD0", {0.5, 0.5, 0.5, 0.5}, "Machine learning scores for prompt"}; + Configurable> mlOutputBkgD0{"mlOutputBkgD0", {0.5, 0.5, 0.5, 0.5}, "Machine learning scores for bkg"}; + Configurable> mlOutputPromptD0bar{"mlOutputPromptD0bar", {0.5, 0.5, 0.5, 0.5}, "Machine learning scores for prompt"}; + Configurable> mlOutputBkgD0bar{"mlOutputBkgD0bar", {0.5, 0.5, 0.5, 0.5}, "Machine learning scores for bkg"}; // signal and sideband region edges, to be defined via json file (initialised to empty) Configurable> signalRegionLeft{"signalRegionLeft", std::vector{vecsignalRegionLeft}, "Inner values of signal region vs pT"}; Configurable> signalRegionRight{"signalRegionRight", std::vector{vecsignalRegionRight}, "Outer values of signal region vs pT"}; @@ -104,7 +104,7 @@ struct HfTaskCorrelationD0Hadrons { Configurable> efficiencyDmeson{"efficiencyDmeson", std::vector{vecEfficiencyDmeson}, "Efficiency values for D meson specie under study"}; Configurable isTowardTransverseAway{"isTowardTransverseAway", false, "Divide into three regions: toward, transverse, and away"}; Configurable leadingParticlePtMin{"leadingParticlePtMin", 0., "Min for leading particle pt"}; - Configurable applyEfficiency{"efficiencyFlagD", 1, "Flag for applying efficiency weights"}; + Configurable applyEfficiency{"applyEfficiency", 1, "Flag for applying efficiency weights"}; HistogramRegistry registry{ "registry", @@ -668,3 +668,4 @@ WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{adaptAnalysisTask(cfgc)}; } + From cdc0dae3b871decb9639f733d7aba59ed07c9b3a Mon Sep 17 00:00:00 2001 From: samrangy Date: Tue, 10 Dec 2024 22:42:20 +0100 Subject: [PATCH 10/11] Fix O2 Linter warnings --- PWGHF/HFC/TableProducer/correlatorD0Hadrons.cxx | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/PWGHF/HFC/TableProducer/correlatorD0Hadrons.cxx b/PWGHF/HFC/TableProducer/correlatorD0Hadrons.cxx index 923f543a841..dea6d32ea02 100644 --- a/PWGHF/HFC/TableProducer/correlatorD0Hadrons.cxx +++ b/PWGHF/HFC/TableProducer/correlatorD0Hadrons.cxx @@ -192,14 +192,15 @@ struct HfCorrelatorD0Hadrons { double massK{0.}; double softPiMass = 0.14543; // pion mass + Q-value of the D*->D0pi decay + HfHelper hfHelper; + SliceCache cache; + Filter collisionFilter = aod::hf_selection_dmeson_collision::dmesonSel == true; Filter trackFilter = requireGlobalTrackWoDCAInFilter() && (nabs(aod::track::eta) < etaTrackMax) && (aod::track::pt > ptTrackMin) && (aod::track::pt < ptTrackMax) && (nabs(aod::track::dcaXY) < dcaXYTrackMax) && (nabs(aod::track::dcaZ) < dcaZTrackMax); Filter d0Filter = (aod::hf_sel_candidate_d0::isSelD0 >= 1) || (aod::hf_sel_candidate_d0::isSelD0bar >= 1); Filter collisionFilterGen = aod::hf_selection_dmeson_collision::dmesonSel == true; Filter particlesFilter = nabs(aod::mcparticle::pdgCode) == static_cast(Pdg::kD0) || ((aod::mcparticle::flags & (uint8_t)o2::aod::mcparticle::enums::PhysicalPrimary) == (uint8_t)o2::aod::mcparticle::enums::PhysicalPrimary); - HfHelper hfHelper; - SliceCache cache; Preslice perCol = aod::hf_cand::collisionId; ConfigurableAxis zPoolBins{"zPoolBins", {VARIABLE_WIDTH, -10.0f, -2.5f, 2.5f, 10.0f}, "z vertex position pools"}; From 0afc62d55bc6751f0921c28a1190f934364e7346 Mon Sep 17 00:00:00 2001 From: ALICE Action Bot Date: Tue, 10 Dec 2024 21:57:04 +0000 Subject: [PATCH 11/11] Please consider the following formatting changes --- PWGHF/HFC/DataModel/CorrelationTables.h | 1 - PWGHF/HFC/TableProducer/correlatorD0Hadrons.cxx | 9 ++++----- PWGHF/HFC/Tasks/taskCorrelationD0Hadrons.cxx | 1 - 3 files changed, 4 insertions(+), 7 deletions(-) diff --git a/PWGHF/HFC/DataModel/CorrelationTables.h b/PWGHF/HFC/DataModel/CorrelationTables.h index 71e983cf194..3b8e68cdbd8 100644 --- a/PWGHF/HFC/DataModel/CorrelationTables.h +++ b/PWGHF/HFC/DataModel/CorrelationTables.h @@ -423,4 +423,3 @@ DECLARE_SOA_TABLE(HfEHadronPair, "AOD", "HFEHADRONPAIR", //! Hfe-Hadrons pairs I } // namespace o2::aod #endif // PWGHF_HFC_DATAMODEL_CORRELATIONTABLES_H_ - diff --git a/PWGHF/HFC/TableProducer/correlatorD0Hadrons.cxx b/PWGHF/HFC/TableProducer/correlatorD0Hadrons.cxx index dea6d32ea02..3500b09dc92 100644 --- a/PWGHF/HFC/TableProducer/correlatorD0Hadrons.cxx +++ b/PWGHF/HFC/TableProducer/correlatorD0Hadrons.cxx @@ -90,7 +90,7 @@ struct HfCorrelatorD0HadronsSelection { HfHelper hfHelper; SliceCache cache; Preslice perCol = aod::hf_cand::collisionId; - + Partition> selectedD0Candidates = aod::hf_sel_candidate_d0::isSelD0 >= selectionFlagD0 || aod::hf_sel_candidate_d0::isSelD0bar >= selectionFlagD0bar; Partition> selectedD0candidatesMc = aod::hf_sel_candidate_d0::isSelD0 >= selectionFlagD0 || aod::hf_sel_candidate_d0::isSelD0bar >= selectionFlagD0bar; @@ -194,15 +194,15 @@ struct HfCorrelatorD0Hadrons { HfHelper hfHelper; SliceCache cache; - + Filter collisionFilter = aod::hf_selection_dmeson_collision::dmesonSel == true; Filter trackFilter = requireGlobalTrackWoDCAInFilter() && (nabs(aod::track::eta) < etaTrackMax) && (aod::track::pt > ptTrackMin) && (aod::track::pt < ptTrackMax) && (nabs(aod::track::dcaXY) < dcaXYTrackMax) && (nabs(aod::track::dcaZ) < dcaZTrackMax); Filter d0Filter = (aod::hf_sel_candidate_d0::isSelD0 >= 1) || (aod::hf_sel_candidate_d0::isSelD0bar >= 1); Filter collisionFilterGen = aod::hf_selection_dmeson_collision::dmesonSel == true; Filter particlesFilter = nabs(aod::mcparticle::pdgCode) == static_cast(Pdg::kD0) || ((aod::mcparticle::flags & (uint8_t)o2::aod::mcparticle::enums::PhysicalPrimary) == (uint8_t)o2::aod::mcparticle::enums::PhysicalPrimary); - + Preslice perCol = aod::hf_cand::collisionId; - + ConfigurableAxis zPoolBins{"zPoolBins", {VARIABLE_WIDTH, -10.0f, -2.5f, 2.5f, 10.0f}, "z vertex position pools"}; ConfigurableAxis multPoolBins{"multPoolBins", {VARIABLE_WIDTH, 0.0f, 2000.0f, 6000.0f, 10000.0f}, "event multiplicity pools (FT0M)"}; ConfigurableAxis multPoolBinsMcGen{"multPoolBinsMcGen", {VARIABLE_WIDTH, 0.0f, 20.0f, 50.0f, 500.0f}, "Mixing bins - MC multiplicity"}; // In MCGen multiplicity is defined by counting tracks @@ -925,4 +925,3 @@ WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) adaptAnalysisTask(cfgc), adaptAnalysisTask(cfgc)}; } - diff --git a/PWGHF/HFC/Tasks/taskCorrelationD0Hadrons.cxx b/PWGHF/HFC/Tasks/taskCorrelationD0Hadrons.cxx index 8829c6684a5..2654984c45a 100644 --- a/PWGHF/HFC/Tasks/taskCorrelationD0Hadrons.cxx +++ b/PWGHF/HFC/Tasks/taskCorrelationD0Hadrons.cxx @@ -668,4 +668,3 @@ WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{adaptAnalysisTask(cfgc)}; } -