diff --git a/Analysis/DataModel/include/AnalysisDataModel/HFSecondaryVertex.h b/Analysis/DataModel/include/AnalysisDataModel/HFSecondaryVertex.h index df8d4b8c48d71..d4a1b16604a81 100644 --- a/Analysis/DataModel/include/AnalysisDataModel/HFSecondaryVertex.h +++ b/Analysis/DataModel/include/AnalysisDataModel/HFSecondaryVertex.h @@ -644,6 +644,26 @@ DECLARE_SOA_TABLE(HfCandProng3MCGen, "AOD", "HFCANDP3MCGEN", //! hf_cand_prong3::OriginMCGen, hf_cand_prong3::FlagMCDecayChanGen); +// definition of columns and tables for D-Dbar correlation pairs +namespace hf_correlation_ddbar +{ +DECLARE_SOA_COLUMN(DeltaPhi, deltaPhi, float); +DECLARE_SOA_COLUMN(DeltaEta, deltaEta, float); +DECLARE_SOA_COLUMN(PtD, ptD, float); +DECLARE_SOA_COLUMN(PtDbar, ptDbar, float); +DECLARE_SOA_COLUMN(MD, mD, float); +DECLARE_SOA_COLUMN(MDbar, mDbar, float); +DECLARE_SOA_COLUMN(SignalStatus, signalStatus, int); +} // namespace hf_correlation_ddbar +DECLARE_SOA_TABLE(DDbarPair, "AOD", "DDBARPAIR", + aod::hf_correlation_ddbar::DeltaPhi, + aod::hf_correlation_ddbar::DeltaEta, + aod::hf_correlation_ddbar::PtD, + aod::hf_correlation_ddbar::PtDbar); +DECLARE_SOA_TABLE(DDbarRecoInfo, "AOD", "DDBARRECOINFO", + aod::hf_correlation_ddbar::MD, + aod::hf_correlation_ddbar::MDbar, + aod::hf_correlation_ddbar::SignalStatus); } // namespace o2::aod #endif // O2_ANALYSIS_HFSECONDARYVERTEX_H_ diff --git a/Analysis/Tasks/PWGHF/CMakeLists.txt b/Analysis/Tasks/PWGHF/CMakeLists.txt index 37517c72640bf..5f782f71b211e 100644 --- a/Analysis/Tasks/PWGHF/CMakeLists.txt +++ b/Analysis/Tasks/PWGHF/CMakeLists.txt @@ -80,9 +80,9 @@ o2_add_dpl_workflow(hf-xic-topkpi-candidate-selector o2_add_dpl_workflow(hf-lc-tok0sp-candidate-selector SOURCES HFLcK0sPCandidateSelector.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisDataModel O2::AnalysisCore O2::DetectorsVertexing O2::AnalysisTasksUtils + PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisDataModel O2::AnalysisCore O2::DetectorsVertexing O2::AnalysisTasksUtils COMPONENT_NAME Analysis) - + o2_add_dpl_workflow(hf-task-d0 SOURCES taskD0.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisDataModel O2::AnalysisCore O2::DetectorsVertexing @@ -125,6 +125,20 @@ o2_add_dpl_workflow(hf-mc-validation o2_add_dpl_workflow(hf-task-lc-tok0sp SOURCES taskLcK0sP.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisDataModel O2::AnalysisCore O2::DetectorsVertexing O2::AnalysisTasksUtils + PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisDataModel O2::AnalysisCore O2::DetectorsVertexing O2::AnalysisTasksUtils COMPONENT_NAME Analysis) +o2_add_dpl_workflow(hf-correlator-d0d0bar + SOURCES HFCorrelatorD0D0bar.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisDataModel O2::AnalysisCore O2::DetectorsVertexing + COMPONENT_NAME Analysis) + +o2_add_dpl_workflow(hf-correlator-dplusdminus + SOURCES HFCorrelatorDplusDminus.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisDataModel O2::AnalysisCore O2::DetectorsVertexing + COMPONENT_NAME Analysis) + +o2_add_dpl_workflow(hf-task-correlation-ddbar + SOURCES taskCorrelationDDbar.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisDataModel O2::AnalysisCore O2::DetectorsVertexing + COMPONENT_NAME Analysis) diff --git a/Analysis/Tasks/PWGHF/HFCorrelatorD0D0bar.cxx b/Analysis/Tasks/PWGHF/HFCorrelatorD0D0bar.cxx new file mode 100644 index 0000000000000..c2f2bf08ff7aa --- /dev/null +++ b/Analysis/Tasks/PWGHF/HFCorrelatorD0D0bar.cxx @@ -0,0 +1,871 @@ +// Copyright CERN and copyright holders of ALICE O2. This software is +// distributed under the terms of the GNU General Public License v3 (GPL +// Version 3), copied verbatim in the file "COPYING". +// +// See http://alice-o2.web.cern.ch/license for full licensing information. +// +// 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 HFCorrelatorD0D0bar.cxx +/// \brief D0-D0bar correlator task - data-like, MC-reco and MC-kine analyses. For ULS and LS pairs +/// +/// \author Fabio Colamaria , INFN Bari + +#include "Framework/AnalysisTask.h" +#include "Framework/HistogramRegistry.h" +#include "AnalysisDataModel/HFSecondaryVertex.h" +#include "AnalysisDataModel/HFCandidateSelectionTables.h" + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; +using namespace o2::aod::hf_cand_prong2; +using namespace o2::aod::hf_correlation_ddbar; +using namespace o2::analysis::hf_cuts_d0_topik; +using namespace o2::constants::math; + +void customize(std::vector& workflowOptions) +{ + ConfigParamSpec optionDoLikeSign{"doLikeSign", VariantType::Bool, false, {"Run Like-Sign analysis."}}; + ConfigParamSpec optionDoMCccbar{"doMCccbar", VariantType::Bool, false, {"Run MC-Gen dedicated tasks."}}; + ConfigParamSpec optionDoMCGen{"doMCGen", VariantType::Bool, false, {"Run MC-Gen dedicated tasks."}}; + ConfigParamSpec optionDoMCRec{"doMCRec", VariantType::Bool, false, {"Run MC-Rec dedicated tasks."}}; + workflowOptions.push_back(optionDoLikeSign); + workflowOptions.push_back(optionDoMCccbar); + workflowOptions.push_back(optionDoMCGen); + workflowOptions.push_back(optionDoMCRec); +} + +#include "Framework/runDataProcessing.h" + +/// +/// Returns deltaPhi value in range [-pi/2., 3.*pi/2], typically used for correlation studies +/// +double getDeltaPhi(double phiD, double phiDbar) +{ + return RecoDecay::constrainAngle(phiDbar - phiD, -o2::constants::math::PI / 2.); +} + +/// definition of variables for D0D0bar pairs vs eta acceptance studies (hDDbarVsEtaCut, in data-like, MC-reco and MC-kine tasks) +const double maxEtaCut = 5.; +const double ptThresholdForMaxEtaCut = 10.; +const double incrementEtaCut = 0.1; +const double incrementPtThreshold = 0.5; +const double epsilon = 1E-5; + +/// D0-D0bar correlation pair builder - for real data and data-like analysis (i.e. reco-level w/o matching request via MC truth) +struct HfCorrelatorD0D0bar { + Produces entryD0D0barPair; + Produces entryD0D0barRecoInfo; + + 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, {{180, 0., 36.}}}}, + {"hPtProng0", "D0,D0bar candidates;prong 0 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, + {"hPtProng1", "D0,D0bar candidates;prong 1 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, + {"hSelectionStatus", "D0,D0bar candidates;selection status;entries", {HistType::kTH1F, {{4, -0.5, 3.5}}}}, + {"hEta", "D0,D0bar candidates;candidate #it{#eta};entries", {HistType::kTH1F, {{100, -5., 5.}}}}, + {"hPhi", "D0,D0bar candidates;candidate #it{#varphi};entries", {HistType::kTH1F, {{32, 0., 2. * o2::constants::math::PI}}}}, + {"hY", "D0,D0bar candidates;candidate #it{#y};entries", {HistType::kTH1F, {{100, -5., 5.}}}}, + {"hDDbarVsEtaCut", "D0,D0bar pairs vs #eta cut;#eta_{max};entries", {HistType::kTH2F, {{(int)(maxEtaCut / incrementEtaCut), 0., maxEtaCut}, {(int)(ptThresholdForMaxEtaCut / incrementPtThreshold), 0., ptThresholdForMaxEtaCut}}}}}}; + + Configurable selectionFlagD0{"selectionFlagD0", 1, "Selection Flag for D0"}; + Configurable selectionFlagD0bar{"selectionFlagD0bar", 1, "Selection Flag for D0bar"}; + Configurable cutYCandMax{"cutYCandMax", -1., "max. cand. rapidity"}; + Configurable cutPtCandMin{"cutPtCandMin", -1., "min. cand. pT"}; + Configurable> bins{"ptBinsForMass", std::vector{o2::analysis::hf_cuts_d0_topik::pTBins_v}, "pT bin limits for candidate mass plots"}; + + void init(o2::framework::InitContext&) + { + registry.add("hMass", "D0,D0bar candidates;inv. mass (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{120, 1.5848, 2.1848}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hMassD0", "D0,D0bar candidates;inv. mass D0 only (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{120, 1.5848, 2.1848}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hMassD0bar", "D0,D0bar candidates;inv. mass D0bar only (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{120, 1.5848, 2.1848}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + } + + Filter filterSelectCandidates = (aod::hf_selcandidate_d0::isSelD0 >= selectionFlagD0 || aod::hf_selcandidate_d0::isSelD0bar >= selectionFlagD0bar); + + void process(aod::Collision const& collision, soa::Filtered> const& candidates) + { + for (auto& candidate1 : candidates) { + if (cutYCandMax >= 0. && std::abs(YD0(candidate1)) > cutYCandMax) { + continue; + } + if (cutPtCandMin >= 0. && candidate1.pt() < cutPtCandMin) { + continue; + } + //check decay channel flag for candidate1 + if (!(candidate1.hfflag() & 1 << DecayType::D0ToPiK)) { + continue; + } + //fill invariant mass plots and generic info from all D0/D0bar candidates + if (candidate1.isSelD0() >= selectionFlagD0) { + registry.fill(HIST("hMass"), InvMassD0(candidate1), candidate1.pt()); + registry.fill(HIST("hMassD0"), InvMassD0(candidate1), candidate1.pt()); + } + if (candidate1.isSelD0bar() >= selectionFlagD0bar) { + registry.fill(HIST("hMass"), InvMassD0bar(candidate1), candidate1.pt()); + registry.fill(HIST("hMassD0bar"), InvMassD0bar(candidate1), candidate1.pt()); + } + 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"), YD0(candidate1)); + registry.fill(HIST("hSelectionStatus"), candidate1.isSelD0bar() + (candidate1.isSelD0() * 2)); + + //D-Dbar correlation dedicated section + //if the candidate is a D0, search for D0bar and evaluate correlations + if (candidate1.isSelD0() < selectionFlagD0) { + continue; + } + for (auto& candidate2 : candidates) { + if (!(candidate2.hfflag() & 1 << DecayType::D0ToPiK)) { //check decay channel flag for candidate2 + continue; + } + if (candidate2.isSelD0bar() < selectionFlagD0bar) { //keep only D0bar candidates passing the selection + continue; + } + //kinematic selection on D0bar candidates + if (cutYCandMax >= 0. && std::abs(YD0(candidate2)) > cutYCandMax) { + continue; + } + if (cutPtCandMin >= 0. && candidate2.pt() < cutPtCandMin) { + continue; + } + //excluding trigger self-correlations (possible in case of both mass hypotheses accepted) + if (candidate1.mRowIndex == candidate2.mRowIndex) { + continue; + } + entryD0D0barPair(getDeltaPhi(candidate2.phi(), candidate1.phi()), + candidate2.eta() - candidate1.eta(), + candidate1.pt(), + candidate2.pt()); + entryD0D0barRecoInfo(InvMassD0(candidate1), + InvMassD0bar(candidate2), + 0); + double etaCut = 0.; + double ptCut = 0.; + do { //fill pairs vs etaCut plot + ptCut = 0.; + etaCut += incrementEtaCut; + do { //fill pairs vs etaCut plot + if (std::abs(candidate1.eta()) < etaCut && std::abs(candidate2.eta()) < etaCut && candidate1.pt() > ptCut && candidate2.pt() > ptCut) { + registry.fill(HIST("hDDbarVsEtaCut"), etaCut - epsilon, ptCut + epsilon); + } + ptCut += incrementPtThreshold; + } while (ptCut < ptThresholdForMaxEtaCut - epsilon); + } while (etaCut < maxEtaCut - epsilon); + //note: candidates selected as both D0 and D0bar are used, and considered in both situation (but not auto-correlated): reflections could play a relevant role. + //another, more restrictive, option, could be to consider only candidates selected with a single option (D0 xor D0bar) + + } // end inner loop (Dbars) + + } //end outer loop + } +}; + +/// D0-D0bar correlation pair builder - for MC reco-level analysis (candidates matched to true signal only, but also the various bkg sources are studied) +struct HfCorrelatorD0D0barMcRec { + + Produces entryD0D0barPair; + Produces entryD0D0barRecoInfo; + + HistogramRegistry registry{ + "registry", + //NOTE: use hMassD0 for trigger normalisation (S*0.955), and hMass2DCorrelationPairs (in final task) for 2D-sideband-subtraction purposes + {{"hPtCandMCRec", "D0,D0bar candidates - MC reco;candidate #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, + {"hPtProng0MCRec", "D0,D0bar candidates - MC reco;prong 0 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, + {"hPtProng1MCRec", "D0,D0bar candidates - MC reco;prong 1 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, + {"hSelectionStatusMCRec", "D0,D0bar candidates - MC reco;selection status;entries", {HistType::kTH1F, {{4, -0.5, 3.5}}}}, + {"hEtaMCRec", "D0,D0bar candidates - MC reco;candidate #it{#eta};entries", {HistType::kTH1F, {{100, -5., 5.}}}}, + {"hPhiMCRec", "D0,D0bar candidates - MC reco;candidate #it{#varphi};entries", {HistType::kTH1F, {{32, 0., 2. * o2::constants::math::PI}}}}, + {"hYMCRec", "D0,D0bar candidates - MC reco;candidate #it{#y};entries", {HistType::kTH1F, {{100, -5., 5.}}}}, + {"hDDbarVsEtaCut", "D0,D0bar pairs vs #eta cut;#eta_{max};entries", {HistType::kTH2F, {{(int)(maxEtaCut / incrementEtaCut), 0., maxEtaCut}, {(int)(ptThresholdForMaxEtaCut / incrementPtThreshold), 0., ptThresholdForMaxEtaCut}}}}}}; + + Configurable selectionFlagD0{"selectionFlagD0", 1, "Selection Flag for D0"}; + Configurable selectionFlagD0bar{"selectionFlagD0bar", 1, "Selection Flag for D0bar"}; + Configurable cutYCandMax{"cutYCandMax", -1., "max. cand. rapidity"}; + Configurable cutPtCandMin{"cutPtCandMin", -1., "min. cand. pT"}; + Configurable> bins{"ptBinsForMass", std::vector{o2::analysis::hf_cuts_d0_topik::pTBins_v}, "pT bin limits for candidate mass plots"}; + + void init(o2::framework::InitContext&) + { + registry.add("hMassD0MCRec", "D0,D0bar candidates - MC reco;inv. mass (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{120, 1.5848, 2.1848}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hMassD0barMCRec", "D0,D0bar candidates - MC reco;inv. mass D0 only (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{120, 1.5848, 2.1848}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + } + + Filter filterSelectCandidates = (aod::hf_selcandidate_d0::isSelD0 >= selectionFlagD0 || aod::hf_selcandidate_d0::isSelD0bar >= selectionFlagD0bar); + + void process(aod::Collision const& collision, soa::Filtered> const& candidates) + { + //MC reco level + bool flagD0Signal = false; + bool flagD0barSignal = false; + for (auto& candidate1 : candidates) { + //check decay channel flag for candidate1 + if (!(candidate1.hfflag() & 1 << DecayType::D0ToPiK)) { + continue; + } + if (cutYCandMax >= 0. && std::abs(YD0(candidate1)) > cutYCandMax) { + continue; + } + if (cutPtCandMin >= 0. && candidate1.pt() < cutPtCandMin) { + continue; + } + if (std::abs(candidate1.flagMCMatchRec()) == 1 << DecayType::D0ToPiK) { + //fill invariant mass plots and generic info from all D0/D0bar candidates + if (candidate1.isSelD0() >= selectionFlagD0 && candidate1.flagMCMatchRec() == 1 << DecayType::D0ToPiK) { //only reco and matched as D0 + registry.fill(HIST("hMassD0MCRec"), InvMassD0(candidate1), candidate1.pt()); + } + if (candidate1.isSelD0bar() >= selectionFlagD0bar && candidate1.flagMCMatchRec() == -(1 << DecayType::D0ToPiK)) { //only reco and matched as D0bar + registry.fill(HIST("hMassD0barMCRec"), InvMassD0bar(candidate1), candidate1.pt()); + } + registry.fill(HIST("hPtCandMCRec"), candidate1.pt()); + registry.fill(HIST("hPtProng0MCRec"), candidate1.ptProng0()); + registry.fill(HIST("hPtProng1MCRec"), candidate1.ptProng1()); + registry.fill(HIST("hEtaMCRec"), candidate1.eta()); + registry.fill(HIST("hPhiMCRec"), candidate1.phi()); + registry.fill(HIST("hYMCRec"), YD0(candidate1)); + registry.fill(HIST("hSelectionStatusMCRec"), candidate1.isSelD0bar() + (candidate1.isSelD0() * 2)); + } + + //D-Dbar correlation dedicated section + //if the candidate is selected ad D0, search for D0bar and evaluate correlations + if (candidate1.isSelD0() < selectionFlagD0) { //discard candidates not selected as D0 in outer loop + continue; + } + flagD0Signal = candidate1.flagMCMatchRec() == 1 << DecayType::D0ToPiK; //flagD0Signal 'true' if candidate1 matched to D0 (particle) + for (auto& candidate2 : candidates) { + if (!(candidate2.hfflag() & 1 << DecayType::D0ToPiK)) { //check decay channel flag for candidate2 + continue; + } + if (candidate2.isSelD0bar() < selectionFlagD0bar) { //discard candidates not selected as D0bar in inner loop + continue; + } + flagD0barSignal = candidate2.flagMCMatchRec() == -(1 << DecayType::D0ToPiK); //flagD0barSignal 'true' if candidate2 matched to D0 (particle) + if (cutYCandMax >= 0. && std::abs(YD0(candidate2)) > cutYCandMax) { + continue; + } + if (cutPtCandMin >= 0. && candidate2.pt() < cutPtCandMin) { + continue; + } + //Excluding trigger self-correlations (possible in case of both mass hypotheses accepted) + if (candidate1.mRowIndex == candidate2.mRowIndex) { + continue; + } + //choice of options (D0/D0bar signal/bkg) + int pairSignalStatus = 0; //0 = bkg/bkg, 1 = bkg/sig, 2 = sig/bkg, 3 = sig/sig + if (flagD0Signal) { + pairSignalStatus += 2; + } + if (flagD0barSignal) { + pairSignalStatus += 1; + } + entryD0D0barPair(getDeltaPhi(candidate2.phi(), candidate1.phi()), + candidate2.eta() - candidate1.eta(), + candidate1.pt(), + candidate2.pt()); + entryD0D0barRecoInfo(InvMassD0(candidate1), + InvMassD0bar(candidate2), + pairSignalStatus); + double etaCut = 0.; + double ptCut = 0.; + do { //fill pairs vs etaCut plot + ptCut = 0.; + etaCut += incrementEtaCut; + do { //fill pairs vs etaCut plot + if (std::abs(candidate1.eta()) < etaCut && std::abs(candidate2.eta()) < etaCut && candidate1.pt() > ptCut && candidate2.pt() > ptCut) { + registry.fill(HIST("hDDbarVsEtaCut"), etaCut - epsilon, ptCut + epsilon); + } + ptCut += incrementPtThreshold; + } while (ptCut < ptThresholdForMaxEtaCut - epsilon); + } while (etaCut < maxEtaCut - epsilon); + } // end inner loop (Dbars) + + } //end outer loop + } +}; + +/// D0-D0bar correlation pair builder - for MC gen-level analysis (no filter/selection, only true signal) +struct HfCorrelatorD0D0barMcGen { + + Produces entryD0D0barPair; + + HistogramRegistry registry{ + "registry", + {{"hMCEvtCount", "Event counter - MC gen;;entries", {HistType::kTH1F, {{1, -0.5, 0.5}}}}, + {"hPtCandMCGen", "D0,D0bar particles - MC gen;particle #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, + {"hEtaMCGen", "D0,D0bar particles - MC gen;particle #it{#eta};entries", {HistType::kTH1F, {{100, -5., 5.}}}}, + {"hPhiMCGen", "D0,D0bar particles - MC gen;particle #it{#varphi};entries", {HistType::kTH1F, {{32, 0., 2. * o2::constants::math::PI}}}}, + {"hYMCGen", "D0,D0bar candidates - MC gen;candidate #it{#y};entries", {HistType::kTH1F, {{100, -5., 5.}}}}, + {"hcountD0D0barPerEvent", "D0,D0bar particles - MC gen;Number per event;entries", {HistType::kTH1F, {{20, 0., 20.}}}}, + {"hDDbarVsEtaCut", "D0,D0bar pairs vs #eta cut of D mesons;#eta_{max};entries", {HistType::kTH2F, {{(int)(maxEtaCut / incrementEtaCut), 0., maxEtaCut}, {(int)(ptThresholdForMaxEtaCut / incrementPtThreshold), 0., ptThresholdForMaxEtaCut}}}}, + {"hDDbarVsDaughterEtaCut", "D0,D0bar pairs vs #eta cut on D daughters;#eta_{max};entries", {HistType::kTH2F, {{(int)(maxEtaCut / incrementEtaCut), 0., maxEtaCut}, {(int)(ptThresholdForMaxEtaCut / incrementPtThreshold), 0., ptThresholdForMaxEtaCut}}}}}}; + + Configurable cutYCandMax{"cutYCandMax", -1., "max. cand. rapidity"}; + Configurable cutPtCandMin{"cutPtCandMin", -1., "min. cand. pT"}; + Configurable> bins{"ptBinsForMass", std::vector{o2::analysis::hf_cuts_d0_topik::pTBins_v}, "pT bin limits for trigger counters"}; + + void init(o2::framework::InitContext&) + { + registry.add("hcountD0triggersMCGen", "D0 trigger particles - MC gen;;N of trigger D0", {HistType::kTH2F, {{1, -0.5, 0.5}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + } + + void process(aod::McCollision const& mccollision, soa::Join const& particlesMC) + { + int counterD0D0bar = 0; + registry.fill(HIST("hMCEvtCount"), 0); + //MC gen level + for (auto& particle1 : particlesMC) { + //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::Code::kD0) { + continue; + } + double yD = RecoDecay::Y(array{particle1.px(), particle1.py(), particle1.pz()}, RecoDecay::getMassPDG(particle1.pdgCode())); + if (cutYCandMax >= 0. && std::abs(yD) > cutYCandMax) { + continue; + } + if (cutPtCandMin >= 0. && particle1.pt() < cutPtCandMin) { + continue; + } + registry.fill(HIST("hPtCandMCGen"), particle1.pt()); + registry.fill(HIST("hEtaMCGen"), particle1.eta()); + registry.fill(HIST("hPhiMCGen"), particle1.phi()); + registry.fill(HIST("hYMCGen"), yD); + counterD0D0bar++; + + //D-Dbar correlation dedicated section + //if it's a D0 particle, search for D0bar and evaluate correlations + if (particle1.pdgCode() != pdg::Code::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("hcountD0triggersMCGen"), 0, particle1.pt()); //to count trigger D0 (for normalisation) + for (auto& particle2 : particlesMC) { + if (particle2.pdgCode() != pdg::Code::kD0bar) { //check that inner particle is D0bar + continue; + } + if (cutYCandMax >= 0. && std::abs(RecoDecay::Y(array{particle2.px(), particle2.py(), particle2.pz()}, RecoDecay::getMassPDG(particle2.pdgCode()))) > cutYCandMax) { + continue; + } + if (cutPtCandMin >= 0. && particle2.pt() < cutPtCandMin) { + continue; + } + entryD0D0barPair(getDeltaPhi(particle2.phi(), particle1.phi()), + particle2.eta() - particle1.eta(), + particle1.pt(), + particle2.pt()); + double etaCut = 0.; + double ptCut = 0.; + + //fill pairs vs etaCut plot + bool rightDecayChannels = false; + if ((std::abs(particle1.flagMCMatchGen()) == 1 << DecayType::D0ToPiK) && (std::abs(particle2.flagMCMatchGen()) == 1 << DecayType::D0ToPiK)) { + rightDecayChannels = true; + } + do { + ptCut = 0.; + etaCut += incrementEtaCut; + do { //fill pairs vs etaCut plot + if (std::abs(particle1.eta()) < etaCut && std::abs(particle2.eta()) < etaCut && particle1.pt() > ptCut && particle2.pt() > ptCut) { //fill with D and Dbar acceptance checks + registry.fill(HIST("hDDbarVsEtaCut"), etaCut - epsilon, ptCut + epsilon); + } + if (rightDecayChannels) { //fill with D and Dbar daughter particls acceptance checks + double etaCandidate1Daughter1 = particlesMC.iteratorAt(particle1.daughter0()).eta(); + double etaCandidate1Daughter2 = particlesMC.iteratorAt(particle1.daughter1()).eta(); + double etaCandidate2Daughter1 = particlesMC.iteratorAt(particle2.daughter0()).eta(); + double etaCandidate2Daughter2 = particlesMC.iteratorAt(particle2.daughter1()).eta(); + if (std::abs(etaCandidate1Daughter1) < etaCut && std::abs(etaCandidate1Daughter2) < etaCut && + std::abs(etaCandidate2Daughter1) < etaCut && std::abs(etaCandidate2Daughter2) < etaCut && + particle1.pt() > ptCut && particle2.pt() > ptCut) { + registry.fill(HIST("hDDbarVsDaughterEtaCut"), etaCut - epsilon, ptCut + epsilon); + } + } + ptCut += incrementPtThreshold; + } while (ptCut < ptThresholdForMaxEtaCut - epsilon); + } while (etaCut < maxEtaCut - epsilon); + } //end inner loop + } //end outer loop + registry.fill(HIST("hcountD0D0barPerEvent"), counterD0D0bar); + } +}; + +/// D0-D0bar correlation pair builder - LIKE SIGN - for real data and data-like analysis (i.e. reco-level w/o matching request via MC truth) +/// NOTE: At the moment, both dPhi-symmetrical correlation pairs (part1-part2 and part2-part1) are filled, +/// since we bin in pT and selecting as trigger the largest pT particle would bias the distributions w.r.t. the ULS case. +struct HfCorrelatorD0D0barLs { + + Produces entryD0D0barPair; + Produces entryD0D0barRecoInfo; + + 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, {{180, 0., 36.}}}}, + {"hPtProng0", "D0,D0bar candidates;prong 0 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, + {"hPtProng1", "D0,D0bar candidates;prong 1 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, + {"hSelectionStatus", "D0,D0bar candidates;selection status;entries", {HistType::kTH1F, {{4, -0.5, 3.5}}}}, + {"hEta", "D0,D0bar candidates;candidate #it{#eta};entries", {HistType::kTH1F, {{100, -5., 5.}}}}, + {"hPhi", "D0,D0bar candidates;candidate #it{#varphi};entries", {HistType::kTH1F, {{32, 0., 2. * o2::constants::math::PI}}}}, + {"hY", "D0,D0bar candidates;candidate #it{#y};entries", {HistType::kTH1F, {{100, -5., 5.}}}}}}; + + Configurable selectionFlagD0{"selectionFlagD0", 1, "Selection Flag for D0"}; + Configurable selectionFlagD0bar{"selectionFlagD0bar", 1, "Selection Flag for D0bar"}; + Configurable cutYCandMax{"cutYCandMax", -1., "max. cand. rapidity"}; + Configurable cutPtCandMin{"cutPtCandMin", -1., "min. cand. pT"}; + Configurable> bins{"ptBinsForMass", std::vector{o2::analysis::hf_cuts_d0_topik::pTBins_v}, "pT bin limits for candidate mass plots"}; + + void init(o2::framework::InitContext&) + { + registry.add("hMass", "D0,D0bar candidates;inv. mass (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{120, 1.5848, 2.1848}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hMassD0", "D0,D0bar candidates;inv. mass D0 only (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{120, 1.5848, 2.1848}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hMassD0bar", "D0,D0bar candidates;inv. mass D0bar only (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{120, 1.5848, 2.1848}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + } + + Filter filterSelectCandidates = (aod::hf_selcandidate_d0::isSelD0 >= selectionFlagD0 || aod::hf_selcandidate_d0::isSelD0bar >= selectionFlagD0bar); + + void process(aod::Collision const& collision, soa::Filtered> const& candidates) + { + for (auto& candidate1 : candidates) { + //check decay channel flag for candidate1 + if (!(candidate1.hfflag() & 1 << DecayType::D0ToPiK)) { + continue; + } + if (cutYCandMax >= 0. && std::abs(YD0(candidate1)) > cutYCandMax) { + continue; + } + if (cutPtCandMin >= 0. && candidate1.pt() < cutPtCandMin) { + continue; + } + //fill invariant mass plots and generic info from all D0/D0bar candidates + if (candidate1.isSelD0() >= selectionFlagD0) { + registry.fill(HIST("hMass"), InvMassD0(candidate1), candidate1.pt()); + registry.fill(HIST("hMassD0"), InvMassD0(candidate1), candidate1.pt()); + } + if (candidate1.isSelD0bar() >= selectionFlagD0bar) { + registry.fill(HIST("hMass"), InvMassD0bar(candidate1), candidate1.pt()); + registry.fill(HIST("hMassD0bar"), InvMassD0bar(candidate1), candidate1.pt()); + } + 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"), YD0(candidate1)); + registry.fill(HIST("hSelectionStatus"), candidate1.isSelD0bar() + (candidate1.isSelD0() * 2)); + + //D-Dbar correlation dedicated section + //For like-sign, first loop on both D0 and D0bars. First candidate is for sure a D0 and D0bars (checked before, so don't re-check anything on it) + for (auto& candidate2 : candidates) { + //check decay channel flag for candidate2 + if (!(candidate2.hfflag() & 1 << DecayType::D0ToPiK)) { + continue; + } + //for the associated, has to have smaller pT, and pass D0sel if trigger passes D0sel, or D0barsel if trigger passes D0barsel + if ((candidate1.isSelD0() >= selectionFlagD0 && candidate2.isSelD0() >= selectionFlagD0) || (candidate1.isSelD0bar() >= selectionFlagD0bar && candidate2.isSelD0bar() >= selectionFlagD0bar)) { + if (cutYCandMax >= 0. && std::abs(YD0(candidate2)) > cutYCandMax) { + continue; + } + if (cutPtCandMin >= 0. && candidate2.pt() < cutPtCandMin) { + continue; + } + //Excluding self-correlations + if (candidate1.mRowIndex == candidate2.mRowIndex) { + continue; + } + entryD0D0barPair(getDeltaPhi(candidate2.phi(), candidate1.phi()), + candidate2.eta() - candidate1.eta(), + candidate1.pt(), + candidate2.pt()); + entryD0D0barRecoInfo(InvMassD0(candidate1), + InvMassD0bar(candidate2), + 0); + } + //note: candidates selected as both D0 and D0bar are used, and considered in both situation (but not auto-correlated): reflections could play a relevant role. + //another, more restrictive, option, could be to consider only candidates selected with a single option (D0 xor D0bar) + } // end inner loop (Dbars) + } //end outer loop + } +}; + +/// D0-D0bar correlation pair builder - LIKE SIGN - for MC reco analysis (data-like but matching to true DO and D0bar) +/// NOTE: At the moment, both dPhi-symmetrical correlation pairs (part1-part2 and part2-part1) are filled, +/// since we bin in pT and selecting as trigger the largest pT particle would bias the distributions w.r.t. the ULS case. +struct HfCorrelatorD0D0barMcRecLs { + + Produces entryD0D0barPair; + Produces entryD0D0barRecoInfo; + + HistogramRegistry registry{ + "registry", + //NOTE: use hMassD0 for trigger normalisation (S*0.955), and hMass2DCorrelationPairs (in final task) for 2D-sideband-subtraction purposes + {{"hPtCandMCRec", "D0,D0bar candidates - MC reco;candidate #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, + {"hPtProng0MCRec", "D0,D0bar candidates - MC reco;prong 0 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, + {"hPtProng1MCRec", "D0,D0bar candidates - MC reco;prong 1 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, + {"hSelectionStatusMCRec", "D0,D0bar candidates - MC reco;selection status;entries", {HistType::kTH1F, {{4, -0.5, 3.5}}}}, + {"hEtaMCRec", "D0,D0bar candidates - MC reco;candidate #it{#eta};entries", {HistType::kTH1F, {{100, -5., 5.}}}}, + {"hPhiMCRec", "D0,D0bar candidates - MC reco;candidate #it{#varphi};entries", {HistType::kTH1F, {{32, 0., 2. * o2::constants::math::PI}}}}, + {"hYMCRec", "D0,D0bar candidates - MC reco;candidate #it{#y};entries", {HistType::kTH1F, {{100, -5., 5.}}}}}}; + + Configurable selectionFlagD0{"selectionFlagD0", 1, "Selection Flag for D0"}; + Configurable selectionFlagD0bar{"selectionFlagD0bar", 1, "Selection Flag for D0bar"}; + Configurable cutYCandMax{"cutYCandMax", -1., "max. cand. rapidity"}; + Configurable cutPtCandMin{"cutPtCandMin", -1., "min. cand. pT"}; + Configurable> bins{"ptBinsForMass", std::vector{o2::analysis::hf_cuts_d0_topik::pTBins_v}, "pT bin limits for candidate mass plots"}; + + void init(o2::framework::InitContext&) + { + registry.add("hMassD0MCRec", "D0,D0bar candidates - MC reco;inv. mass (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{120, 1.5848, 2.1848}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hMassD0barMCRec", "D0,D0bar candidates - MC reco;inv. mass D0 only (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{120, 1.5848, 2.1848}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + } + + Filter filterSelectCandidates = (aod::hf_selcandidate_d0::isSelD0 >= selectionFlagD0 || aod::hf_selcandidate_d0::isSelD0bar >= selectionFlagD0bar); + + void process(aod::Collision const& collision, soa::Filtered> const& candidates) + { + //MC reco level + for (auto& candidate1 : candidates) { + //check decay channel flag for candidate1 + if (!(candidate1.hfflag() & 1 << DecayType::D0ToPiK)) { + continue; + } + if (cutYCandMax >= 0. && std::abs(YD0(candidate1)) > cutYCandMax) { + continue; + } + if (cutPtCandMin >= 0. && candidate1.pt() < cutPtCandMin) { + continue; + } + if (std::abs(candidate1.flagMCMatchRec()) == 1 << DecayType::D0ToPiK) { + //fill invariant mass plots and generic info from all D0/D0bar candidates + if (candidate1.isSelD0() >= selectionFlagD0 && candidate1.flagMCMatchRec() == DecayType::D0ToPiK) { //only reco and matched as D0 + registry.fill(HIST("hMassD0MCRec"), InvMassD0(candidate1)); + } + if (candidate1.isSelD0bar() >= selectionFlagD0bar && candidate1.flagMCMatchRec() == DecayType::D0ToPiK) { //only reco and matched as D0bar + registry.fill(HIST("hMassD0barMCRec"), InvMassD0bar(candidate1)); + } + registry.fill(HIST("hPtCandMCRec"), candidate1.pt()); + registry.fill(HIST("hPtProng0MCRec"), candidate1.ptProng0()); + registry.fill(HIST("hPtProng1MCRec"), candidate1.ptProng1()); + registry.fill(HIST("hEtaMCRec"), candidate1.eta()); + registry.fill(HIST("hPhiMCRec"), candidate1.phi()); + registry.fill(HIST("hYMCRec"), YD0(candidate1)); + registry.fill(HIST("hSelectionStatusMCRec"), candidate1.isSelD0bar() + (candidate1.isSelD0() * 2)); + + //D-Dbar correlation dedicated section + //For like-sign, first loop on both D0 and D0bars. First candidate is for sure a D0 and D0bars (looping on filtered) and was already matched, so don't re-check anything on it) + for (auto& candidate2 : candidates) { + //check decay channel flag for candidate2 + if (!(candidate2.hfflag() & 1 << DecayType::D0ToPiK)) { + continue; + } + bool conditionLSForD0 = (candidate1.isSelD0() >= selectionFlagD0bar && candidate1.flagMCMatchRec() == 1 << DecayType::D0ToPiK) && (candidate2.isSelD0() >= selectionFlagD0bar && candidate2.flagMCMatchRec() == 1 << DecayType::D0ToPiK); + bool conditionLSForD0bar = (candidate1.isSelD0bar() >= selectionFlagD0bar && candidate1.flagMCMatchRec() == -(1 << DecayType::D0ToPiK)) && (candidate2.isSelD0bar() >= selectionFlagD0bar && candidate2.flagMCMatchRec() == -(1 << DecayType::D0ToPiK)); + if (conditionLSForD0 || conditionLSForD0bar) { //LS pair (of D0 or of D0bar) + pt2= 0. && std::abs(YD0(candidate2)) > cutYCandMax) { + continue; + } + if (cutPtCandMin >= 0. && candidate2.pt() < cutPtCandMin) { + continue; + } + //Excluding self-correlations + if (candidate1.mRowIndex == candidate2.mRowIndex) { + continue; + } + entryD0D0barPair(getDeltaPhi(candidate2.phi(), candidate1.phi()), + candidate2.eta() - candidate1.eta(), + candidate1.pt(), + candidate2.pt()); + entryD0D0barRecoInfo(InvMassD0(candidate1), + InvMassD0bar(candidate2), + 0); //for LS studies we set a dummy 0 for pairSignalStatus (there are no more the usual 4 possible combinations) + + } //end inner if (MC match) + } // end inner loop (Dbars) + } //end outer if (MC match) + } //end outer loop + } +}; + +/// D0-D0bar correlation pair builder - for MC gen-level analysis, like sign particles +/// NOTE: At the moment, both dPhi-symmetrical correlation pairs (part1-part2 and part2-part1) are filled, +/// since we bin in pT and selecting as trigger the largest pT particle would bias the distributions w.r.t. the ULS case. +struct HfCorrelatorD0D0barMcGenLs { + + Produces entryD0D0barPair; + + HistogramRegistry registry{ + "registry", + {{"hMCEvtCount", "Event counter - MC gen;;entries", {HistType::kTH1F, {{1, -0.5, 0.5}}}}, + {"hPtCandMCGen", "D0,D0bar particles - MC gen;particle #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, + {"hEtaMCGen", "D0,D0bar particles - MC gen;particle #it{#eta};entries", {HistType::kTH1F, {{100, -5., 5.}}}}, + {"hPhiMCGen", "D0,D0bar particles - MC gen;particle #it{#varphi};entries", {HistType::kTH1F, {{32, 0., 2. * o2::constants::math::PI}}}}, + {"hYMCGen", "D0,D0bar candidates - MC gen;candidate #it{#y};entries", {HistType::kTH1F, {{100, -5., 5.}}}}, + {"hcountD0D0barPerEvent", "D0,D0bar particles - MC gen;Number per event;entries", {HistType::kTH1F, {{20, 0., 20.}}}}}}; + + Configurable cutYCandMax{"cutYCandMax", -1., "max. cand. rapidity"}; + Configurable cutPtCandMin{"cutPtCandMin", -1., "min. cand. pT"}; + Configurable> bins{"ptBinsForMass", std::vector{o2::analysis::hf_cuts_d0_topik::pTBins_v}, "pT bin limits for trigger counters"}; + + void init(o2::framework::InitContext&) + { + registry.add("hcountD0triggersMCGen", "D0 trigger particles - MC gen;;N of trigger D0", {HistType::kTH2F, {{1, -0.5, 0.5}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + } + + void process(aod::McCollision const& mccollision, soa::Join const& particlesMC) + { + int counterD0D0bar = 0; + registry.fill(HIST("hMCEvtCount"), 0); + //MC gen level + for (auto& particle1 : particlesMC) { + //check if the particle is D0 or D0bar (both can be trigger) - NOTE: decay channel is not probed! + if (std::abs(particle1.pdgCode()) != pdg::Code::kD0) { + continue; + } + double yD = RecoDecay::Y(array{particle1.px(), particle1.py(), particle1.pz()}, RecoDecay::getMassPDG(particle1.pdgCode())); + if (cutYCandMax >= 0. && std::abs(yD) > cutYCandMax) { + continue; + } + if (cutPtCandMin >= 0. && particle1.pt() < cutPtCandMin) { + continue; + } + + registry.fill(HIST("hPtCandMCGen"), particle1.pt()); + registry.fill(HIST("hEtaMCGen"), particle1.eta()); + registry.fill(HIST("hPhiMCGen"), particle1.phi()); + registry.fill(HIST("hYMCGen"), yD); + counterD0D0bar++; + //D-Dbar correlation dedicated section + //if it's D0, search for D0bar and evaluate correlations. + registry.fill(HIST("hcountD0triggersMCGen"), 0, particle1.pt()); //to count trigger D0 (normalisation) + for (auto& particle2 : particlesMC) { + if (std::abs(particle2.pdgCode()) != pdg::Code::kD0) { //check that associated is a D0/D0bar (both are fine) + continue; + } + if (cutYCandMax >= 0. && std::abs(RecoDecay::Y(array{particle2.px(), particle2.py(), particle2.pz()}, RecoDecay::getMassPDG(particle2.pdgCode()))) > cutYCandMax) { + continue; + } + if (cutPtCandMin >= 0. && particle2.pt() < cutPtCandMin) { + continue; + } + if (particle2.pdgCode() == particle1.pdgCode()) { //like-sign condition (both 421 or both -421) and pT_Trig>pT_assoc + //Excluding self-correlations + if (particle1.mRowIndex == particle2.mRowIndex) { + continue; + } + entryD0D0barPair(getDeltaPhi(particle2.phi(), particle1.phi()), + particle2.eta() - particle1.eta(), + particle1.pt(), + particle2.pt()); + } + } // end inner loop (Dbars) + } //end outer loop + registry.fill(HIST("hcountD0D0barPerEvent"), counterD0D0bar); + } +}; + +/// c-cbar correlator table builder - for MC gen-level analysis +struct HfCorrelatorCCbarMcGen { + + Produces entryccbarPair; + + HistogramRegistry registry{ + "registry", + {{"hMCEvtCount", "Event counter - MC gen;;entries", {HistType::kTH1F, {{1, -0.5, 0.5}}}}, + {"hPtCandMCGen", "c,cbar particles - MC gen;particle #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, + {"hEtaMCGen", "c,cbar particles - MC gen;particle #it{#eta};entries", {HistType::kTH1F, {{100, -5., 5.}}}}, + {"hYMCGen", "c,cbar candidates - MC gen;candidate #it{#y};entries", {HistType::kTH1F, {{100, -5., 5.}}}}, + {"hPhiMCGen", "c,cbar particles - MC gen;particle #it{#varphi};entries", {HistType::kTH1F, {{32, 0., 2. * o2::constants::math::PI}}}}, + {"hcountCCbarPerEvent", "c,cbar particles - MC gen;Number per event;entries", {HistType::kTH1F, {{20, 0., 20.}}}}, + {"hcountCCbarPerEventPreEtaCut", "c,cbar particles - MC gen;Number per event pre #eta cut;entries", {HistType::kTH1F, {{20, 0., 20.}}}}}}; + + Configurable cutYCandMax{"cutYCandMax", -1., "max. cand. rapidity"}; + Configurable cutPtCandMin{"cutPtCandMin", -1., "min. cand. pT"}; + Configurable> bins{"ptBinsForMass", std::vector{o2::analysis::hf_cuts_d0_topik::pTBins_v}, "pT bin limits for trigger counters"}; + + void init(o2::framework::InitContext&) + { + registry.add("hcountCtriggersMCGen", "c trigger particles - MC gen;;N of trigger c quark", {HistType::kTH2F, {{1, -0.5, 0.5}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + } + + void process(aod::McCollision const& mccollision, soa::Join const& particlesMC) + { + registry.fill(HIST("hMCEvtCount"), 0); + int counterccbar = 0, counterccbarPreEtasel = 0; + + //loop over particles at MC gen level + for (auto& particle1 : particlesMC) { + if (std::abs(particle1.pdgCode()) != PDG_t::kCharm) { //search c or cbar particles + continue; + } + int partMothPDG = particlesMC.iteratorAt(particle1.mother0()).pdgCode(); + //check whether mothers of quark c/cbar are still '4'/'-4' particles - in that case the c/cbar quark comes from its own fragmentation, skip it + if (partMothPDG == particle1.pdgCode()) { + continue; + } + counterccbarPreEtasel++; //count c or cbar (before kinematic selection) + double yC = RecoDecay::Y(array{particle1.px(), particle1.py(), particle1.pz()}, RecoDecay::getMassPDG(particle1.pdgCode())); + if (cutYCandMax >= 0. && std::abs(yC) > cutYCandMax) { + continue; + } + if (cutPtCandMin >= 0. && particle1.pt() < cutPtCandMin) { + continue; + } + registry.fill(HIST("hPtCandMCGen"), particle1.pt()); + registry.fill(HIST("hEtaMCGen"), particle1.eta()); + registry.fill(HIST("hPhiMCGen"), particle1.phi()); + registry.fill(HIST("hYMCGen"), yC); + counterccbar++; //count if c or cbar don't come from themselves during fragmentation (after kinematic selection) + + //c-cbar correlation dedicated section + //if it's c, search for cbar and evaluate correlations. + if (particle1.pdgCode() != PDG_t::kCharm) { + continue; + } + registry.fill(HIST("hcountCtriggersMCGen"), 0, particle1.pt()); //to count trigger c quark (for normalisation) + + for (auto& particle2 : particlesMC) { + if (particle2.pdgCode() != PDG_t::kCharmBar) { //check that inner particle is a cbar + continue; + } + if (cutYCandMax >= 0. && std::abs(RecoDecay::Y(array{particle2.px(), particle2.py(), particle2.pz()}, RecoDecay::getMassPDG(particle2.pdgCode()))) > cutYCandMax) { + continue; + } + if (cutPtCandMin >= 0. && particle2.pt() < cutPtCandMin) { + continue; + } + //check whether mothers of quark cbar (from associated loop) are still '-4' particles - in that case the cbar quark comes from its own fragmentation, skip it + if (particlesMC.iteratorAt(particle2.mother0()).pdgCode() == PDG_t::kCharmBar) { + continue; + } + entryccbarPair(getDeltaPhi(particle2.phi(), particle1.phi()), + particle2.eta() - particle1.eta(), + particle1.pt(), + particle2.pt()); + } // end inner loop + } //end outer loop + registry.fill(HIST("hcountCCbarPerEvent"), counterccbar); + registry.fill(HIST("hcountCCbarPerEventPreEtaCut"), counterccbarPreEtasel); + } +}; + +/// c-cbar correlator table builder - for MC gen-level analysis - Like Sign +struct HfCorrelatorCCbarMcGenLs { + + Produces entryccbarPair; + + HistogramRegistry registry{ + "registry", + {{"hMCEvtCount", "Event counter - MC gen;;entries", {HistType::kTH1F, {{1, -0.5, 0.5}}}}, + {"hPtCandMCGen", "c,cbar particles - MC gen;particle #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, + {"hEtaMCGen", "c,cbar particles - MC gen;particle #it{#eta};entries", {HistType::kTH1F, {{100, -5., 5.}}}}, + {"hYMCGen", "c,cbar candidates - MC gen;candidate #it{#y};entries", {HistType::kTH1F, {{100, -5., 5.}}}}, + {"hPhiMCGen", "c,cbar particles - MC gen;particle #it{#varphi};entries", {HistType::kTH1F, {{32, 0., 2. * o2::constants::math::PI}}}}, + {"hcountCCbarPerEvent", "c,cbar particles - MC gen;Number per event;entries", {HistType::kTH1F, {{20, 0., 20.}}}}, + {"hcountCCbarPerEventPreEtaCut", "c,cbar particles - MC gen;Number per event pre #eta cut;entries", {HistType::kTH1F, {{20, 0., 20.}}}}}}; + + Configurable cutYCandMax{"cutYCandMax", -1., "max. cand. rapidity"}; + Configurable cutPtCandMin{"cutPtCandMin", -1., "min. cand. pT"}; + Configurable> bins{"ptBinsForMass", std::vector{o2::analysis::hf_cuts_d0_topik::pTBins_v}, "pT bin limits for trigger counters"}; + + void init(o2::framework::InitContext&) + { + registry.add("hcountCtriggersMCGen", "c trigger particles - MC gen;;N of trigger c quark", {HistType::kTH2F, {{1, -0.5, 0.5}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + } + + void process(aod::McCollision const& mccollision, soa::Join const& particlesMC) + { + registry.fill(HIST("hMCEvtCount"), 0); + int counterccbar = 0, counterccbarPreEtasel = 0; + + //loop over particles at MC gen level + for (auto& particle1 : particlesMC) { + if (std::abs(particle1.pdgCode()) != PDG_t::kCharm) { //search c or cbar particles + continue; + } + int partMothPDG = particlesMC.iteratorAt(particle1.mother0()).pdgCode(); + //check whether mothers of quark c/cbar are still '4'/'-4' particles - in that case the c/cbar quark comes from its own fragmentation, skip it + if (partMothPDG == particle1.pdgCode()) { + continue; + } + counterccbarPreEtasel++; //count c or cbar (before kinematic selection) + double yC = RecoDecay::Y(array{particle1.px(), particle1.py(), particle1.pz()}, RecoDecay::getMassPDG(particle1.pdgCode())); + if (cutYCandMax >= 0. && std::abs(yC) > cutYCandMax) { + continue; + } + if (cutPtCandMin >= 0. && particle1.pt() < cutPtCandMin) { + continue; + } + registry.fill(HIST("hPtCandMCGen"), particle1.pt()); + registry.fill(HIST("hEtaMCGen"), particle1.eta()); + registry.fill(HIST("hPhiMCGen"), particle1.phi()); + registry.fill(HIST("hYMCGen"), yC); + counterccbar++; //count if c or cbar don't come from themselves during fragmentation (after kinematic selection) + + //c-cbar correlation dedicated section + registry.fill(HIST("hcountCtriggersMCGen"), 0, particle1.pt()); //to count trigger c quark (for normalisation) + + for (auto& particle2 : particlesMC) { + if (std::abs(particle2.pdgCode()) != PDG_t::kCharm) { //search c or cbar for associated particles + continue; + } + if (cutYCandMax >= 0. && std::abs(RecoDecay::Y(array{particle2.px(), particle2.py(), particle2.pz()}, RecoDecay::getMassPDG(particle2.pdgCode()))) > cutYCandMax) { + continue; + } + if (cutPtCandMin >= 0. && particle2.pt() < cutPtCandMin) { + continue; + } + if (particle2.pdgCode() == particle1.pdgCode()) { + //check whether mothers of quark cbar (from associated loop) are still '-4' particles - in that case the cbar quark comes from its own fragmentation, skip it + if (particlesMC.iteratorAt(particle2.mother0()).pdgCode() == particle2.pdgCode()) { + continue; + } + //Excluding self-correlations + if (particle1.mRowIndex == particle2.mRowIndex) { + continue; + } + entryccbarPair(getDeltaPhi(particle2.phi(), particle1.phi()), + particle2.eta() - particle1.eta(), + particle1.pt(), + particle2.pt()); + } // end outer if (check PDG associate) + } // end inner loop + } //end outer loop + registry.fill(HIST("hcountCCbarPerEvent"), counterccbar); + registry.fill(HIST("hcountCCbarPerEventPreEtaCut"), counterccbarPreEtasel); + } +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + WorkflowSpec workflow{}; + const bool doMCccbar = cfgc.options().get("doMCccbar"); + const bool doMCGen = cfgc.options().get("doMCGen"); + const bool doMCRec = cfgc.options().get("doMCRec"); + const bool doLikeSign = cfgc.options().get("doLikeSign"); + if (!doLikeSign) { //unlike-sign analyses + if (doMCGen) { //MC-Gen analysis + workflow.push_back(adaptAnalysisTask(cfgc)); + } else if (doMCRec) { //MC-Reco analysis + workflow.push_back(adaptAnalysisTask(cfgc)); + } else if (doMCccbar) { //MC-Reco analysis + workflow.push_back(adaptAnalysisTask(cfgc)); + } else { //data analysis + workflow.push_back(adaptAnalysisTask(cfgc)); + } + } else { //like-sign analyses + if (doMCGen) { //MC-Gen analysis + workflow.push_back(adaptAnalysisTask(cfgc)); + } else if (doMCRec) { //MC-Reco analysis + workflow.push_back(adaptAnalysisTask(cfgc)); + } else if (doMCccbar) { //MC-Reco analysis + workflow.push_back(adaptAnalysisTask(cfgc)); + } else { //data analysis + workflow.push_back(adaptAnalysisTask(cfgc)); + } + } + + return workflow; +} diff --git a/Analysis/Tasks/PWGHF/HFCorrelatorDplusDminus.cxx b/Analysis/Tasks/PWGHF/HFCorrelatorDplusDminus.cxx new file mode 100644 index 0000000000000..2a52f4d3bebf0 --- /dev/null +++ b/Analysis/Tasks/PWGHF/HFCorrelatorDplusDminus.cxx @@ -0,0 +1,890 @@ +// Copyright CERN and copyright holders of ALICE O2. This software is +// distributed under the terms of the GNU General Public License v3 (GPL +// Version 3), copied verbatim in the file "COPYING". +// +// See http://alice-o2.web.cern.ch/license for full licensing information. +// +// 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 HFCorrelatorDplusDminus.cxx +/// \brief Dplus-Dminus correlator task - data-like, MC-reco and MC-kine analyses. For ULS and LS pairs +/// +/// \author Fabio Colamaria , INFN Bari + +#include "Framework/AnalysisTask.h" +#include "Framework/HistogramRegistry.h" +#include "AnalysisDataModel/HFSecondaryVertex.h" +#include "AnalysisDataModel/HFCandidateSelectionTables.h" + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; +using namespace o2::aod::hf_cand_prong3; +using namespace o2::aod::hf_correlation_ddbar; +using namespace o2::analysis::hf_cuts_dplus_topikpi; +using namespace o2::constants::math; + +void customize(std::vector& workflowOptions) +{ + ConfigParamSpec optionDoLikeSign{"doLikeSign", VariantType::Bool, false, {"Run Like-Sign analysis."}}; + ConfigParamSpec optionDoMCccbar{"doMCccbar", VariantType::Bool, false, {"Run MC-Gen dedicated tasks."}}; + ConfigParamSpec optionDoMCGen{"doMCGen", VariantType::Bool, false, {"Run MC-Gen dedicated tasks."}}; + ConfigParamSpec optionDoMCRec{"doMCRec", VariantType::Bool, false, {"Run MC-Rec dedicated tasks."}}; + workflowOptions.push_back(optionDoLikeSign); + workflowOptions.push_back(optionDoMCccbar); + workflowOptions.push_back(optionDoMCGen); + workflowOptions.push_back(optionDoMCRec); +} + +#include "Framework/runDataProcessing.h" + +/// +/// Returns deltaPhi value in range [-pi/2., 3.*pi/2], typically used for correlation studies +/// +double getDeltaPhi(double phiD, double phiDbar) +{ + return RecoDecay::constrainAngle(phiDbar - phiD, -o2::constants::math::PI / 2.); +} + +/// definition of variables for DplusDminus pairs vs eta acceptance studies (hDDbarVsEtaCut, in data-like, MC-reco and MC-kine tasks) +const double maxEtaCut = 5.; +const double ptThresholdForMaxEtaCut = 10.; +const double incrementEtaCut = 0.1; +const double incrementPtThreshold = 0.5; +const double epsilon = 1E-5; + +/// Dplus-Dminus correlation pair builder - for real data and data-like analysis (i.e. reco-level w/o matching request via MC truth) +struct HfCorrelatorDplusDminus { + Produces entryDplusDminusPair; + Produces entryDplusDminusRecoInfo; + + HistogramRegistry registry{ + "registry", + //NOTE: use hMassDplus for trigger normalisation (S*0.955), and hMass2DCorrelationPairs (in final task) for 2D-sideband-subtraction purposes + {{"hPtCand", "Dplus,Dminus candidates;candidate #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, + {"hPtProng0", "Dplus,Dminus candidates;prong 0 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, + {"hPtProng1", "Dplus,Dminus candidates;prong 1 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, + {"hPtProng2", "Dplus,Dminus candidates;prong 2 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, + {"hSelectionStatus", "Dplus,Dminus candidates;selection status;entries", {HistType::kTH1F, {{2, -0.5, 1.5}}}}, + {"hEta", "Dplus,Dminus candidates;candidate #it{#eta};entries", {HistType::kTH1F, {{100, -5., 5.}}}}, + {"hPhi", "Dplus,Dminus candidates;candidate #it{#varphi};entries", {HistType::kTH1F, {{32, 0., 2. * o2::constants::math::PI}}}}, + {"hY", "Dplus,Dminus candidates;candidate #it{#y};entries", {HistType::kTH1F, {{100, -5., 5.}}}}, + {"hDDbarVsEtaCut", "Dplus,Dminus pairs vs #eta cut;#eta_{max};entries", {HistType::kTH2F, {{(int)(maxEtaCut / incrementEtaCut), 0., maxEtaCut}, {(int)(ptThresholdForMaxEtaCut / incrementPtThreshold), 0., ptThresholdForMaxEtaCut}}}}}}; + + Configurable selectionFlagDplus{"selectionFlagDplus", 1, "Selection Flag for Dplus,Dminus"}; + Configurable cutYCandMax{"cutYCandMax", -1., "max. cand. rapidity"}; + Configurable cutPtCandMin{"cutPtCandMin", -1., "min. cand. pT"}; + Configurable> bins{"ptBinsForMass", std::vector{o2::analysis::hf_cuts_dplus_topikpi::pTBins_v}, "pT bin limits for candidate mass plots"}; + + void init(o2::framework::InitContext&) + { + registry.add("hMass", "Dplus,Dminus candidates;inv. mass (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{120, 1.5848, 2.1848}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hMassDplus", "Dplus,Dminus candidates;inv. mass Dplus only (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{120, 1.5848, 2.1848}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hMassDminus", "Dplus,Dminus candidates;inv. mass Dminus only (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{120, 1.5848, 2.1848}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + } + + Filter filterSelectCandidates = (aod::hf_selcandidate_dplus::isSelDplusToPiKPi >= selectionFlagDplus); + + void process(aod::Collision const& collision, soa::Filtered> const& candidates, aod::BigTracks const& tracks) + { + for (auto& candidate1 : candidates) { + if (cutYCandMax >= 0. && std::abs(YDPlus(candidate1)) > cutYCandMax) { + continue; + } + if (cutPtCandMin >= 0. && candidate1.pt() < cutPtCandMin) { + continue; + } + //check decay channel flag for candidate1 + if (!(candidate1.hfflag() & 1 << DecayType::DPlusToPiKPi)) { //probably dummy since already selected? not sure... + continue; + } + int outerParticleSign = 1; //Dplus + auto outerSecondTrack = candidate1.index1_as(); + if (outerSecondTrack.sign() == 1) { + outerParticleSign = -1; //Dminus (second daughter track is positive) + } + + //fill invariant mass plots and generic info from all Dplus/Dminus candidates + if (outerParticleSign == 1) { + registry.fill(HIST("hMass"), InvMassDPlus(candidate1), candidate1.pt()); + registry.fill(HIST("hMassDplus"), InvMassDPlus(candidate1), candidate1.pt()); + } else { + registry.fill(HIST("hMass"), InvMassDPlus(candidate1), candidate1.pt()); + registry.fill(HIST("hMassDminus"), InvMassDPlus(candidate1), candidate1.pt()); + } + registry.fill(HIST("hPtCand"), candidate1.pt()); + registry.fill(HIST("hPtProng0"), candidate1.ptProng0()); + registry.fill(HIST("hPtProng1"), candidate1.ptProng1()); + registry.fill(HIST("hPtProng2"), candidate1.ptProng2()); + registry.fill(HIST("hEta"), candidate1.eta()); + registry.fill(HIST("hPhi"), candidate1.phi()); + registry.fill(HIST("hY"), YDPlus(candidate1)); + registry.fill(HIST("hSelectionStatus"), candidate1.isSelDplusToPiKPi()); + + //D-Dbar correlation dedicated section + //if the candidate is a Dplus, search for Dminus and evaluate correlations + if (outerParticleSign != 1) { + continue; + } + for (auto& candidate2 : candidates) { + //check decay channel flag for candidate2 + if (!(candidate2.hfflag() & 1 << DecayType::DPlusToPiKPi)) { //probably dummy since already selected? not sure... + continue; + } + auto innerSecondTrack = candidate2.index1_as(); + if (innerSecondTrack.sign() != 1) { //keep only Dminus (with second daughter track positive) + continue; + } + if (cutYCandMax >= 0. && std::abs(YDPlus(candidate2)) > cutYCandMax) { + continue; + } + if (cutPtCandMin >= 0. && candidate2.pt() < cutPtCandMin) { + continue; + } + entryDplusDminusPair(getDeltaPhi(candidate2.phi(), candidate1.phi()), + candidate2.eta() - candidate1.eta(), + candidate1.pt(), + candidate2.pt()); + entryDplusDminusRecoInfo(InvMassDPlus(candidate1), + InvMassDPlus(candidate2), + 0); + double etaCut = 0.; + double ptCut = 0.; + do { //fill pairs vs etaCut plot + ptCut = 0.; + etaCut += incrementEtaCut; + do { //fill pairs vs etaCut plot + if (std::abs(candidate1.eta()) < etaCut && std::abs(candidate2.eta()) < etaCut && candidate1.pt() > ptCut && candidate2.pt() > ptCut) { + registry.fill(HIST("hDDbarVsEtaCut"), etaCut - epsilon, ptCut + epsilon); + } + ptCut += incrementPtThreshold; + } while (ptCut < ptThresholdForMaxEtaCut - epsilon); + } while (etaCut < maxEtaCut - epsilon); + } // end inner loop (Dminus) + } //end outer loop (Dplus) + } +}; + +/// Dplus-Dminus correlation pair builder - for MC reco-level analysis (candidates matched to true signal only, but also the various bkg sources are studied) +struct HfCorrelatorDplusDminusMcRec { + Produces entryDplusDminusPair; + Produces entryDplusDminusRecoInfo; + + HistogramRegistry registry{ + "registry", + //NOTE: use hMassDplus for trigger normalisation (S*0.955), and hMass2DCorrelationPairs (in final task) for 2D-sideband-subtraction purposes + {{"hPtCandMCRec", "Dplus,Dminus candidates - MC reco;candidate #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, + {"hPtProng0MCRec", "Dplus,Dminus candidates - MC reco;prong 0 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, + {"hPtProng1MCRec", "Dplus,Dminus candidates - MC reco;prong 1 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, + {"hPtProng2MCRec", "Dplus,Dminus candidates - MC reco;prong 2 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, + {"hSelectionStatusMCRec", "Dplus,Dminus candidates - MC reco;selection status;entries", {HistType::kTH1F, {{2, -0.5, 1.5}}}}, + {"hEtaMCRec", "Dplus,Dminus candidates - MC reco;candidate #it{#eta};entries", {HistType::kTH1F, {{100, -5., 5.}}}}, + {"hPhiMCRec", "Dplus,Dminus candidates - MC reco;candidate #it{#varphi};entries", {HistType::kTH1F, {{32, 0., 2. * o2::constants::math::PI}}}}, + {"hYMCRec", "Dplus,Dminus candidates - MC reco;candidate #it{#y};entries", {HistType::kTH1F, {{100, -5., 5.}}}}, + {"hDDbarVsEtaCut", "Dplus,Dminus pairs vs #eta cut;#eta_{max};entries", {HistType::kTH2F, {{(int)(maxEtaCut / incrementEtaCut), 0., maxEtaCut}, {(int)(ptThresholdForMaxEtaCut / incrementPtThreshold), 0., ptThresholdForMaxEtaCut}}}}}}; + + Configurable selectionFlagDplus{"selectionFlagDplus", 1, "Selection Flag for Dplus,Dminus"}; + Configurable cutYCandMax{"cutYCandMax", -1., "max. cand. rapidity"}; + Configurable cutPtCandMin{"cutPtCandMin", -1., "min. cand. pT"}; + Configurable> bins{"ptBinsForMass", std::vector{o2::analysis::hf_cuts_dplus_topikpi::pTBins_v}, "pT bin limits for candidate mass plots"}; + + void init(o2::framework::InitContext&) + { + registry.add("hMassDplusMCRec", "Dplus,Dminus candidates - MC reco;inv. mass (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{120, 1.5848, 2.1848}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hMassDminusMCRec", "Dplus,Dminus candidates - MC reco;inv. mass D0 only (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{120, 1.5848, 2.1848}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + } + + Filter filterSelectCandidates = (aod::hf_selcandidate_dplus::isSelDplusToPiKPi >= selectionFlagDplus); + + void process(aod::Collision const& collision, soa::Filtered> const& candidates, aod::BigTracks const& tracks) + { + //MC reco level + bool flagDplusSignal = false; + bool flagDminusSignal = false; + for (auto& candidate1 : candidates) { + //check decay channel flag for candidate1 + if (!(candidate1.hfflag() & 1 << DecayType::DPlusToPiKPi)) { + continue; + } + if (cutYCandMax >= 0. && std::abs(YDPlus(candidate1)) > cutYCandMax) { + continue; + } + if (cutPtCandMin >= 0. && candidate1.pt() < cutPtCandMin) { + continue; + } + int outerParticleSign = 1; //Dplus + auto outerSecondTrack = candidate1.index1_as(); + if (outerSecondTrack.sign() == 1) { + outerParticleSign = -1; //Dminus (second daughter track is positive) + } + if (std::abs(candidate1.flagMCMatchRec()) == 1 << DecayType::DPlusToPiKPi) { + //fill invariant mass plots and generic info from all Dplus/Dminus candidates + if (outerParticleSign == 1) { //matched as Dplus + registry.fill(HIST("hMassDplusMCRec"), InvMassDPlus(candidate1), candidate1.pt()); + } else { //matched as Dminus + registry.fill(HIST("hMassDminusMCRec"), InvMassDPlus(candidate1), candidate1.pt()); + } + registry.fill(HIST("hPtCandMCRec"), candidate1.pt()); + registry.fill(HIST("hPtProng0MCRec"), candidate1.ptProng0()); + registry.fill(HIST("hPtProng1MCRec"), candidate1.ptProng1()); + registry.fill(HIST("hPtProng2MCRec"), candidate1.ptProng2()); + registry.fill(HIST("hEtaMCRec"), candidate1.eta()); + registry.fill(HIST("hPhiMCRec"), candidate1.phi()); + registry.fill(HIST("hYMCRec"), YDPlus(candidate1)); + registry.fill(HIST("hSelectionStatusMCRec"), candidate1.isSelDplusToPiKPi()); + } + + //D-Dbar correlation dedicated section + if (outerParticleSign == -1) { + continue; //reject Dminus in outer loop + } + flagDplusSignal = std::abs(candidate1.flagMCMatchRec()) == 1 << DecayType::DPlusToPiKPi; //flagDplusSignal 'true' if candidate1 matched to Dplus + for (auto& candidate2 : candidates) { + if (!(candidate2.hfflag() & 1 << DecayType::DPlusToPiKPi)) { //check decay channel flag for candidate2 + continue; + } + auto innerSecondTrack = candidate2.index1_as(); + if (innerSecondTrack.sign() != 1) { //keep only Dminus (with second daughter track positive) + continue; + } + flagDminusSignal = std::abs(candidate2.flagMCMatchRec()) == 1 << DecayType::DPlusToPiKPi; //flagDminusSignal 'true' if candidate2 matched to Dminus + if (cutYCandMax >= 0. && std::abs(YDPlus(candidate2)) > cutYCandMax) { + continue; + } + if (cutPtCandMin >= 0. && candidate2.pt() < cutPtCandMin) { + continue; + } + //choice of options (Dplus/Dminus signal/bkg) + int pairSignalStatus = 0; //0 = bkg/bkg, 1 = bkg/sig, 2 = sig/bkg, 3 = sig/sig + if (flagDplusSignal) { + pairSignalStatus += 2; + } + if (flagDminusSignal) { + pairSignalStatus += 1; + } + entryDplusDminusPair(getDeltaPhi(candidate2.phi(), candidate1.phi()), + candidate2.eta() - candidate1.eta(), + candidate1.pt(), + candidate2.pt()); + entryDplusDminusRecoInfo(InvMassDPlus(candidate1), + InvMassDPlus(candidate2), + pairSignalStatus); + double etaCut = 0.; + double ptCut = 0.; + do { //fill pairs vs etaCut plot + ptCut = 0.; + etaCut += incrementEtaCut; + do { //fill pairs vs etaCut plot + if (std::abs(candidate1.eta()) < etaCut && std::abs(candidate2.eta()) < etaCut && candidate1.pt() > ptCut && candidate2.pt() > ptCut) { + registry.fill(HIST("hDDbarVsEtaCut"), etaCut - epsilon, ptCut + epsilon); + } + ptCut += incrementPtThreshold; + } while (ptCut < ptThresholdForMaxEtaCut - epsilon); + } while (etaCut < maxEtaCut - epsilon); + } // end inner loop (Dbars) + + } //end outer loop + } +}; + +/// Dplus-Dminus correlation pair builder - for MC gen-level analysis (no filter/selection, only true signal) +struct HfCorrelatorDplusDminusMcGen { + + Produces entryDplusDminusPair; + + HistogramRegistry registry{ + "registry", + {{"hMCEvtCount", "Event counter - MC gen;;entries", {HistType::kTH1F, {{1, -0.5, 0.5}}}}, + {"hPtCandMCGen", "Dplus,Dminus particles - MC gen;particle #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, + {"hEtaMCGen", "Dplus,Dminus particles - MC gen;particle #it{#eta};entries", {HistType::kTH1F, {{100, -5., 5.}}}}, + {"hPhiMCGen", "Dplus,Dminus particles - MC gen;particle #it{#varphi};entries", {HistType::kTH1F, {{32, 0., 2. * o2::constants::math::PI}}}}, + {"hYMCGen", "Dplus,Dminus candidates - MC gen;candidate #it{#y};entries", {HistType::kTH1F, {{100, -5., 5.}}}}, + {"hcountDplusDminusPerEvent", "Dplus,Dminus particles - MC gen;Number per event;entries", {HistType::kTH1F, {{20, 0., 20.}}}}, + {"hDDbarVsEtaCut", "Dplus,Dminus pairs vs #eta cut of D mesons;#eta_{max};entries", {HistType::kTH2F, {{(int)(maxEtaCut / incrementEtaCut), 0., maxEtaCut}, {(int)(ptThresholdForMaxEtaCut / incrementPtThreshold), 0., ptThresholdForMaxEtaCut}}}}, + {"hDDbarVsDaughterEtaCut", "Dplus,Dminus pairs vs #eta cut on D daughters;#eta_{max};entries", {HistType::kTH2F, {{(int)(maxEtaCut / incrementEtaCut), 0., maxEtaCut}, {(int)(ptThresholdForMaxEtaCut / incrementPtThreshold), 0., ptThresholdForMaxEtaCut}}}}}}; + + Configurable cutYCandMax{"cutYCandMax", -1., "max. cand. rapidity"}; + Configurable cutPtCandMin{"cutPtCandMin", -1., "min. cand. pT"}; + Configurable> bins{"ptBinsForMass", std::vector{o2::analysis::hf_cuts_dplus_topikpi::pTBins_v}, "pT bin limits for candidate mass plots"}; + + void init(o2::framework::InitContext&) + { + registry.add("hcountDplustriggersMCGen", "Dplus trigger particles - MC gen;;N of trigger D0", {HistType::kTH2F, {{1, -0.5, 0.5}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + } + + void process(aod::McCollision const& mccollision, soa::Join const& particlesMC) + { + int counterDplusDminus = 0; + registry.fill(HIST("hMCEvtCount"), 0); + //MC gen level + for (auto& particle1 : particlesMC) { + //check if the particle is Dplus or Dminus (for general plot filling and selection, so both cases are fine) - NOTE: decay channel is not probed! + if (std::abs(particle1.pdgCode()) != pdg::Code::kDPlus) { + continue; + } + double yD = RecoDecay::Y(array{particle1.px(), particle1.py(), particle1.pz()}, RecoDecay::getMassPDG(particle1.pdgCode())); + if (cutYCandMax >= 0. && std::abs(yD) > cutYCandMax) { + continue; + } + if (cutPtCandMin >= 0. && particle1.pt() < cutPtCandMin) { + continue; + } + registry.fill(HIST("hPtCandMCGen"), particle1.pt()); + registry.fill(HIST("hEtaMCGen"), particle1.eta()); + registry.fill(HIST("hPhiMCGen"), particle1.phi()); + registry.fill(HIST("hYMCGen"), yD); + counterDplusDminus++; + + //D-Dbar correlation dedicated section + //if it's a Dplus particle, search for Dminus and evaluate correlations + if (particle1.pdgCode() != pdg::Code::kDPlus) { //just checking the particle PDG, not the decay channel (differently from Reco: you have a BR factor btw such levels!) + continue; + } + registry.fill(HIST("hcountDplustriggersMCGen"), 0, particle1.pt()); //to count trigger Dplus (for normalisation) + for (auto& particle2 : particlesMC) { + if (particle2.pdgCode() != -pdg::Code::kDPlus) { //check that inner particle is a Dminus + continue; + } + if (cutYCandMax >= 0. && std::abs(RecoDecay::Y(array{particle2.px(), particle2.py(), particle2.pz()}, RecoDecay::getMassPDG(particle2.pdgCode()))) > cutYCandMax) { + continue; + } + if (cutPtCandMin >= 0. && particle2.pt() < cutPtCandMin) { + continue; + } + entryDplusDminusPair(getDeltaPhi(particle2.phi(), particle1.phi()), + particle2.eta() - particle1.eta(), + particle1.pt(), + particle2.pt()); + double etaCut = 0.; + double ptCut = 0.; + + //fill pairs vs etaCut plot + bool rightDecayChannels = false; + if ((std::abs(particle1.flagMCMatchGen()) == 1 << DecayType::DPlusToPiKPi) && (std::abs(particle2.flagMCMatchGen()) == 1 << DecayType::DPlusToPiKPi)) { + rightDecayChannels = true; + } + do { + ptCut = 0.; + etaCut += incrementEtaCut; + do { //fill pairs vs etaCut plot + if (std::abs(particle1.eta()) < etaCut && std::abs(particle2.eta()) < etaCut && particle1.pt() > ptCut && particle2.pt() > ptCut) { + registry.fill(HIST("hDDbarVsEtaCut"), etaCut - epsilon, ptCut + epsilon); + } + if (rightDecayChannels) { //fill with D and Dbar daughter particls acceptance checks + double etaCandidate1Daughter1 = particlesMC.iteratorAt(particle1.daughter0()).eta(); + double etaCandidate1Daughter2 = particlesMC.iteratorAt(particle1.daughter0() + 1).eta(); + double etaCandidate1Daughter3 = particlesMC.iteratorAt(particle1.daughter0() + 2).eta(); + double etaCandidate2Daughter1 = particlesMC.iteratorAt(particle2.daughter0()).eta(); + double etaCandidate2Daughter2 = particlesMC.iteratorAt(particle2.daughter0() + 1).eta(); + double etaCandidate2Daughter3 = particlesMC.iteratorAt(particle2.daughter0() + 2).eta(); + if (std::abs(etaCandidate1Daughter1) < etaCut && std::abs(etaCandidate1Daughter2) < etaCut && std::abs(etaCandidate1Daughter3) < etaCut && + std::abs(etaCandidate2Daughter1) < etaCut && std::abs(etaCandidate2Daughter2) < etaCut && std::abs(etaCandidate2Daughter3) < etaCut && + particle1.pt() > ptCut && particle2.pt() > ptCut) { + registry.fill(HIST("hDDbarVsDaughterEtaCut"), etaCut - epsilon, ptCut + epsilon); + } + } + ptCut += incrementPtThreshold; + } while (ptCut < ptThresholdForMaxEtaCut - epsilon); + } while (etaCut < maxEtaCut - epsilon); + } //end inner loop + } //end outer loop + registry.fill(HIST("hcountDplusDminusPerEvent"), counterDplusDminus); + } +}; + +/// Dplus-Dminus correlation pair builder - LIKE SIGN - for real data and data-like analysis (i.e. reco-level w/o matching request via MC truth) +/// NOTE: At the moment, both dPhi-symmetrical correlation pairs (part1-part2 and part2-part1) are filled, since we bin in pT and selecting as trigger the largest pT particle would bias the distributions w.r.t. the ULS case. +struct HfCorrelatorDplusDminusLs { + + Produces entryDplusDminusPair; + Produces entryDplusDminusRecoInfo; + + HistogramRegistry registry{ + "registry", + //NOTE: use hMassDplus for trigger normalisation (S*0.955), and hMass2DCorrelationPairs (in final task) for 2D-sideband-subtraction purposes + {{"hPtCand", "Dplus,Dminus candidates;candidate #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, + {"hPtProng0", "Dplus,Dminus candidates;prong 0 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, + {"hPtProng1", "Dplus,Dminus candidates;prong 1 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, + {"hPtProng2", "Dplus,Dminus candidates;prong 2 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, + {"hSelectionStatus", "Dplus,Dminus candidates;selection status;entries", {HistType::kTH1F, {{2, -0.5, 1.5}}}}, + {"hEta", "Dplus,Dminus candidates;candidate #it{#eta};entries", {HistType::kTH1F, {{100, -5., 5.}}}}, + {"hPhi", "Dplus,Dminus candidates;candidate #it{#varphi};entries", {HistType::kTH1F, {{32, 0., 2. * o2::constants::math::PI}}}}, + {"hY", "Dplus,Dminus candidates;candidate #it{#y};entries", {HistType::kTH1F, {{100, -5., 5.}}}}}}; + + Configurable selectionFlagDplus{"dSelectivonFlagDplus", 1, "Selection Flag for Dplus,Dminus"}; + Configurable cutYCandMax{"cutYCandMax", -1., "max. cand. rapidity"}; + Configurable cutPtCandMin{"cutPtCandMin", -1., "min. cand. pT"}; + Configurable> bins{"ptBinsForMass", std::vector{o2::analysis::hf_cuts_dplus_topikpi::pTBins_v}, "pT bin limits for candidate mass plots"}; + + void init(o2::framework::InitContext&) + { + registry.add("hMass", "Dplus,Dminus candidates;inv. mass (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{120, 1.5848, 2.1848}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hMassDplus", "Dplus,Dminus candidates;inv. mass Dplus only (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{120, 1.5848, 2.1848}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hMassDminus", "Dplus,Dminus candidates;inv. mass Dminus only (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{120, 1.5848, 2.1848}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + } + + Filter filterSelectCandidates = (aod::hf_selcandidate_dplus::isSelDplusToPiKPi >= selectionFlagDplus); + + void process(aod::Collision const& collision, soa::Filtered> const& candidates, aod::BigTracks const& tracks) + { + for (auto& candidate1 : candidates) { + if (cutYCandMax >= 0. && std::abs(YDPlus(candidate1)) > cutYCandMax) { + continue; + } + if (cutPtCandMin >= 0. && candidate1.pt() < cutPtCandMin) { + continue; + } + //check decay channel flag for candidate1 + if (!(candidate1.hfflag() & 1 << DecayType::DPlusToPiKPi)) { //probably dummy since already selected? not sure... + continue; + } + int outerParticleSign = 1; //Dplus + auto outerSecondTrack = candidate1.index1_as(); + if (outerSecondTrack.sign() == 1) { + outerParticleSign = -1; //Dminus (second daughter track is positive) + } + //fill invariant mass plots and generic info from all Dplus/Dminus candidates + if (outerParticleSign == 1) { + registry.fill(HIST("hMass"), InvMassDPlus(candidate1), candidate1.pt()); + registry.fill(HIST("hMassDplus"), InvMassDPlus(candidate1), candidate1.pt()); + } else { + registry.fill(HIST("hMass"), InvMassDPlus(candidate1), candidate1.pt()); + registry.fill(HIST("hMassDminus"), InvMassDPlus(candidate1), candidate1.pt()); + } + registry.fill(HIST("hPtCand"), candidate1.pt()); + registry.fill(HIST("hPtProng0"), candidate1.ptProng0()); + registry.fill(HIST("hPtProng1"), candidate1.ptProng1()); + registry.fill(HIST("hPtProng2"), candidate1.ptProng2()); + registry.fill(HIST("hEta"), candidate1.eta()); + registry.fill(HIST("hPhi"), candidate1.phi()); + registry.fill(HIST("hY"), YDPlus(candidate1)); + registry.fill(HIST("hSelectionStatus"), candidate1.isSelDplusToPiKPi()); + + //D-Dbar correlation dedicated section + //For like-sign, first loop on both Dplus and Dminus. First candidate is for sure a Dplus/Dminus (checked before, so don't re-check anything on it) + for (auto& candidate2 : candidates) { + //check decay channel flag for candidate2 + if (!(candidate2.hfflag() & 1 << DecayType::DPlusToPiKPi)) { + continue; + } + //check if inner particle is Dplus or Dminus + int innerParticleSign = 1; //Dplus + auto innerSecondTrack = candidate2.index1_as(); + if (innerSecondTrack.sign() == 1) { + innerParticleSign = -1; //Dminus (second daughter track is positive) + } + //for the associated, has to have smaller pT, and has to be Dplus if outer is Dplu, or Dminus if outer is Dminus + if (innerParticleSign == outerParticleSign) { + if (cutYCandMax >= 0. && std::abs(YDPlus(candidate2)) > cutYCandMax) { + continue; + } + if (cutPtCandMin >= 0. && candidate2.pt() < cutPtCandMin) { + continue; + } + //Excluding self-correlations + if (candidate1.mRowIndex == candidate2.mRowIndex) { + continue; + } + entryDplusDminusPair(getDeltaPhi(candidate2.phi(), candidate1.phi()), + candidate2.eta() - candidate1.eta(), + candidate1.pt(), + candidate2.pt()); + entryDplusDminusRecoInfo(InvMassDPlus(candidate1), + InvMassDPlus(candidate2), + 0); + } + } // end inner loop + } //end outer loop + } +}; + +/// Dplus-Dminus correlation pair builder - LIKE SIGN - for MC reco analysis (data-like but matching to true Dplus and Dminus) +/// NOTE: At the moment, both dPhi-symmetrical correlation pairs (part1-part2 and part2-part1) are filled, since we bin in pT and selecting as trigger the largest pT particle would bias the distributions w.r.t. the ULS case. +struct HfCorrelatorDplusDminusMcRecLs { + Produces entryDplusDminusPair; + Produces entryDplusDminusRecoInfo; + + HistogramRegistry registry{ + "registry", + //NOTE: use hMassD0 for trigger normalisation (S*0.955), and hMass2DCorrelationPairs (in final task) for 2D-sideband-subtraction purposes + {{"hPtCandMCRec", "Dplus,Dminus candidates - MC reco;candidate #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, + {"hPtProng0MCRec", "Dplus,Dminus candidates - MC reco;prong 0 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, + {"hPtProng1MCRec", "Dplus,Dminus candidates - MC reco;prong 1 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, + {"hPtProng2MCRec", "Dplus,Dminus candidates - MC reco;prong 2 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, + {"hSelectionStatusMCRec", "Dplus,Dminus candidates - MC reco;selection status;entries", {HistType::kTH1F, {{2, -0.5, 1.5}}}}, + {"hEtaMCRec", "Dplus,Dminus candidates - MC reco;candidate #it{#eta};entries", {HistType::kTH1F, {{100, -5., 5.}}}}, + {"hPhiMCRec", "Dplus,Dminus candidates - MC reco;candidate #it{#varphi};entries", {HistType::kTH1F, {{32, 0., 2. * o2::constants::math::PI}}}}, + {"hYMCRec", "Dplus,Dminus candidates - MC reco;candidate #it{#y};entries", {HistType::kTH1F, {{100, -5., 5.}}}}}}; + + Configurable selectionFlagDplus{"selectionFlagDplus", 1, "Selection Flag for Dplus,Dminus"}; + Configurable cutYCandMax{"cutYCandMax", -1., "max. cand. rapidity"}; + Configurable cutPtCandMin{"cutPtCandMin", -1., "min. cand. pT"}; + Configurable> bins{"ptBinsForMass", std::vector{o2::analysis::hf_cuts_dplus_topikpi::pTBins_v}, "pT bin limits for candidate mass plots"}; + + void init(o2::framework::InitContext&) + { + registry.add("hMassDplusMCRec", "Dplus,Dminus candidates - MC reco;inv. mass (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{120, 1.5848, 2.1848}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hMassDminusMCRec", "Dplus,Dminus candidates - MC reco;inv. mass D0 only (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{120, 1.5848, 2.1848}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + } + + Filter filterSelectCandidates = (aod::hf_selcandidate_dplus::isSelDplusToPiKPi >= selectionFlagDplus); + + void process(aod::Collision const& collision, soa::Filtered> const& candidates, aod::BigTracks const& tracks) + { + //MC reco level + for (auto& candidate1 : candidates) { + //check decay channel flag for candidate1 + if (!(candidate1.hfflag() & 1 << DecayType::DPlusToPiKPi)) { + continue; + } + if (cutYCandMax >= 0. && std::abs(YDPlus(candidate1)) > cutYCandMax) { + continue; + } + if (cutPtCandMin >= 0. && candidate1.pt() < cutPtCandMin) { + continue; + } + int outerParticleSign = 1; //Dplus + auto outerSecondTrack = candidate1.index1_as(); + if (outerSecondTrack.sign() == 1) { + outerParticleSign = -1; //Dminus (second daughter track is positive) + } + if (std::abs(candidate1.flagMCMatchRec()) != 1 << DecayType::DPlusToPiKPi) { //keep only Dplus/Dminus matched candidates + continue; + } + //fill invariant mass plots and generic info from all Dplus/Dminus candidates + if (outerParticleSign == 1) { //matched as Dplus + registry.fill(HIST("hMassDplusMCRec"), InvMassDPlus(candidate1), candidate1.pt()); + } else { //matched as Dminus + registry.fill(HIST("hMassDminusMCRec"), InvMassDPlus(candidate1), candidate1.pt()); + } + registry.fill(HIST("hPtCandMCRec"), candidate1.pt()); + registry.fill(HIST("hPtProng0MCRec"), candidate1.ptProng0()); + registry.fill(HIST("hPtProng1MCRec"), candidate1.ptProng1()); + registry.fill(HIST("hPtProng2MCRec"), candidate1.ptProng2()); + registry.fill(HIST("hEtaMCRec"), candidate1.eta()); + registry.fill(HIST("hPhiMCRec"), candidate1.phi()); + registry.fill(HIST("hYMCRec"), YDPlus(candidate1)); + registry.fill(HIST("hSelectionStatusMCRec"), candidate1.isSelDplusToPiKPi()); + + //D-Dbar correlation dedicated section + //For like-sign, first loop on both Dplus and Dminus. First candidate is for sure a Dplus/Dminus (looping on filtered) and was already matched, so don't re-check anything on it) + for (auto& candidate2 : candidates) { + //check decay channel flag for candidate2 + if (!(candidate2.hfflag() & 1 << DecayType::DPlusToPiKPi)) { + continue; + } + int innerParticleSign = 1; //Dplus + auto innerSecondTrack = candidate2.index1_as(); + if (innerSecondTrack.sign() == 1) { + innerParticleSign = -1; //Dminus (second daughter track is positive) + } + if (std::abs(candidate2.flagMCMatchRec()) != 1 << DecayType::DPlusToPiKPi) { //reject fake candidates + continue; + } + if (outerParticleSign == innerParticleSign) { //LS pair (of Dplus or Dminus) + pt2= 0. && std::abs(YDPlus(candidate2)) > cutYCandMax) { + continue; + } + if (cutPtCandMin >= 0. && candidate2.pt() < cutPtCandMin) { + continue; + } + //Excluding self-correlations (in principle not possible due to the '<' condition, but could rounding break it?) + if (candidate1.mRowIndex == candidate2.mRowIndex) { + continue; + } + entryDplusDminusPair(getDeltaPhi(candidate2.phi(), candidate1.phi()), + candidate2.eta() - candidate1.eta(), + candidate1.pt(), + candidate2.pt()); + entryDplusDminusRecoInfo(InvMassDPlus(candidate1), + InvMassDPlus(candidate2), + 0); //for LS studies we set a dummy 0 for pairSignalStatus (there are no more the usual 4 possible combinations) + + } //end inner if (MC match) + } // end inner loop + } //end outer loop + } +}; + +/// Dplus-Dminus correlation pair builder - for MC gen-level analysis, like sign particles +/// NOTE: At the moment, both dPhi-symmetrical correlation pairs (part1-part2 and part2-part1) are filled, since we bin in pT and selecting as trigger the largest pT particle would bias the distributions w.r.t. the ULS case. +struct HfCorrelatorDplusDminusMcGenLs { + + Produces entryDplusDminusPair; + + HistogramRegistry registry{ + "registry", + {{"hMCEvtCount", "Event counter - MC gen;;entries", {HistType::kTH1F, {{1, -0.5, 0.5}}}}, + {"hPtCandMCGen", "Dplus,Dminus particles - MC gen;particle #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, + {"hEtaMCGen", "Dplus,Dminus particles - MC gen;particle #it{#eta};entries", {HistType::kTH1F, {{100, -5., 5.}}}}, + {"hPhiMCGen", "Dplus,Dminus particles - MC gen;particle #it{#varphi};entries", {HistType::kTH1F, {{32, 0., 2. * o2::constants::math::PI}}}}, + {"hYMCGen", "Dplus,Dminus candidates - MC gen;candidate #it{#y};entries", {HistType::kTH1F, {{100, -5., 5.}}}}, + {"hcountDplusDminusPerEvent", "Dplus,Dminus particles - MC gen;Number per event;entries", {HistType::kTH1F, {{20, 0., 20.}}}}}}; + + Configurable cutYCandMax{"cutYCandMax", -1., "max. cand. rapidity"}; + Configurable cutPtCandMin{"cutPtCandMin", -1., "min. cand. pT"}; + Configurable> bins{"ptBinsForMass", std::vector{o2::analysis::hf_cuts_dplus_topikpi::pTBins_v}, "pT bin limits for trigger counters"}; + + void init(o2::framework::InitContext&) + { + registry.add("hcountDplustriggersMCGen", "Dplus trigger particles - MC gen;;N of trigger D0", {HistType::kTH2F, {{1, -0.5, 0.5}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + } + + void process(aod::McCollision const& mccollision, soa::Join const& particlesMC) + { + int counterDplusDminus = 0; + registry.fill(HIST("hMCEvtCount"), 0); + //MC gen level + for (auto& particle1 : particlesMC) { + //check if the particle is Dplus or Dminus (both can be trigger) - NOTE: decay channel is not probed! + if (std::abs(particle1.pdgCode()) != pdg::Code::kDPlus) { + continue; + } + double yD = RecoDecay::Y(array{particle1.px(), particle1.py(), particle1.pz()}, RecoDecay::getMassPDG(particle1.pdgCode())); + if (cutYCandMax >= 0. && std::abs(yD) > cutYCandMax) { + continue; + } + if (cutPtCandMin >= 0. && particle1.pt() < cutPtCandMin) { + continue; + } + + registry.fill(HIST("hPtCandMCGen"), particle1.pt()); + registry.fill(HIST("hEtaMCGen"), particle1.eta()); + registry.fill(HIST("hPhiMCGen"), particle1.phi()); + registry.fill(HIST("hYMCGen"), yD); + counterDplusDminus++; + + //D-Dbar correlation dedicated section + //if it's Dplus/Dminus, search for LS pair companions and evaluate correlations. + registry.fill(HIST("hcountDplustriggersMCGen"), 0, particle1.pt()); //to count trigger Dplus (normalisation) + for (auto& particle2 : particlesMC) { + if (std::abs(particle2.pdgCode()) != pdg::Code::kDPlus) { //check that associated is a Dplus/Dminus (both are fine) + continue; + } + if (cutYCandMax >= 0. && std::abs(RecoDecay::Y(array{particle2.px(), particle2.py(), particle2.pz()}, RecoDecay::getMassPDG(particle2.pdgCode()))) > cutYCandMax) { + continue; + } + if (cutPtCandMin >= 0. && particle2.pt() < cutPtCandMin) { + continue; + } + if (particle2.pdgCode() == particle1.pdgCode()) { //like-sign condition (both 411 or both -411) and pT_Trig>pT_assoc + //Excluding self-correlations (in principle not possible due to the '<' condition, but could rounding break it?) + if (particle1.mRowIndex == particle2.mRowIndex) { + continue; + } + entryDplusDminusPair(getDeltaPhi(particle2.phi(), particle1.phi()), + particle2.eta() - particle1.eta(), + particle1.pt(), + particle2.pt()); + } + } // end inner loop + } //end outer loop + registry.fill(HIST("hcountDplusDminusPerEvent"), counterDplusDminus); + } +}; + +/// c-cbar correlator table builder - for MC gen-level analysis +struct HfCorrelatorCCbarMcGen { + + Produces entryccbarPair; + + HistogramRegistry registry{ + "registry", + {{"hMCEvtCount", "Event counter - MC gen;;entries", {HistType::kTH1F, {{1, -0.5, 0.5}}}}, + {"hPtCandMCGen", "c,cbar particles - MC gen;particle #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, + {"hEtaMCGen", "c,cbar particles - MC gen;particle #it{#eta};entries", {HistType::kTH1F, {{100, -5., 5.}}}}, + {"hYMCGen", "c,cbar candidates - MC gen;candidate #it{#y};entries", {HistType::kTH1F, {{100, -5., 5.}}}}, + {"hPhiMCGen", "c,cbar particles - MC gen;particle #it{#varphi};entries", {HistType::kTH1F, {{32, 0., 2. * o2::constants::math::PI}}}}, + {"hcountCCbarPerEvent", "c,cbar particles - MC gen;Number per event;entries", {HistType::kTH1F, {{20, 0., 20.}}}}, + {"hcountCCbarPerEventPreEtaCut", "c,cbar particles - MC gen;Number per event pre #eta cut;entries", {HistType::kTH1F, {{20, 0., 20.}}}}}}; + + Configurable cutYCandMax{"cutYCandMax", -1., "max. cand. rapidity"}; + Configurable cutPtCandMin{"cutPtCandMin", -1., "min. cand. pT"}; + Configurable> bins{"ptBinsForMass", std::vector{o2::analysis::hf_cuts_d0_topik::pTBins_v}, "pT bin limits for trigger counters"}; + + void init(o2::framework::InitContext&) + { + registry.add("hcountCtriggersMCGen", "c trigger particles - MC gen;;N of trigger c quark", {HistType::kTH2F, {{1, -0.5, 0.5}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + } + + void process(aod::McCollision const& mccollision, soa::Join const& particlesMC) + { + registry.fill(HIST("hMCEvtCount"), 0); + int counterccbar = 0, counterccbarPreEtasel = 0; + + //loop over particles at MC gen level + for (auto& particle1 : particlesMC) { + if (std::abs(particle1.pdgCode()) != PDG_t::kCharm) { //search c or cbar particles + continue; + } + int partMothPDG = particlesMC.iteratorAt(particle1.mother0()).pdgCode(); + //check whether mothers of quark c/cbar are still '4'/'-4' particles - in that case the c/cbar quark comes from its own fragmentation, skip it + if (partMothPDG == particle1.pdgCode()) { + continue; + } + counterccbarPreEtasel++; //count c or cbar (before kinematic selection) + double yC = RecoDecay::Y(array{particle1.px(), particle1.py(), particle1.pz()}, RecoDecay::getMassPDG(particle1.pdgCode())); + if (cutYCandMax >= 0. && std::abs(yC) > cutYCandMax) { + continue; + } + if (cutPtCandMin >= 0. && particle1.pt() < cutPtCandMin) { + continue; + } + registry.fill(HIST("hPtCandMCGen"), particle1.pt()); + registry.fill(HIST("hEtaMCGen"), particle1.eta()); + registry.fill(HIST("hPhiMCGen"), particle1.phi()); + registry.fill(HIST("hYMCGen"), yC); + counterccbar++; //count if c or cbar don't come from themselves during fragmentation (after kinematic selection) + + //c-cbar correlation dedicated section + //if it's c, search for cbar and evaluate correlations. + if (particle1.pdgCode() != PDG_t::kCharm) { + continue; + } + registry.fill(HIST("hcountCtriggersMCGen"), 0, particle1.pt()); //to count trigger c quark (for normalisation) + + for (auto& particle2 : particlesMC) { + if (particle2.pdgCode() != PDG_t::kCharmBar) { + continue; + } + if (cutYCandMax >= 0. && std::abs(RecoDecay::Y(array{particle2.px(), particle2.py(), particle2.pz()}, RecoDecay::getMassPDG(particle2.pdgCode()))) > cutYCandMax) { + continue; + } + if (cutPtCandMin >= 0. && particle2.pt() < cutPtCandMin) { + continue; + } + //check whether mothers of quark cbar (from associated loop) are still '-4' particles - in that case the cbar quark comes from its own fragmentation, skip it + if (particlesMC.iteratorAt(particle2.mother0()).pdgCode() == PDG_t::kCharmBar) { + continue; + } + entryccbarPair(getDeltaPhi(particle2.phi(), particle1.phi()), + particle2.eta() - particle1.eta(), + particle1.pt(), + particle2.pt()); + } // end inner loop + } //end outer loop + registry.fill(HIST("hcountCCbarPerEvent"), counterccbar); + registry.fill(HIST("hcountCCbarPerEventPreEtaCut"), counterccbarPreEtasel); + } +}; + +/// c-cbar correlator table builder - for MC gen-level analysis - Like Sign +struct HfCorrelatorCCbarMcGenLs { + + Produces entryccbarPair; + + HistogramRegistry registry{ + "registry", + {{"hMCEvtCount", "Event counter - MC gen;;entries", {HistType::kTH1F, {{1, -0.5, 0.5}}}}, + {"hPtCandMCGen", "c,cbar particles - MC gen;particle #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, + {"hEtaMCGen", "c,cbar particles - MC gen;particle #it{#eta};entries", {HistType::kTH1F, {{100, -5., 5.}}}}, + {"hYMCGen", "c,cbar candidates - MC gen;candidate #it{#y};entries", {HistType::kTH1F, {{100, -5., 5.}}}}, + {"hPhiMCGen", "c,cbar particles - MC gen;particle #it{#varphi};entries", {HistType::kTH1F, {{32, 0., 2. * o2::constants::math::PI}}}}, + {"hcountCCbarPerEvent", "c,cbar particles - MC gen;Number per event;entries", {HistType::kTH1F, {{20, 0., 20.}}}}, + {"hcountCCbarPerEventPreEtaCut", "c,cbar particles - MC gen;Number per event pre #eta cut;entries", {HistType::kTH1F, {{20, 0., 20.}}}}}}; + + Configurable cutYCandMax{"cutYCandMax", -1., "max. cand. rapidity"}; + Configurable cutPtCandMin{"cutPtCandMin", -1., "min. cand. pT"}; + Configurable> bins{"ptBinsForMass", std::vector{o2::analysis::hf_cuts_d0_topik::pTBins_v}, "pT bin limits for trigger counters"}; + + void init(o2::framework::InitContext&) + { + registry.add("hcountCtriggersMCGen", "c trigger particles - MC gen;;N of trigger c quark", {HistType::kTH2F, {{1, -0.5, 0.5}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + } + + void process(aod::McCollision const& mccollision, soa::Join const& particlesMC) + { + registry.fill(HIST("hMCEvtCount"), 0); + int counterccbar = 0, counterccbarPreEtasel = 0; + + //loop over particles at MC gen level + for (auto& particle1 : particlesMC) { + if (std::abs(particle1.pdgCode()) != PDG_t::kCharm) { //search c or cbar particles + continue; + } + int partMothPDG = particlesMC.iteratorAt(particle1.mother0()).pdgCode(); + //check whether mothers of quark c/cbar are still '4'/'-4' particles - in that case the c/cbar quark comes from its own fragmentation, skip it + if (partMothPDG == particle1.pdgCode()) { + continue; + } + counterccbarPreEtasel++; //count c or cbar (before kinematic selection) + double yC = RecoDecay::Y(array{particle1.px(), particle1.py(), particle1.pz()}, RecoDecay::getMassPDG(particle1.pdgCode())); + if (cutYCandMax >= 0. && std::abs(yC) > cutYCandMax) { + continue; + } + if (cutPtCandMin >= 0. && particle1.pt() < cutPtCandMin) { + continue; + } + registry.fill(HIST("hPtCandMCGen"), particle1.pt()); + registry.fill(HIST("hEtaMCGen"), particle1.eta()); + registry.fill(HIST("hPhiMCGen"), particle1.phi()); + registry.fill(HIST("hYMCGen"), yC); + counterccbar++; //count if c or cbar don't come from themselves during fragmentation (after kinematic selection) + + //c-cbar correlation dedicated section + registry.fill(HIST("hcountCtriggersMCGen"), 0, particle1.pt()); //to count trigger c quark (for normalisation) + + for (auto& particle2 : particlesMC) { + if (std::abs(particle2.pdgCode()) != PDG_t::kCharm) { //search c or cbar for associated particles + continue; + } + if (cutYCandMax >= 0. && std::abs(RecoDecay::Y(array{particle2.px(), particle2.py(), particle2.pz()}, RecoDecay::getMassPDG(particle2.pdgCode()))) > cutYCandMax) { + continue; + } + if (cutPtCandMin >= 0. && particle2.pt() < cutPtCandMin) { + continue; + } + if (particle2.pdgCode() == particle1.pdgCode()) { + //check whether mothers of quark cbar (from associated loop) are still '-4' particles - in that case the cbar quark comes from its own fragmentation, skip it + if (particlesMC.iteratorAt(particle2.mother0()).pdgCode() == particle2.pdgCode()) { + continue; + } + //Excluding self-correlations (in principle not possible due to the '<' condition, but could rounding break it?) + if (particle1.mRowIndex == particle2.mRowIndex) { + continue; + } + entryccbarPair(getDeltaPhi(particle2.phi(), particle1.phi()), + particle2.eta() - particle1.eta(), + particle1.pt(), + particle2.pt()); + } // end outer if (check PDG associate) + } // end inner loop + } //end outer loop + registry.fill(HIST("hcountCCbarPerEvent"), counterccbar); + registry.fill(HIST("hcountCCbarPerEventPreEtaCut"), counterccbarPreEtasel); + } +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + WorkflowSpec workflow{}; + const bool doMCccbar = cfgc.options().get("doMCccbar"); + const bool doMCGen = cfgc.options().get("doMCGen"); + const bool doMCRec = cfgc.options().get("doMCRec"); + const bool doLikeSign = cfgc.options().get("doLikeSign"); + if (!doLikeSign) { //unlike-sign analyses + if (doMCGen) { //MC-Gen analysis + workflow.push_back(adaptAnalysisTask(cfgc)); + } else if (doMCRec) { //MC-Reco analysis + workflow.push_back(adaptAnalysisTask(cfgc)); + } else if (doMCccbar) { //MC-Reco analysis + workflow.push_back(adaptAnalysisTask(cfgc)); + } else { //data analysis + workflow.push_back(adaptAnalysisTask(cfgc)); + } + } else { //like-sign analyses + if (doMCGen) { //MC-Gen analysis + workflow.push_back(adaptAnalysisTask(cfgc)); + } else if (doMCRec) { //MC-Reco analysis + workflow.push_back(adaptAnalysisTask(cfgc)); + } else if (doMCccbar) { //MC-Reco analysis + workflow.push_back(adaptAnalysisTask(cfgc)); + } else { //data analysis + workflow.push_back(adaptAnalysisTask(cfgc)); + } + } + + return workflow; +} diff --git a/Analysis/Tasks/PWGHF/taskCorrelationDDbar.cxx b/Analysis/Tasks/PWGHF/taskCorrelationDDbar.cxx new file mode 100644 index 0000000000000..b7f0b9bf72e21 --- /dev/null +++ b/Analysis/Tasks/PWGHF/taskCorrelationDDbar.cxx @@ -0,0 +1,391 @@ +// Copyright CERN and copyright holders of ALICE O2. This software is +// distributed under the terms of the GNU General Public License v3 (GPL +// Version 3), copied verbatim in the file "COPYING". +// +// See http://alice-o2.web.cern.ch/license for full licensing information. +// +// 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 taskCorrelationDDbar.cxx +/// \brief D-Dbar analysis task - data-like, MC-reco and MC-kine analyses. For ULS and LS pairs +/// +/// \author Fabio Colamaria , INFN Bari + +#include "Framework/AnalysisTask.h" +#include "Framework/HistogramRegistry.h" +#include "AnalysisDataModel/HFSecondaryVertex.h" +#include "AnalysisDataModel/HFCandidateSelectionTables.h" + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; +using namespace o2::aod::hf_cand_prong2; +using namespace o2::aod::hf_correlation_ddbar; +using namespace o2::analysis::hf_cuts_d0_topik; +using namespace o2::constants::math; + +namespace o2::aod +{ +using DDbarPairFull = soa::Join; +} // namespace o2::aod + +void customize(std::vector& workflowOptions) +{ + ConfigParamSpec optionDoMCGen{"doMCGen", VariantType::Bool, false, {"Run MC-Gen dedicated tasks."}}; + ConfigParamSpec optionDoMCRec{"doMCRec", VariantType::Bool, false, {"Run MC-Rec dedicated tasks."}}; + workflowOptions.push_back(optionDoMCGen); + workflowOptions.push_back(optionDoMCRec); +} + +#include "Framework/runDataProcessing.h" + +/// +/// Returns deltaPhi value in range [-pi/2., 3.*pi/2], typically used for correlation studies +/// +double getDeltaPhi(double phiD, double phiDbar) +{ + return RecoDecay::constrainAngle(phiDbar - phiD, -o2::constants::math::PI / 2.); +} + +/// +/// Returns phi of candidate/particle evaluated from x and y components of segment connecting primary and secondary vertices +/// +double evaluatePhiByVertex(double xVertex1, double xVertex2, double yVertex1, double yVertex2) +{ + return RecoDecay::Phi(xVertex2 - xVertex1, yVertex2 - yVertex1); +} +/* +/// definition of axes for THnSparse (dPhi,dEta,pTD,pTDbar) - note: last two axis are dummy, will be replaced later +const int nBinsSparse[4] = {32,120,10,10}; +const double binMinSparse[4] = {-o2::constants::math::PI / 2.,-6.,0.,0.; //is the minimum for all the bins +const double binMaxSparse[4] = {3. * o2::constants::math::PI / 2.,6.,10.,10.}; //is the maximum for all the bins +*/ + +// string definitions, used for histogram axis labels +const TString stringPtD = "#it{p}_{T}^{D} (GeV/#it{c});"; +const TString stringPtDbar = "#it{p}_{T}^{Dbar} (GeV/#it{c});"; +const TString stringDeltaPt = "#it{p}_{T}^{Dbar}-#it{p}_{T}^{D} (GeV/#it{c});"; +const TString stringDeltaPtMaxMin = "#it{p}_{T}^{max}-#it{p}_{T}^{min} (GeV/#it{c});"; +const TString stringDeltaEta = "#it{#eta}^{Dbar}-#it{#eta}^{D};"; +const TString stringDeltaPhi = "#it{#varphi}^{Dbar}-#it{#varphi}^{D} (rad);"; +const TString stringDDbar = "D,Dbar candidates "; +const TString stringSignal = "signal region;"; +const TString stringSideband = "sidebands;"; +const TString stringMCParticles = "MC gen - D,Dbar particles;"; +const TString stringMCReco = "MC reco - D,Dbar candidates "; + +/// D-Dbar 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 +struct HfTaskCorrelationDDbar { + + HistogramRegistry registry{ + "registry", + //NOTE: use hMassD0 (from correlator task) for normalisation, and hMass2DCorrelationPairs for 2D-sideband-subtraction purposes + {{"hMass2DCorrelationPairs", stringDDbar + "2D;inv. mass D (GeV/#it{c}^{2});inv. mass Dbar (GeV/#it{c}^{2});" + stringPtD + stringPtDbar + "entries", {HistType::kTHnSparseD, {{200, 1.6, 2.1}, {200, 1.6, 2.1}, {10, 0., 10.}, {10, 0., 10.}}}}, //note: axes 3 and 4 (the pT) are updated in the init() + {"hDeltaEtaPtIntSignalRegion", stringDDbar + stringSignal + stringDeltaEta + "entries", {HistType::kTH1F, {{200, -10., 10.}}}}, + {"hDeltaPhiPtIntSignalRegion", stringDDbar + stringSignal + stringDeltaPhi + "entries", {HistType::kTH1F, {{32, -o2::constants::math::PI / 2., 3. * o2::constants::math::PI / 2.}}}}, + {"hCorrel2DPtIntSignalRegion", stringDDbar + stringSignal + stringDeltaPhi + stringDeltaEta + "entries", {HistType::kTH2F, {{32, -o2::constants::math::PI / 2., 3. * o2::constants::math::PI / 2.}, {200, -10., 10.}}}}, + {"hCorrel2DVsPtSignalRegion", stringDDbar + stringSignal + stringDeltaPhi + stringDeltaEta + stringPtD + stringPtDbar + "entries", {HistType::kTHnSparseD, {{32, -o2::constants::math::PI / 2., 3. * o2::constants::math::PI / 2.}, {120, -6., 6.}, {10, 0., 10.}, {10, 0., 10.}}}}, //note: axes 3 and 4 (the pT) are updated in the init() + {"hDeltaPtDDbarSignalRegion", stringDDbar + stringSignal + stringDeltaPt + "entries", {HistType::kTH1F, {{144, -36., 36.}}}}, + {"hDeltaPtMaxMinSignalRegion", stringDDbar + stringSignal + stringDeltaPtMaxMin + "entries", {HistType::kTH1F, {{72, 0., 36.}}}}, + {"hDeltaEtaPtIntSidebands", stringDDbar + stringSideband + stringDeltaEta + "entries", {HistType::kTH1F, {{200, -10., 10.}}}}, + {"hDeltaPhiPtIntSidebands", stringDDbar + stringSideband + stringDeltaPhi + "entries", {HistType::kTH1F, {{32, -o2::constants::math::PI / 2., 3. * o2::constants::math::PI / 2.}}}}, + {"hCorrel2DPtIntSidebands", stringDDbar + stringSideband + stringDeltaPhi + stringDeltaEta + "entries", {HistType::kTH2F, {{32, -o2::constants::math::PI / 2., 3. * o2::constants::math::PI / 2.}, {200, -10., 10.}}}}, + {"hCorrel2DVsPtSidebands", stringDDbar + stringSideband + stringDeltaPhi + stringDeltaEta + stringPtD + stringPtDbar + "entries", {HistType::kTHnSparseD, {{32, -o2::constants::math::PI / 2., 3. * o2::constants::math::PI / 2.}, {120, -6., 6.}, {10, 0., 10.}, {10, 0., 10.}}}}, //note: axes 3 and 4 (the pT) are updated in the init() + {"hDeltaPtDDbarSidebands", stringDDbar + stringSideband + stringDeltaPt + "entries", {HistType::kTH1F, {{144, -36., 36.}}}}, + {"hDeltaPtMaxMinSidebands", stringDDbar + stringSideband + stringDeltaPtMaxMin + "entries", {HistType::kTH1F, {{72, 0., 36.}}}}}}; + + //pT ranges for correlation plots: the default values are those embedded in hf_cuts_d0_topik (i.e. the mass pT bins), but can be redefined via json files + Configurable> binsCorrelations{"ptBinsForCorrelations", std::vector{o2::analysis::hf_cuts_d0_topik::pTBins_v}, "pT bin limits for correlation plots"}; + //signal and sideband region edges, to be defined via json file (initialised to empty) + Configurable> signalRegionInner{"signalRegionInner", std::vector(), "Inner values of signal region vs pT"}; + Configurable> signalRegionOuter{"signalRegionOuter", std::vector(), "Outer values of signal region vs pT"}; + Configurable> sidebandLeftInner{"sidebandLeftInner", std::vector(), "Inner values of left sideband vs pT"}; + Configurable> sidebandLeftOuter{"sidebandLeftOuter", std::vector(), "Outer values of left sideband vs pT"}; + Configurable> sidebandRightInner{"sidebandRightInner", std::vector(), "Inner values of right sideband vs pT"}; + Configurable> sidebandRightOuter{"sidebandRightOuter", std::vector(), "Outer values of right sideband vs pT"}; + + void init(o2::framework::InitContext&) + { + // redefinition of pT axes for THnSparse holding correlation entries + int nBinspTaxis = binsCorrelations->size() - 1; + const double* valuespTaxis = binsCorrelations->data(); + + for (int i = 2; i <= 3; i++) { + registry.get(HIST("hMass2DCorrelationPairs"))->GetAxis(i)->Set(nBinspTaxis, valuespTaxis); + registry.get(HIST("hCorrel2DVsPtSignalRegion"))->GetAxis(i)->Set(nBinspTaxis, valuespTaxis); + registry.get(HIST("hCorrel2DVsPtSidebands"))->GetAxis(i)->Set(nBinspTaxis, valuespTaxis); + } + } + + void process(aod::DDbarPairFull const& pairEntries) + { + for (auto& pairEntry : pairEntries) { + //define variables for widely used quantities + double deltaPhi = pairEntry.deltaPhi(); + double deltaEta = pairEntry.deltaEta(); + double ptD = pairEntry.ptD(); + double ptDbar = pairEntry.ptDbar(); + double massD = pairEntry.mD(); + double massDbar = pairEntry.mDbar(); + + //fill 2D invariant mass plots + registry.fill(HIST("hMass2DCorrelationPairs"), massD, massDbar, ptD, ptDbar); + + //reject entries outside pT ranges of interest + int pTBinD = o2::analysis::findBin(binsCorrelations, ptD); + int pTBinDbar = o2::analysis::findBin(binsCorrelations, ptDbar); + if (pTBinD == -1 || pTBinDbar == -1) { //at least one particle outside accepted pT range + continue; + } + + //check if correlation entry belongs to signal region, sidebands or is outside both, and fill correlation plots + if (massD > signalRegionInner->at(pTBinD) && massD < signalRegionOuter->at(pTBinD) && massDbar > signalRegionInner->at(pTBinDbar) && massDbar < signalRegionOuter->at(pTBinDbar)) { + //in signal region + registry.fill(HIST("hCorrel2DVsPtSignalRegion"), deltaPhi, deltaEta, ptD, ptDbar); + registry.fill(HIST("hCorrel2DPtIntSignalRegion"), deltaPhi, deltaEta); + registry.fill(HIST("hDeltaEtaPtIntSignalRegion"), deltaEta); + registry.fill(HIST("hDeltaPhiPtIntSignalRegion"), deltaPhi); + registry.fill(HIST("hDeltaPtDDbarSignalRegion"), ptDbar - ptD); + registry.fill(HIST("hDeltaPtMaxMinSignalRegion"), std::abs(ptDbar - ptD)); + } + + if ((massD > sidebandLeftInner->at(pTBinD) && massD < sidebandLeftOuter->at(pTBinD) && massDbar > sidebandLeftInner->at(pTBinDbar) && massDbar < sidebandRightOuter->at(pTBinDbar)) || + (massD > sidebandRightInner->at(pTBinD) && massD < sidebandRightOuter->at(pTBinD) && massDbar > sidebandLeftInner->at(pTBinDbar) && massDbar < sidebandRightOuter->at(pTBinDbar)) || + (massD > sidebandLeftInner->at(pTBinD) && massD < sidebandRightOuter->at(pTBinD) && massDbar > sidebandLeftInner->at(pTBinDbar) && massDbar < sidebandLeftOuter->at(pTBinDbar)) || + (massD > sidebandLeftInner->at(pTBinD) && massD < sidebandRightOuter->at(pTBinD) && massDbar > sidebandRightInner->at(pTBinDbar) && massDbar < sidebandRightOuter->at(pTBinDbar))) { + //in sideband region + registry.fill(HIST("hCorrel2DVsPtSidebands"), deltaPhi, deltaEta, ptD, ptDbar); + registry.fill(HIST("hCorrel2DPtIntSidebands"), deltaPhi, deltaEta); + registry.fill(HIST("hDeltaEtaPtIntSidebands"), deltaEta); + registry.fill(HIST("hDeltaPhiPtIntSidebands"), deltaPhi); + registry.fill(HIST("hDeltaPtDDbarSidebands"), ptDbar - ptD); + registry.fill(HIST("hDeltaPtMaxMinSidebands"), std::abs(ptDbar - ptD)); + } + } //end loop + } +}; + +/// D-Dbar correlation pair filling task, from pair tables - for MC reco-level analysis (candidates matched to true signal only, but also bkg sources are studied) +/// Works on both USL and LS analyses pair tables +struct HfTaskCorrelationDDbarMcRec { + + HistogramRegistry registry{ + "registry", + //NOTE: use hMassD0 (from correlator task) for normalisation, and hMass2DCorrelationPairs for 2D-sideband-subtraction purposes + {{"hMass2DCorrelationPairsMCRecSigSig", stringDDbar + "2D SigSig - MC reco;inv. mass D (GeV/#it{c}^{2});inv. mass Dbar (GeV/#it{c}^{2});" + stringPtD + stringPtDbar + "entries", {HistType::kTHnSparseD, {{200, 1.6, 2.1}, {200, 1.6, 2.1}, {10, 0., 10.}, {10, 0., 10.}}}}, //note: axes 3 and 4 (the pT) are updated in the init() + {"hMass2DCorrelationPairsMCRecSigBkg", stringDDbar + "2D SigBkg - MC reco;inv. mass D (GeV/#it{c}^{2});inv. mass Dbar (GeV/#it{c}^{2});" + stringPtD + stringPtDbar + "entries", {HistType::kTHnSparseD, {{200, 1.6, 2.1}, {200, 1.6, 2.1}, {10, 0., 10.}, {10, 0., 10.}}}}, //note: axes 3 and 4 (the pT) are updated in the init() + {"hMass2DCorrelationPairsMCRecBkgSig", stringDDbar + "2D BkgSig - MC reco;inv. mass D (GeV/#it{c}^{2});inv. mass Dbar (GeV/#it{c}^{2});" + stringPtD + stringPtDbar + "entries", {HistType::kTHnSparseD, {{200, 1.6, 2.1}, {200, 1.6, 2.1}, {10, 0., 10.}, {10, 0., 10.}}}}, //note: axes 3 and 4 (the pT) are updated in the init() + {"hMass2DCorrelationPairsMCRecBkgBkg", stringDDbar + "2D BkgBkg - MC reco;inv. mass D (GeV/#it{c}^{2});inv. mass Dbar (GeV/#it{c}^{2});" + stringPtD + stringPtDbar + "entries", {HistType::kTHnSparseD, {{200, 1.6, 2.1}, {200, 1.6, 2.1}, {10, 0., 10.}, {10, 0., 10.}}}}, //note: axes 3 and 4 (the pT) are updated in the init() + {"hDeltaEtaPtIntSignalRegionMCRec", stringMCReco + stringSignal + stringDeltaEta + "entries", {HistType::kTH1F, {{200, -10., 10.}}}}, + {"hDeltaPhiPtIntSignalRegionMCRec", stringMCReco + stringSignal + stringDeltaPhi + "entries", {HistType::kTH1F, {{32, -o2::constants::math::PI / 2., 3. * o2::constants::math::PI / 2.}}}}, + {"hCorrel2DPtIntSignalRegionMCRec", stringMCReco + stringSignal + stringDeltaPhi + stringDeltaEta + "entries", {HistType::kTH2F, {{32, -o2::constants::math::PI / 2., 3. * o2::constants::math::PI / 2.}, {200, -10., 10.}}}}, + {"hDeltaPtDDbarSignalRegionMCRec", stringMCReco + stringSignal + stringDeltaPt + "entries", {HistType::kTH1F, {{144, -36., 36.}}}}, + {"hDeltaPtMaxMinSignalRegionMCRec", stringMCReco + stringSignal + stringDeltaPtMaxMin + "entries", {HistType::kTH1F, {{72, 0., 36.}}}}, + {"hCorrel2DVsPtSignalRegionMCRecSigSig", stringMCReco + "SigSig" + stringSignal + stringDeltaPhi + stringDeltaEta + stringPtD + stringPtDbar + "entries", {HistType::kTHnSparseD, {{32, -o2::constants::math::PI / 2., 3. * o2::constants::math::PI / 2.}, {120, -6., 6.}, {10, 0., 10.}, {10, 0., 10.}}}}, //note: axes 3 and 4 (the pT) are updated in the init() + {"hCorrel2DVsPtSignalRegionMCRecSigBkg", stringMCReco + "SigBkg" + stringSignal + +stringDeltaPhi + stringDeltaEta + stringPtD + stringPtDbar + "entries", {HistType::kTHnSparseD, {{32, -o2::constants::math::PI / 2., 3. * o2::constants::math::PI / 2.}, {120, -6., 6.}, {10, 0., 10.}, {10, 0., 10.}}}}, //note: axes 3 and 4 (the pT) are updated in the init() + {"hCorrel2DVsPtSignalRegionMCRecBkgSig", stringMCReco + "BkgSig" + stringSignal + +stringDeltaPhi + stringDeltaEta + stringPtD + stringPtDbar + "entries", {HistType::kTHnSparseD, {{32, -o2::constants::math::PI / 2., 3. * o2::constants::math::PI / 2.}, {120, -6., 6.}, {10, 0., 10.}, {10, 0., 10.}}}}, //note: axes 3 and 4 (the pT) are updated in the init() + {"hCorrel2DVsPtSignalRegionMCRecBkgBkg", stringMCReco + "BkgBkg" + stringSignal + +stringDeltaPhi + stringDeltaEta + stringPtD + stringPtDbar + "entries", {HistType::kTHnSparseD, {{32, -o2::constants::math::PI / 2., 3. * o2::constants::math::PI / 2.}, {120, -6., 6.}, {10, 0., 10.}, {10, 0., 10.}}}}, //note: axes 3 and 4 (the pT) are updated in the init() + {"hDeltaEtaPtIntSidebandsMCRec", stringMCReco + stringSideband + stringDeltaEta + "entries", {HistType::kTH1F, {{200, -10., 10.}}}}, + {"hDeltaPhiPtIntSidebandsMCRec", stringMCReco + stringSideband + stringDeltaPhi + "entries", {HistType::kTH1F, {{32, -o2::constants::math::PI / 2., 3. * o2::constants::math::PI / 2.}}}}, + {"hCorrel2DPtIntSidebandsMCRec", stringMCReco + stringSideband + stringDeltaPhi + stringDeltaEta + "entries", {HistType::kTH2F, {{32, -o2::constants::math::PI / 2., 3. * o2::constants::math::PI / 2.}, {200, -10., 10.}}}}, + {"hDeltaPtDDbarSidebandsMCRec", stringMCReco + stringSideband + stringDeltaPt + "entries", {HistType::kTH1F, {{144, -36., 36.}}}}, + {"hDeltaPtMaxMinSidebandsMCRec", stringMCReco + stringSideband + stringDeltaPtMaxMin + "entries", {HistType::kTH1F, {{72, 0., 36.}}}}, + {"hCorrel2DVsPtSidebandsMCRecSigSig", stringMCReco + "SigSig" + stringSideband + stringDeltaPhi + stringDeltaEta + stringPtD + stringPtDbar + "entries", {HistType::kTHnSparseD, {{32, -o2::constants::math::PI / 2., 3. * o2::constants::math::PI / 2.}, {120, -6., 6.}, {10, 0., 10.}, {10, 0., 10.}}}}, //note: axes 3 and 4 (the pT) are updated in the init() - should be empty, kept for cross-check and debug + {"hCorrel2DVsPtSidebandsMCRecSigBkg", stringMCReco + "SigBkg" + stringSideband + stringDeltaPhi + stringDeltaEta + stringPtD + stringPtDbar + "entries", {HistType::kTHnSparseD, {{32, -o2::constants::math::PI / 2., 3. * o2::constants::math::PI / 2.}, {120, -6., 6.}, {10, 0., 10.}, {10, 0., 10.}}}}, //note: axes 3 and 4 (the pT) are updated in the init() + {"hCorrel2DVsPtSidebandsMCRecBkgSig", stringMCReco + "BkgSig" + stringSideband + stringDeltaPhi + stringDeltaEta + stringPtD + stringPtDbar + "entries", {HistType::kTHnSparseD, {{32, -o2::constants::math::PI / 2., 3. * o2::constants::math::PI / 2.}, {120, -6., 6.}, {10, 0., 10.}, {10, 0., 10.}}}}, //note: axes 3 and 4 (the pT) are updated in the init() + {"hCorrel2DVsPtSidebandsMCRecBkgBkg", stringMCReco + "BkgBkg" + stringSideband + stringDeltaPhi + stringDeltaEta + stringPtD + stringPtDbar + "entries", {HistType::kTHnSparseD, {{32, -o2::constants::math::PI / 2., 3. * o2::constants::math::PI / 2.}, {120, -6., 6.}, {10, 0., 10.}, {10, 0., 10.}}}}}}; //note: axes 3 and 4 (the pT) are updated in the init() + + //pT ranges for correlation plots: the default values are those embedded in hf_cuts_d0_topik (i.e. the mass pT bins), but can be redefined via json files + Configurable> binsCorrelations{"ptBinsForCorrelations", std::vector{o2::analysis::hf_cuts_d0_topik::pTBins_v}, "pT bin limits for correlation plots"}; + //signal and sideband region edges, to be defined via json file (initialised to empty) + Configurable> signalRegionInner{"signalRegionInner", std::vector(), "Inner values of signal region vs pT"}; + Configurable> signalRegionOuter{"signalRegionOuter", std::vector(), "Outer values of signal region vs pT"}; + Configurable> sidebandLeftInner{"sidebandLeftInner", std::vector(), "Inner values of left sideband vs pT"}; + Configurable> sidebandLeftOuter{"sidebandLeftOuter", std::vector(), "Outer values of left sideband vs pT"}; + Configurable> sidebandRightInner{"sidebandRightInner", std::vector(), "Inner values of right sideband vs pT"}; + Configurable> sidebandRightOuter{"sidebandRightOuter", std::vector(), "Outer values of right sideband vs pT"}; + + void init(o2::framework::InitContext&) + { + // redefinition of pT axes for THnSparse holding correlation entries + int nBinspTaxis = binsCorrelations->size() - 1; + const double* valuespTaxis = binsCorrelations->data(); + + for (int i = 2; i <= 3; i++) { + registry.get(HIST("hMass2DCorrelationPairsMCRecSigSig"))->GetAxis(i)->Set(nBinspTaxis, valuespTaxis); + registry.get(HIST("hMass2DCorrelationPairsMCRecSigBkg"))->GetAxis(i)->Set(nBinspTaxis, valuespTaxis); + registry.get(HIST("hMass2DCorrelationPairsMCRecBkgSig"))->GetAxis(i)->Set(nBinspTaxis, valuespTaxis); + registry.get(HIST("hMass2DCorrelationPairsMCRecBkgBkg"))->GetAxis(i)->Set(nBinspTaxis, valuespTaxis); + registry.get(HIST("hCorrel2DVsPtSignalRegionMCRecSigSig"))->GetAxis(i)->Set(nBinspTaxis, valuespTaxis); + registry.get(HIST("hCorrel2DVsPtSignalRegionMCRecSigBkg"))->GetAxis(i)->Set(nBinspTaxis, valuespTaxis); + registry.get(HIST("hCorrel2DVsPtSignalRegionMCRecBkgSig"))->GetAxis(i)->Set(nBinspTaxis, valuespTaxis); + registry.get(HIST("hCorrel2DVsPtSignalRegionMCRecBkgBkg"))->GetAxis(i)->Set(nBinspTaxis, valuespTaxis); + registry.get(HIST("hCorrel2DVsPtSidebandsMCRecSigSig"))->GetAxis(i)->Set(nBinspTaxis, valuespTaxis); + registry.get(HIST("hCorrel2DVsPtSidebandsMCRecSigBkg"))->GetAxis(i)->Set(nBinspTaxis, valuespTaxis); + registry.get(HIST("hCorrel2DVsPtSidebandsMCRecBkgSig"))->GetAxis(i)->Set(nBinspTaxis, valuespTaxis); + registry.get(HIST("hCorrel2DVsPtSidebandsMCRecBkgBkg"))->GetAxis(i)->Set(nBinspTaxis, valuespTaxis); + } + } + + void process(aod::DDbarPairFull const& pairEntries) + { + for (auto& pairEntry : pairEntries) { + //define variables for widely used quantities + double deltaPhi = pairEntry.deltaPhi(); + double deltaEta = pairEntry.deltaEta(); + double ptD = pairEntry.ptD(); + double ptDbar = pairEntry.ptDbar(); + double massD = pairEntry.mD(); + double massDbar = pairEntry.mDbar(); + + //reject entries outside pT ranges of interest + int pTBinD = o2::analysis::findBin(binsCorrelations, ptD); + int pTBinDbar = o2::analysis::findBin(binsCorrelations, ptDbar); + if (pTBinD == -1 || pTBinDbar == -1) { //at least one particle outside accepted pT range + continue; + } + + //fill 2D invariant mass plots + switch (pairEntry.signalStatus()) { + case 0: //D Bkg, Dbar Bkg + registry.fill(HIST("hMass2DCorrelationPairsMCRecBkgBkg"), massD, massDbar, ptD, ptDbar); + break; + case 1: //D Bkg, Dbar Sig + registry.fill(HIST("hMass2DCorrelationPairsMCRecBkgSig"), massD, massDbar, ptD, ptDbar); + break; + case 2: //D Sig, Dbar Bkg + registry.fill(HIST("hMass2DCorrelationPairsMCRecSigBkg"), massD, massDbar, ptD, ptDbar); + break; + case 3: //D Sig, Dbar Sig + registry.fill(HIST("hMass2DCorrelationPairsMCRecSigSig"), massD, massDbar, ptD, ptDbar); + break; + default: //should not happen for MC reco + break; + } + + //check if correlation entry belongs to signal region, sidebands or is outside both, and fill correlation plots + if (massD > signalRegionInner->at(pTBinD) && massD < signalRegionOuter->at(pTBinD) && massDbar > signalRegionInner->at(pTBinDbar) && massDbar < signalRegionOuter->at(pTBinDbar)) { + //in signal region + registry.fill(HIST("hCorrel2DPtIntSignalRegionMCRec"), deltaPhi, deltaEta); + registry.fill(HIST("hDeltaEtaPtIntSignalRegionMCRec"), deltaEta); + registry.fill(HIST("hDeltaPhiPtIntSignalRegionMCRec"), deltaPhi); + registry.fill(HIST("hDeltaPtDDbarSignalRegionMCRec"), ptDbar - ptD); + registry.fill(HIST("hDeltaPtMaxMinSignalRegionMCRec"), std::abs(ptDbar - ptD)); + switch (pairEntry.signalStatus()) { + case 0: //D Bkg, Dbar Bkg + registry.fill(HIST("hCorrel2DVsPtSignalRegionMCRecBkgBkg"), deltaPhi, deltaEta, ptD, ptDbar); + break; + case 1: //D Bkg, Dbar Sig + registry.fill(HIST("hCorrel2DVsPtSignalRegionMCRecBkgSig"), deltaPhi, deltaEta, ptD, ptDbar); + break; + case 2: //D Sig, Dbar Bkg + registry.fill(HIST("hCorrel2DVsPtSignalRegionMCRecSigBkg"), deltaPhi, deltaEta, ptD, ptDbar); + break; + case 3: //D Sig, Dbar Sig + registry.fill(HIST("hCorrel2DVsPtSignalRegionMCRecSigSig"), deltaPhi, deltaEta, ptD, ptDbar); + break; + default: //should not happen for MC reco + break; + } + } + + if ((massD > sidebandLeftInner->at(pTBinD) && massD < sidebandLeftOuter->at(pTBinD) && massDbar > sidebandLeftInner->at(pTBinDbar) && massDbar < sidebandRightOuter->at(pTBinDbar)) || + (massD > sidebandRightInner->at(pTBinD) && massD < sidebandRightOuter->at(pTBinD) && massDbar > sidebandLeftInner->at(pTBinDbar) && massDbar < sidebandRightOuter->at(pTBinDbar)) || + (massD > sidebandLeftInner->at(pTBinD) && massD < sidebandRightOuter->at(pTBinD) && massDbar > sidebandLeftInner->at(pTBinDbar) && massDbar < sidebandLeftOuter->at(pTBinDbar)) || + (massD > sidebandLeftInner->at(pTBinD) && massD < sidebandRightOuter->at(pTBinD) && massDbar > sidebandRightInner->at(pTBinDbar) && massDbar < sidebandRightOuter->at(pTBinDbar))) { + //in sideband region + registry.fill(HIST("hCorrel2DPtIntSidebandsMCRec"), deltaPhi, deltaEta); + registry.fill(HIST("hDeltaEtaPtIntSidebandsMCRec"), deltaEta); + registry.fill(HIST("hDeltaPhiPtIntSidebandsMCRec"), deltaPhi); + registry.fill(HIST("hDeltaPtDDbarSidebandsMCRec"), ptDbar - ptD); + registry.fill(HIST("hDeltaPtMaxMinSidebandsMCRec"), std::abs(ptDbar - ptD)); + switch (pairEntry.signalStatus()) { + case 0: //D Bkg, Dbar Bkg + registry.fill(HIST("hCorrel2DVsPtSidebandsMCRecBkgBkg"), deltaPhi, deltaEta, ptD, ptDbar); + break; + case 1: //D Bkg, Dbar Sig + registry.fill(HIST("hCorrel2DVsPtSidebandsMCRecBkgSig"), deltaPhi, deltaEta, ptD, ptDbar); + break; + case 2: //D Sig, Dbar Bkg + registry.fill(HIST("hCorrel2DVsPtSidebandsMCRecSigBkg"), deltaPhi, deltaEta, ptD, ptDbar); + break; + case 3: //D Sig, Dbar Sig + registry.fill(HIST("hCorrel2DVsPtSidebandsMCRecSigSig"), deltaPhi, deltaEta, ptD, ptDbar); + break; + default: //should not happen for MC reco + break; + } + } + } //end loop + } +}; + +/// D-Dbar correlation pair filling task, from pair tables - for MC gen-level analysis (no filter/selection, only true signal) - Ok for both USL and LS analyses +/// Works on both USL and LS analyses pair tables (and if tables are filled with quark pairs as well) +struct HfTaskCorrelationDDbarMcGen { + + HistogramRegistry registry{ + "registry", + {{"hDeltaEtaPtIntMCGen", stringMCParticles + stringDeltaEta + "entries", {HistType::kTH1F, {{200, -10., 10.}}}}, + {"hDeltaPhiPtIntMCGen", stringMCParticles + stringDeltaPhi + "entries", {HistType::kTH1F, {{32, -o2::constants::math::PI / 2., 3. * o2::constants::math::PI / 2.}}}}, + {"hCorrel2DPtIntMCGen", stringMCParticles + stringDeltaPhi + stringDeltaEta + "entries", {HistType::kTH2F, {{32, -o2::constants::math::PI / 2., 3. * o2::constants::math::PI / 2.}, {200, -10., 10.}}}}, + {"hCorrel2DVsPtMCGen", stringMCParticles + stringDeltaPhi + stringDeltaEta + stringPtD + stringPtDbar + "entries", {HistType::kTHnSparseD, {{32, -o2::constants::math::PI / 2., 3. * o2::constants::math::PI / 2.}, {120, -6., 6.}, {10, 0., 10.}, {10, 0., 10.}}}}, //note: axes 3 and 4 (the pT) are updated in the init() + {"hDeltaPtDDbarMCGen", stringMCParticles + stringDeltaPt + "entries", {HistType::kTH1F, {{144, -36., 36.}}}}, + {"hDeltaPtMaxMinMCGen", stringMCParticles + stringDeltaPtMaxMin + "entries", {HistType::kTH1F, {{72, 0., 36.}}}}}}; + + //pT ranges for correlation plots: the default values are those embedded in hf_cuts_d0_topik (i.e. the mass pT bins), but can be redefined via json files + Configurable> binsCorrelations{"ptBinsForCorrelations", std::vector{o2::analysis::hf_cuts_d0_topik::pTBins_v}, "pT bin limits for correlation plots"}; + + void init(o2::framework::InitContext&) + { + // redefinition of pT axes for THnSparse holding correlation entries + int nBinspTaxis = binsCorrelations->size() - 1; + const double* valuespTaxis = binsCorrelations->data(); + + for (int i = 2; i <= 3; i++) { + registry.get(HIST("hCorrel2DVsPtMCGen"))->GetAxis(i)->Set(nBinspTaxis, valuespTaxis); + } + } + + void process(aod::DDbarPair const& pairEntries) + { + for (auto& pairEntry : pairEntries) { + //define variables for widely used quantities + double deltaPhi = pairEntry.deltaPhi(); + double deltaEta = pairEntry.deltaEta(); + double ptD = pairEntry.ptD(); + double ptDbar = pairEntry.ptDbar(); + + //reject entries outside pT ranges of interest + if (o2::analysis::findBin(binsCorrelations, ptD) == -1 || o2::analysis::findBin(binsCorrelations, ptDbar) == -1) { + continue; + } + + registry.fill(HIST("hCorrel2DVsPtMCGen"), deltaPhi, deltaEta, ptD, ptDbar); + registry.fill(HIST("hCorrel2DPtIntMCGen"), deltaPhi, deltaEta); + registry.fill(HIST("hDeltaEtaPtIntMCGen"), deltaEta); + registry.fill(HIST("hDeltaPhiPtIntMCGen"), deltaPhi); + registry.fill(HIST("hDeltaPtDDbarMCGen"), ptDbar - ptD); + registry.fill(HIST("hDeltaPtMaxMinMCGen"), std::abs(ptDbar - ptD)); + } //end loop + } +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + WorkflowSpec workflow{}; + //MC-based tasks + const bool doMCGen = cfgc.options().get("doMCGen"); + const bool doMCRec = cfgc.options().get("doMCRec"); + if (doMCGen) { //MC-Gen analysis + workflow.push_back(adaptAnalysisTask(cfgc)); + } else if (doMCRec) { //MC-Rec analysis + workflow.push_back(adaptAnalysisTask(cfgc)); + } else { //data analysis + workflow.push_back(adaptAnalysisTask(cfgc)); + } + return workflow; +}