diff --git a/Analysis/DataModel/include/AnalysisDataModel/HFSecondaryVertex.h b/Analysis/DataModel/include/AnalysisDataModel/HFSecondaryVertex.h index 8c62fd41b43ba..03f26ad75c0f2 100644 --- a/Analysis/DataModel/include/AnalysisDataModel/HFSecondaryVertex.h +++ b/Analysis/DataModel/include/AnalysisDataModel/HFSecondaryVertex.h @@ -487,6 +487,26 @@ DECLARE_SOA_TABLE(HfCandProng3MCGen, "AOD", "HFCANDP3MCGEN", hf_cand_prong3::OriginMCGen, hf_cand_prong3::FlagMCDecayChanGen); +// definition of columns and tables for D0-D0bar correlation pairs +namespace hf_d0d0bar_correlation +{ +DECLARE_SOA_COLUMN(DeltaPhi, deltaPhi, float); +DECLARE_SOA_COLUMN(DeltaEta, deltaEta, float); +DECLARE_SOA_COLUMN(PtD0, ptD0, float); +DECLARE_SOA_COLUMN(PtD0bar, ptD0bar, float); +DECLARE_SOA_COLUMN(MD0, mD0, float); +DECLARE_SOA_COLUMN(MD0bar, mD0bar, float); +DECLARE_SOA_COLUMN(SignalStatus, signalStatus, int); +} // namespace hf_d0d0bar_correlation +DECLARE_SOA_TABLE(D0D0barPair, "AOD", "D0D0BARPAIR", + aod::hf_d0d0bar_correlation::DeltaPhi, + aod::hf_d0d0bar_correlation::DeltaEta, + aod::hf_d0d0bar_correlation::PtD0, + aod::hf_d0d0bar_correlation::PtD0bar); +DECLARE_SOA_TABLE(D0D0barRecoInfo, "AOD", "D0D0BARRECOINFO", + aod::hf_d0d0bar_correlation::MD0, + aod::hf_d0d0bar_correlation::MD0bar, + aod::hf_d0d0bar_correlation::SignalStatus); } // namespace o2::aod #endif // O2_ANALYSIS_HFSECONDARYVERTEX_H_ diff --git a/Analysis/Tasks/PWGHF/CMakeLists.txt b/Analysis/Tasks/PWGHF/CMakeLists.txt index 6fb285562c45d..58e51f2421987 100644 --- a/Analysis/Tasks/PWGHF/CMakeLists.txt +++ b/Analysis/Tasks/PWGHF/CMakeLists.txt @@ -108,6 +108,16 @@ o2_add_dpl_workflow(hf-task-x PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisDataModel O2::AnalysisCore O2::DetectorsVertexing COMPONENT_NAME Analysis) +o2_add_dpl_workflow(hf-d0d0bar-correlator + SOURCES D0D0barCorrelator.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisDataModel O2::AnalysisCore O2::DetectorsVertexing + COMPONENT_NAME Analysis) + +o2_add_dpl_workflow(hf-task-d0d0bar-correlation + SOURCES taskD0D0barCorrelation.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisDataModel O2::AnalysisCore O2::DetectorsVertexing + COMPONENT_NAME Analysis) + o2_add_dpl_workflow(hf-mc-validation SOURCES HFMCValidation.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisDataModel O2::AnalysisCore O2::DetectorsVertexing diff --git a/Analysis/Tasks/PWGHF/D0D0barCorrelator.cxx b/Analysis/Tasks/PWGHF/D0D0barCorrelator.cxx new file mode 100644 index 0000000000000..6ea037906b20e --- /dev/null +++ b/Analysis/Tasks/PWGHF/D0D0barCorrelator.cxx @@ -0,0 +1,848 @@ +// 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 taskD0D0barCorrelation.cxx +/// \brief D0-D0bar 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 "AnalysisCore/HFSelectorCuts.h" +#include "AnalysisDataModel/HFSecondaryVertex.h" +#include "AnalysisDataModel/HFCandidateSelectionTables.h" + +using namespace o2; +using namespace o2::framework; +using namespace o2::aod::hf_cand_prong2; +using namespace o2::analysis::hf_cuts_d0_topik; +using namespace o2::framework::expressions; +using namespace o2::constants::math; +using namespace o2::aod::hf_d0d0bar_correlation; + +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 D0D0barCorrelator { + 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, {{5, -0.5, 4.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 dSelectionFlagD0{"dSelectionFlagD0", 1, "Selection Flag for D0"}; + Configurable dSelectionFlagD0bar{"dSelectionFlagD0bar", 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 >= dSelectionFlagD0 || aod::hf_selcandidate_d0::isSelD0bar >= dSelectionFlagD0bar); + + 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. && std::abs(candidate1.pt()) < cutPtCandMin) { + continue; + } + //check decay channel flag for candidate1 + if (!(candidate1.hfflag() & 1 << D0ToPiK)) { + continue; + } + //fill invariant mass plots and generic info from all D0/D0bar candidates + if (candidate1.isSelD0() >= dSelectionFlagD0) { + registry.fill(HIST("hMass"), InvMassD0(candidate1), candidate1.pt()); + registry.fill(HIST("hMassD0"), InvMassD0(candidate1), candidate1.pt()); + } + if (candidate1.isSelD0bar() >= dSelectionFlagD0bar) { + 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.isSelD0() + (candidate1.isSelD0bar() * 2)); + + //D-Dbar correlation dedicated section + //if the candidate is a D0, search for D0bar and evaluate correlations + if (candidate1.isSelD0() >= dSelectionFlagD0) { + for (auto& candidate2 : candidates) { + //check decay channel flag for candidate2 + if (!(candidate2.hfflag() & 1 << D0ToPiK)) { + continue; + } + if (candidate2.isSelD0bar() >= dSelectionFlagD0bar) { + if (cutYCandMax >= 0. && std::abs(YD0(candidate2)) > cutYCandMax) { + continue; + } + if (cutPtCandMin >= 0. && std::abs(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 D0D0barCorrelatorMCRec { + + 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, {{5, -0.5, 4.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 dSelectionFlagD0{"dSelectionFlagD0", 1, "Selection Flag for D0"}; + Configurable dSelectionFlagD0bar{"dSelectionFlagD0bar", 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 >= dSelectionFlagD0 || aod::hf_selcandidate_d0::isSelD0bar >= dSelectionFlagD0bar); + + void process(aod::Collision const& collision, soa::Filtered> const& candidates) + { + //MC reco level + bool flagD0Signal = kFALSE; + bool flagD0barSignal = kFALSE; + for (auto& candidate1 : candidates) { + //check decay channel flag for candidate1 + if (!(candidate1.hfflag() & 1 << D0ToPiK)) { + continue; + } + if (cutYCandMax >= 0. && std::abs(YD0(candidate1)) > cutYCandMax) { + continue; + } + if (cutPtCandMin >= 0. && std::abs(candidate1.pt()) < cutPtCandMin) { + continue; + } + if (std::abs(candidate1.flagMCMatchRec()) == 1 << D0ToPiK) { + //fill invariant mass plots and generic info from all D0/D0bar candidates + if (candidate1.isSelD0() >= dSelectionFlagD0 && candidate1.flagMCMatchRec() == 1 << D0ToPiK) { //only reco and matched as D0 + registry.fill(HIST("hMassD0MCRec"), InvMassD0(candidate1), candidate1.pt()); + } + if (candidate1.isSelD0bar() >= dSelectionFlagD0bar && candidate1.flagMCMatchRec() == -1 << 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.isSelD0() + (candidate1.isSelD0bar() * 2)); + } + + //D-Dbar correlation dedicated section + //if the candidate is selected ad D0, search for D0bar and evaluate correlations + if (candidate1.isSelD0() < dSelectionFlagD0) { //discard candidates not selected as D0 in outer loop + continue; + } + if (candidate1.flagMCMatchRec() == 1 << D0ToPiK) { //candidate matched to D0 (particle) + flagD0Signal = kTRUE; + } else { //candidate of bkg, wrongly selected as D0 + flagD0Signal = kFALSE; + } + for (auto& candidate2 : candidates) { + if (!(candidate2.hfflag() & 1 << D0ToPiK)) { //check decay channel flag for candidate2 + continue; + } + if (candidate2.isSelD0bar() < dSelectionFlagD0bar) { //discard candidates not selected as D0bar in inner loop + continue; + } + if (candidate2.flagMCMatchRec() == -1 << D0ToPiK) { //candidate matched to D0bar (antiparticle) + flagD0barSignal = kTRUE; + } else { //candidate of bkg, wrongly selected as D0bar + flagD0barSignal = kFALSE; + } + if (cutYCandMax >= 0. && std::abs(YD0(candidate2)) > cutYCandMax) { + continue; + } + if (cutPtCandMin >= 0. && std::abs(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 D0D0barCorrelatorMCGen { + + 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;#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()) != 421) { + continue; + } + if (cutYCandMax >= 0. && std::abs(RecoDecay::Y(array{particle1.px(), particle1.py(), particle1.pz()}, RecoDecay::getMassPDG(particle1.pdgCode()))) > cutYCandMax) { + continue; + } + if (cutPtCandMin >= 0. && std::abs(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"), RecoDecay::Y(array{particle1.px(), particle1.py(), particle1.pz()}, RecoDecay::getMassPDG(particle1.pdgCode()))); + counterD0D0bar++; + + //D-Dbar correlation dedicated section + //if it's a D0 particle, search for D0bar and evaluate correlations + if (particle1.pdgCode() == 421) { //just checking the particle PDG, not the decay channel (differently from Reco: you have a BR factor btw such levels!) + registry.fill(HIST("hcountD0triggersMCGen"), 0, particle1.pt()); //to count trigger D0 (for normalisation) + for (auto& particle2 : particlesMC) { + if (particle2.pdgCode() == -421) { + if (cutYCandMax >= 0. && std::abs(RecoDecay::Y(array{particle2.px(), particle2.py(), particle2.pz()}, RecoDecay::getMassPDG(particle2.pdgCode()))) > cutYCandMax) { + continue; + } + if (cutPtCandMin >= 0. && std::abs(particle2.pt()) < cutPtCandMin) { + continue; + } + entryD0D0barPair(getDeltaPhi(particle2.phi(), particle1.phi()), + particle2.eta() - particle1.eta(), + particle1.pt(), + particle2.pt()); + 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(particle1.eta()) < etaCut && std::abs(particle2.eta()) < etaCut && particle1.pt() > ptCut && particle2.pt() > ptCut) + registry.fill(HIST("hDDbarVsEtaCut"), etaCut - epsilon, ptCut + epsilon); + ptCut += incrementPtThreshold; + } while (ptCut < ptThresholdForMaxEtaCut - epsilon); + } while (etaCut < maxEtaCut - epsilon); + } // end D0bar check + } //end inner loop + } //end D0 check + + } //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) +struct D0D0barCorrelatorLS { + + 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, {{5, -0.5, 4.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 dSelectionFlagD0{"dSelectionFlagD0", 1, "Selection Flag for D0"}; + Configurable dSelectionFlagD0bar{"dSelectionFlagD0bar", 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 >= dSelectionFlagD0 || aod::hf_selcandidate_d0::isSelD0bar >= dSelectionFlagD0bar); + + void process(aod::Collision const& collision, soa::Filtered> const& candidates) + { + for (auto& candidate1 : candidates) { + //check decay channel flag for candidate1 + if (!(candidate1.hfflag() & 1 << D0ToPiK)) { + continue; + } + if (cutYCandMax >= 0. && std::abs(YD0(candidate1)) > cutYCandMax) { + continue; + } + if (cutPtCandMin >= 0. && std::abs(candidate1.pt()) < cutPtCandMin) { + continue; + } + //fill invariant mass plots and generic info from all D0/D0bar candidates + if (candidate1.isSelD0() >= dSelectionFlagD0) { + registry.fill(HIST("hMass"), InvMassD0(candidate1), candidate1.pt()); + registry.fill(HIST("hMassD0"), InvMassD0(candidate1), candidate1.pt()); + } + if (candidate1.isSelD0bar() >= dSelectionFlagD0bar) { + 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.isSelD0() + (candidate1.isSelD0bar() * 2)); + + double ptParticle1 = candidate1.pt(); //trigger particle is the largest-pT one + + //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 << D0ToPiK)) { + continue; + } + //for the associated, has to have smaller pT, and pass D0sel if trigger passes D0sel, or D0barsel if trigger passes D0barsel + if (candidate2.pt() < ptParticle1 && ((candidate1.isSelD0() >= dSelectionFlagD0 && candidate2.isSelD0() >= dSelectionFlagD0) || (candidate1.isSelD0bar() >= dSelectionFlagD0bar && candidate2.isSelD0bar() >= dSelectionFlagD0bar))) { + if (cutYCandMax >= 0. && std::abs(YD0(candidate2)) > cutYCandMax) { + continue; + } + if (cutPtCandMin >= 0. && std::abs(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; + } + 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) +struct D0D0barCorrelatorMCRecLS { + + 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, {{5, -0.5, 4.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 dSelectionFlagD0{"dSelectionFlagD0", 1, "Selection Flag for D0"}; + Configurable dSelectionFlagD0bar{"dSelectionFlagD0bar", 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 >= dSelectionFlagD0 || aod::hf_selcandidate_d0::isSelD0bar >= dSelectionFlagD0bar); + + 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 << D0ToPiK)) { + continue; + } + if (cutYCandMax >= 0. && std::abs(YD0(candidate1)) > cutYCandMax) { + continue; + } + if (cutPtCandMin >= 0. && std::abs(candidate1.pt()) < cutPtCandMin) { + continue; + } + if (std::abs(candidate1.flagMCMatchRec()) == 1 << D0ToPiK) { + //fill invariant mass plots and generic info from all D0/D0bar candidates + if (candidate1.isSelD0() >= dSelectionFlagD0 && candidate1.flagMCMatchRec() == D0ToPiK) { //only reco and matched as D0 + registry.fill(HIST("hMassD0MCRec"), InvMassD0(candidate1)); + } + if (candidate1.isSelD0bar() >= dSelectionFlagD0bar && candidate1.flagMCMatchRec() == 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)); + + double ptParticle1 = candidate1.pt(); //trigger particle is the largest pT one + + //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 << D0ToPiK)) { + continue; + } + bool conditionLSForD0 = (candidate1.isSelD0() >= dSelectionFlagD0bar && candidate1.flagMCMatchRec() == 1 << D0ToPiK) && (candidate2.isSelD0() >= dSelectionFlagD0bar && candidate2.flagMCMatchRec() == 1 << D0ToPiK); + bool conditionLSForD0bar = (candidate1.isSelD0bar() >= dSelectionFlagD0bar && candidate1.flagMCMatchRec() == -1 << D0ToPiK) && (candidate2.isSelD0bar() >= dSelectionFlagD0bar && candidate2.flagMCMatchRec() == -1 << D0ToPiK); + if (candidate2.pt() < ptParticle1 && (conditionLSForD0 || conditionLSForD0bar)) { //LS pair (of D0 or of D0bar) + pt2= 0. && std::abs(YD0(candidate2)) > cutYCandMax) { + continue; + } + if (cutPtCandMin >= 0. && std::abs(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; + } + 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 +struct D0D0barCorrelatorMCGenLS { + + 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()) != 421) { + continue; + } + if (cutYCandMax >= 0. && std::abs(RecoDecay::Y(array{particle1.px(), particle1.py(), particle1.pz()}, RecoDecay::getMassPDG(particle1.pdgCode()))) > cutYCandMax) { + continue; + } + if (cutPtCandMin >= 0. && std::abs(particle1.pt()) < cutPtCandMin) { + continue; + } + + double ptParticle1 = particle1.pt(); //trigger particle is the largest pT one + + registry.fill(HIST("hPtCandMCGen"), particle1.pt()); + registry.fill(HIST("hEtaMCGen"), particle1.eta()); + registry.fill(HIST("hPhiMCGen"), particle1.phi()); + registry.fill(HIST("hYMCGen"), RecoDecay::Y(array{particle1.px(), particle1.py(), particle1.pz()}, RecoDecay::getMassPDG(particle1.pdgCode()))); + 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()) != 421) { //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. && std::abs(particle2.pt()) < cutPtCandMin) { + continue; + } + if (particle2.pt() < ptParticle1 && particle2.pdgCode() == particle1.pdgCode()) { //like-sign condition (both 421 or both -421) 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; + } + 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 CCbarCorrelatorMCGen { + + Produces entryD0D0barPair; + + 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}}}}, + {"hcountD0D0barPerEvent", "D0,cbar 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", "c 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) + { + 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()) != 4) { //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) + if (cutYCandMax >= 0. && std::abs(RecoDecay::Y(array{particle1.px(), particle1.py(), particle1.pz()}, RecoDecay::getMassPDG(particle1.pdgCode()))) > cutYCandMax) { + continue; + } + if (cutPtCandMin >= 0. && std::abs(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"), RecoDecay::Y(array{particle1.px(), particle1.py(), particle1.pz()}, RecoDecay::getMassPDG(particle1.pdgCode()))); + 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() == 4) { + + registry.fill(HIST("hcountD0triggersMCGen"), 0, particle1.pt()); //to count trigger c quark (for normalisation) + + for (auto& particle2 : particlesMC) { + if (cutYCandMax >= 0. && std::abs(RecoDecay::Y(array{particle2.px(), particle2.py(), particle2.pz()}, RecoDecay::getMassPDG(particle2.pdgCode()))) > cutYCandMax) { + continue; + } + if (cutPtCandMin >= 0. && std::abs(particle2.pt()) < cutPtCandMin) { + continue; + } + if (particle2.pdgCode() == -4) { + //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() == -4) { + continue; + } + entryD0D0barPair(getDeltaPhi(particle2.phi(), particle1.phi()), + particle2.eta() - particle1.eta(), + particle1.pt(), + particle2.pt()); + } // end outer if (check cbar) + } // end inner loop + } //end outer if (check c) + } //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 CCbarCorrelatorMCGenLS { + + Produces entryD0D0barPair; + + 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}}}}, + {"hcountD0D0barPerEvent", "D0,cbar 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", "c 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) + { + 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()) != 4) { //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) + if (cutYCandMax >= 0. && std::abs(RecoDecay::Y(array{particle1.px(), particle1.py(), particle1.pz()}, RecoDecay::getMassPDG(particle1.pdgCode()))) > cutYCandMax) { + continue; + } + if (cutPtCandMin >= 0. && std::abs(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"), RecoDecay::Y(array{particle1.px(), particle1.py(), particle1.pz()}, RecoDecay::getMassPDG(particle1.pdgCode()))); + counterccbar++; //count if c or cbar don't come from themselves during fragmentation (after kinematic selection) + + //c-cbar correlation dedicated section + double ptParticle1 = particle1.pt(); //trigger particle is the largest pT one + registry.fill(HIST("hcountD0triggersMCGen"), 0, particle1.pt()); //to count trigger c quark (for normalisation) + + for (auto& particle2 : particlesMC) { + if (std::abs(particle2.pdgCode()) != 4) { //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. && std::abs(particle2.pt()) < cutPtCandMin) { + continue; + } + if (particle2.pt() < ptParticle1 && 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; + } + entryD0D0barPair(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, TaskName{"d0d0bar-correlator-mc-gen"})); + } else if (doMCRec) { //MC-Reco analysis + workflow.push_back(adaptAnalysisTask(cfgc, TaskName{"d0d0bar-correlator-mc-rec"})); + } else if (doMCccbar) { //MC-Reco analysis + workflow.push_back(adaptAnalysisTask(cfgc, TaskName{"ccbar-correlator-mc-gen"})); + } else { //data analysis + workflow.push_back(adaptAnalysisTask(cfgc, TaskName{"d0d0bar-correlator"})); + } + } else { //ike-sign analyses + if (doMCGen) { //MC-Gen analysis + workflow.push_back(adaptAnalysisTask(cfgc, TaskName{"d0d0bar-correlator-mc-gen-ls"})); + } else if (doMCRec) { //MC-Reco analysis + workflow.push_back(adaptAnalysisTask(cfgc, TaskName{"d0d0bar-correlator-mc-rec-ls"})); + } else if (doMCccbar) { //MC-Reco analysis + workflow.push_back(adaptAnalysisTask(cfgc, TaskName{"ccbar-correlator-mc-gen-ls"})); + } else { //data analysis + workflow.push_back(adaptAnalysisTask(cfgc, TaskName{"d0d0bar-correlator-ls"})); + } + } + + return workflow; +} diff --git a/Analysis/Tasks/PWGHF/taskD0D0barCorrelation.cxx b/Analysis/Tasks/PWGHF/taskD0D0barCorrelation.cxx new file mode 100644 index 0000000000000..c09db03c40085 --- /dev/null +++ b/Analysis/Tasks/PWGHF/taskD0D0barCorrelation.cxx @@ -0,0 +1,479 @@ +// 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 taskD0D0barCorrelation.cxx +/// \brief D0-D0bar 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 "AnalysisCore/HFSelectorCuts.h" +#include "AnalysisDataModel/HFSecondaryVertex.h" +#include "AnalysisDataModel/HFCandidateSelectionTables.h" + +using namespace o2; +using namespace o2::framework; +using namespace o2::aod::hf_cand_prong2; +using namespace o2::analysis::hf_cuts_d0_topik; +using namespace o2::framework::expressions; +using namespace o2::constants::math; +using namespace o2::aod::hf_d0d0bar_correlation; + +namespace o2::aod +{ +using D0D0barPairFull = 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 deltaPhi value in range [-pi, pi], for resolution distributions +/// +double getDeltaPhiForResolution(double phiD, double phiDbar) +{ + return RecoDecay::constrainAngle(phiDbar - phiD, -o2::constants::math::PI); +} + +/// +/// 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 +*/ + +/// D0-D0bar 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 TaskD0D0barCorrelation { + + HistogramRegistry registry{ + "registry", + //NOTE: use hMassD0 (from correlator task) for normalisation, and hMass2DCorrelationPairs for 2D-sideband-subtraction purposes + {{"hMass2DCorrelationPairs", "D0,D0bar candidates 2D;inv. mass D0 (#pi K) (GeV/#it{c}^{2});inv. mass D0bar (#pi K) (GeV/#it{c}^{2});#it{p}_{T}^{D0};#it{p}_{T}^{D0bar};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", "D0,D0bar candidates signal region;#it{#eta}^{D0bar}-#it{#eta}^{D0};entries", {HistType::kTH1F, {{200, -10., 10.}}}}, + {"hDeltaPhiPtIntSignalRegion", "D0,D0bar candidates signal region;#it{#varphi}^{D0bar}-#it{#varphi}^{D0};entries", {HistType::kTH1F, {{32, -o2::constants::math::PI / 2., 3. * o2::constants::math::PI / 2.}}}}, + {"hCorrel2DPtIntSignalRegion", "D0,D0bar candidates signal region;#it{#varphi}^{D0bar}-#it{#varphi}^{D0};#it{#eta}^{D0bar}-#it{#eta}^{D0};entries", {HistType::kTH2F, {{32, -o2::constants::math::PI / 2., 3. * o2::constants::math::PI / 2.}, {200, -10., 10.}}}}, + {"hCorrel2DVsPtSignalRegion", "D0,D0bar candidates signal region;#it{#varphi}^{D0bar}-#it{#varphi}^{D0};#it{#eta}^{D0bar}-#it{#eta}^{D0};#it{p}_{T}^{D0};#it{p}_{T}^{D0bar};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", "D0,D0bar candidates signal region;#it{p}_{T}^{D0bar}-#it{p}_{T}^{D0};entries", {HistType::kTH1F, {{144, -36., 36.}}}}, + {"hDeltaPtMaxMinSignalRegion", "D0,D0bar candidates signal region;#it{p}_{T}^{max}-#it{p}_{T}^{min};entries", {HistType::kTH1F, {{72, 0., 36.}}}}, + {"hDeltaEtaPtIntSidebands", "D0,D0bar candidates sidebands;#it{#eta}^{D0bar}-#it{#eta}^{D0};entries", {HistType::kTH1F, {{200, -10., 10.}}}}, + {"hDeltaPhiPtIntSidebands", "D0,D0bar candidates sidebands;#it{#varphi}^{D0bar}-#it{#varphi}^{D0};entries", {HistType::kTH1F, {{32, -o2::constants::math::PI / 2., 3. * o2::constants::math::PI / 2.}}}}, + {"hCorrel2DPtIntSidebands", "D0,D0bar candidates sidebands;#it{#varphi}^{D0bar}-#it{#varphi}^{D0};#it{#eta}^{D0bar}-#it{#eta}^{D0};entries", {HistType::kTH2F, {{32, -o2::constants::math::PI / 2., 3. * o2::constants::math::PI / 2.}, {200, -10., 10.}}}}, + {"hCorrel2DVsPtSidebands", "D0,D0bar candidates sidebands;#it{#varphi}^{D0bar}-#it{#varphi}^{D0};#it{#eta}^{D0bar}-#it{#eta}^{D0};#it{p}_{T}^{D0};#it{p}_{T}^{D0bar};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", "D0,D0bar candidates sidebands;#it{p}_{T}^{D0bar}-#it{p}_{T}^{D0};entries", {HistType::kTH1F, {{144, -36., 36.}}}}, + {"hDeltaPtMaxMinSidebands", "D0,D0bar candidates sidebands;#it{p}_{T}^{max}-#it{p}_{T}^{min};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::D0D0barPairFull const& pairEntries) + { + for (auto& pairEntry : pairEntries) { + //define variables for widely used quantities + double deltaPhi = pairEntry.deltaPhi(); + double deltaEta = pairEntry.deltaEta(); + double ptD0 = pairEntry.ptD0(); + double ptD0bar = pairEntry.ptD0bar(); + + //reject entries outside pT ranges of interest + double minPtAllowed = binsCorrelations->at(0); + double maxPtAllowed = binsCorrelations->at(binsCorrelations->size() - 1); + if (ptD0 < minPtAllowed || ptD0bar < minPtAllowed || ptD0 > maxPtAllowed || ptD0bar > maxPtAllowed) { + continue; + } + + //fill 2D invariant mass plots + registry.fill(HIST("hMass2DCorrelationPairs"), pairEntry.mD0(), pairEntry.mD0bar(), ptD0, ptD0bar); + + //check if correlation entry belongs to signal region, sidebands or is outside both, and fill correlation plots + int pTBinD0 = o2::analysis::findBin(binsCorrelations, ptD0); + int pTBinD0bar = o2::analysis::findBin(binsCorrelations, ptD0bar); + + if (pairEntry.mD0() > signalRegionInner->at(pTBinD0) && pairEntry.mD0() < signalRegionOuter->at(pTBinD0) && pairEntry.mD0bar() > signalRegionInner->at(pTBinD0bar) && pairEntry.mD0bar() < signalRegionOuter->at(pTBinD0bar)) { + //in signal region + registry.fill(HIST("hCorrel2DVsPtSignalRegion"), deltaPhi, deltaEta, ptD0, ptD0bar); + registry.fill(HIST("hCorrel2DPtIntSignalRegion"), deltaPhi, deltaEta); + registry.fill(HIST("hDeltaEtaPtIntSignalRegion"), deltaEta); + registry.fill(HIST("hDeltaPhiPtIntSignalRegion"), deltaPhi); + registry.fill(HIST("hDeltaPtDDbarSignalRegion"), ptD0bar - ptD0); + registry.fill(HIST("hDeltaPtMaxMinSignalRegion"), std::abs(ptD0bar - ptD0)); + } + + if ((pairEntry.mD0() > sidebandLeftInner->at(pTBinD0) && pairEntry.mD0() < sidebandLeftOuter->at(pTBinD0) && pairEntry.mD0bar() > sidebandLeftInner->at(pTBinD0bar) && pairEntry.mD0bar() < sidebandRightOuter->at(pTBinD0bar)) || + (pairEntry.mD0() > sidebandRightInner->at(pTBinD0) && pairEntry.mD0() < sidebandRightOuter->at(pTBinD0) && pairEntry.mD0bar() > sidebandLeftInner->at(pTBinD0bar) && pairEntry.mD0bar() < sidebandRightOuter->at(pTBinD0bar)) || + (pairEntry.mD0() > sidebandLeftInner->at(pTBinD0) && pairEntry.mD0() < sidebandRightOuter->at(pTBinD0) && pairEntry.mD0bar() > sidebandLeftInner->at(pTBinD0bar) && pairEntry.mD0bar() < sidebandLeftOuter->at(pTBinD0bar)) || + (pairEntry.mD0() > sidebandLeftInner->at(pTBinD0) && pairEntry.mD0() < sidebandRightOuter->at(pTBinD0) && pairEntry.mD0bar() > sidebandRightInner->at(pTBinD0bar) && pairEntry.mD0bar() < sidebandRightOuter->at(pTBinD0bar))) { + //in sideband region + registry.fill(HIST("hCorrel2DVsPtSidebands"), deltaPhi, deltaEta, ptD0, ptD0bar); + registry.fill(HIST("hCorrel2DPtIntSidebands"), deltaPhi, deltaEta); + registry.fill(HIST("hDeltaEtaPtIntSidebands"), deltaEta); + registry.fill(HIST("hDeltaPhiPtIntSidebands"), deltaPhi); + registry.fill(HIST("hDeltaPtDDbarSidebands"), ptD0bar - ptD0); + registry.fill(HIST("hDeltaPtMaxMinSidebands"), std::abs(ptD0bar - ptD0)); + } + } //end loop + } +}; + +/// D0-D0bar 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 TaskD0D0barCorrelationMCRec { + + HistogramRegistry registry{ + "registry", + //NOTE: use hMassD0 (from correlator task) for normalisation, and hMass2DCorrelationPairs for 2D-sideband-subtraction purposes + {{"hMass2DCorrelationPairsMCRecSigSig", "D0,D0bar candidates 2D SigSig - MC reco;inv. mass D0 (#pi K) (GeV/#it{c}^{2});inv. mass D0bar (#pi K) (GeV/#it{c}^{2});#it{p}_{T}^{D0};#it{p}_{T}^{D0bar};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", "D0,D0bar candidates 2D SigBkg - MC reco;inv. mass D0 (#pi K) (GeV/#it{c}^{2});inv. mass D0bar (#pi K) (GeV/#it{c}^{2});#it{p}_{T}^{D0};#it{p}_{T}^{D0bar};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", "D0,D0bar candidates 2D BkgSig - MC reco;inv. mass D0 (#pi K) (GeV/#it{c}^{2});inv. mass D0bar (#pi K) (GeV/#it{c}^{2});#it{p}_{T}^{D0};#it{p}_{T}^{D0bar};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", "D0,D0bar candidates 2D BkgBkg - MC reco;inv. mass D0 (#pi K) (GeV/#it{c}^{2});inv. mass D0bar (#pi K) (GeV/#it{c}^{2});#it{p}_{T}^{D0};#it{p}_{T}^{D0bar};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", "D0,D0bar candidates signal region - MC reco;#it{#eta}^{D0bar}-#it{#eta}^{D0};entries", {HistType::kTH1F, {{200, -10., 10.}}}}, + {"hDeltaPhiPtIntSignalRegionMCRec", "D0,D0bar candidates signal region - MC reco;#it{#varphi}^{D0bar}-#it{#varphi}^{D0};entries", {HistType::kTH1F, {{32, -o2::constants::math::PI / 2., 3. * o2::constants::math::PI / 2.}}}}, + {"hCorrel2DPtIntSignalRegionMCRec", "D0,D0bar candidates signal region - MC reco;#it{#varphi}^{D0bar}-#it{#varphi}^{D0};#it{#eta}^{D0bar}-#it{#eta}^{D0};entries", {HistType::kTH2F, {{32, -o2::constants::math::PI / 2., 3. * o2::constants::math::PI / 2.}, {200, -10., 10.}}}}, + {"hDeltaPtDDbarSignalRegionMCRec", "D0,D0bar candidates signal region - MC reco;#it{p}_{T}^{D0bar}-#it{p}_{T}^{D0};entries", {HistType::kTH1F, {{144, -36., 36.}}}}, + {"hDeltaPtMaxMinSignalRegionMCRec", "D0,D0bar candidates signal region - MC reco;#it{p}_{T}^{max}-#it{p}_{T}^{min};entries", {HistType::kTH1F, {{72, 0., 36.}}}}, + {"hCorrel2DVsPtSignalRegionMCRecSigSig", "D0,D0bar candidates signal region SigSig - MC reco;#it{#varphi}^{D0bar}-#it{#varphi}^{D0};#it{#eta}^{D0bar}-#it{#eta}^{D0};#it{p}_{T}^{D0};#it{p}_{T}^{D0bar};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", "D0,D0bar candidates signal region SigBkg - MC reco;#it{#varphi}^{D0bar}-#it{#varphi}^{D0};#it{#eta}^{D0bar}-#it{#eta}^{D0};#it{p}_{T}^{D0};#it{p}_{T}^{D0bar};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", "D0,D0bar candidates signal region BkgSig - MC reco;#it{#varphi}^{D0bar}-#it{#varphi}^{D0};#it{#eta}^{D0bar}-#it{#eta}^{D0};#it{p}_{T}^{D0};#it{p}_{T}^{D0bar};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", "D0,D0bar candidates signal region BkgBkg - MC reco;#it{#varphi}^{D0bar}-#it{#varphi}^{D0};#it{#eta}^{D0bar}-#it{#eta}^{D0};#it{p}_{T}^{D0};#it{p}_{T}^{D0bar};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", "D0,D0bar candidates sidebands - MC reco;#it{#eta}^{D0bar}-#it{#eta}^{D0};entries", {HistType::kTH1F, {{200, -10., 10.}}}}, + {"hDeltaPhiPtIntSidebandsMCRec", "D0,D0bar candidates sidebands - MC reco;#it{#varphi}^{D0bar}-#it{#varphi}^{D0};entries", {HistType::kTH1F, {{32, -o2::constants::math::PI / 2., 3. * o2::constants::math::PI / 2.}}}}, + {"hCorrel2DPtIntSidebandsMCRec", "D0,D0bar candidates sidebands - MC reco;#it{#varphi}^{D0bar}-#it{#varphi}^{D0};#it{#eta}^{D0bar}-#it{#eta}^{D0};entries", {HistType::kTH2F, {{32, -o2::constants::math::PI / 2., 3. * o2::constants::math::PI / 2.}, {200, -10., 10.}}}}, + {"hCorrel2DVsPtSidebandsMCRecSigSig", "D0,D0bar candidates sidebands SigSig - MC reco;#it{#varphi}^{D0bar}-#it{#varphi}^{D0};#it{#eta}^{D0bar}-#it{#eta}^{D0};#it{p}_{T}^{D0};#it{p}_{T}^{D0bar};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", "D0,D0bar candidates sidebands SigBkg - MC reco;#it{#varphi}^{D0bar}-#it{#varphi}^{D0};#it{#eta}^{D0bar}-#it{#eta}^{D0};#it{p}_{T}^{D0};#it{p}_{T}^{D0bar};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", "D0,D0bar candidates sidebands BkgSig - MC reco;#it{#varphi}^{D0bar}-#it{#varphi}^{D0};#it{#eta}^{D0bar}-#it{#eta}^{D0};#it{p}_{T}^{D0};#it{p}_{T}^{D0bar};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", "D0,D0bar candidates sidebands BkgBkg - MC reco;#it{#varphi}^{D0bar}-#it{#varphi}^{D0};#it{#eta}^{D0bar}-#it{#eta}^{D0};#it{p}_{T}^{D0};#it{p}_{T}^{D0bar};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() + {"hDeltaPtDDbarSidebandsMCRec", "D0,D0bar candidates signal region - MC reco;#it{p}_{T}^{D0bar}-#it{p}_{T}^{D0};entries", {HistType::kTH1F, {{144, -36., 36.}}}}, + {"hDeltaPtMaxMinSidebandsMCRec", "D0,D0bar candidates signal region - MC reco;#it{p}_{T}^{max}-#it{p}_{T}^{min};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("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::D0D0barPairFull const& pairEntries) + { + for (auto& pairEntry : pairEntries) { + //define variables for widely used quantities + double deltaPhi = pairEntry.deltaPhi(); + double deltaEta = pairEntry.deltaEta(); + double ptD0 = pairEntry.ptD0(); + double ptD0bar = pairEntry.ptD0bar(); + + //reject entries outside pT ranges of interest + double minPtAllowed = binsCorrelations->at(0); + double maxPtAllowed = binsCorrelations->at(binsCorrelations->size() - 1); + if (ptD0 < minPtAllowed || ptD0bar < minPtAllowed || ptD0 > maxPtAllowed || ptD0bar > maxPtAllowed) { + continue; + } + + //fill 2D invariant mass plots + switch (pairEntry.signalStatus()) { + case 0: //D0 Bkg, D0bar Bkg + registry.fill(HIST("hMass2DCorrelationPairsMCRecBkgBkg"), pairEntry.mD0(), pairEntry.mD0bar(), ptD0, ptD0bar); + break; + case 1: //D0 Bkg, D0bar Sig + registry.fill(HIST("hMass2DCorrelationPairsMCRecBkgSig"), pairEntry.mD0(), pairEntry.mD0bar(), ptD0, ptD0bar); + break; + case 2: //D0 Sig, D0bar Bkg + registry.fill(HIST("hMass2DCorrelationPairsMCRecSigBkg"), pairEntry.mD0(), pairEntry.mD0bar(), ptD0, ptD0bar); + break; + case 3: //D0 Sig, D0bar Sig + registry.fill(HIST("hMass2DCorrelationPairsMCRecSigSig"), pairEntry.mD0(), pairEntry.mD0bar(), ptD0, ptD0bar); + 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 + int pTBinD0 = o2::analysis::findBin(binsCorrelations, ptD0); + int pTBinD0bar = o2::analysis::findBin(binsCorrelations, ptD0bar); + + if (pairEntry.mD0() > signalRegionInner->at(pTBinD0) && pairEntry.mD0() < signalRegionOuter->at(pTBinD0) && pairEntry.mD0bar() > signalRegionInner->at(pTBinD0bar) && pairEntry.mD0bar() < signalRegionOuter->at(pTBinD0bar)) { + //in signal region + registry.fill(HIST("hCorrel2DPtIntSignalRegionMCRec"), deltaPhi, deltaEta); + registry.fill(HIST("hDeltaEtaPtIntSignalRegionMCRec"), deltaEta); + registry.fill(HIST("hDeltaPhiPtIntSignalRegionMCRec"), deltaPhi); + registry.fill(HIST("hDeltaPtDDbarSignalRegionMCRec"), ptD0bar - ptD0); + registry.fill(HIST("hDeltaPtMaxMinSignalRegionMCRec"), std::abs(ptD0bar - ptD0)); + switch (pairEntry.signalStatus()) { + case 0: //D0 Bkg, D0bar Bkg + registry.fill(HIST("hCorrel2DVsPtSignalRegionMCRecBkgBkg"), deltaPhi, deltaEta, ptD0, ptD0bar); + break; + case 1: //D0 Bkg, D0bar Sig + registry.fill(HIST("hCorrel2DVsPtSignalRegionMCRecBkgSig"), deltaPhi, deltaEta, ptD0, ptD0bar); + break; + case 2: //D0 Sig, D0bar Bkg + registry.fill(HIST("hCorrel2DVsPtSignalRegionMCRecSigBkg"), deltaPhi, deltaEta, ptD0, ptD0bar); + break; + case 3: //D0 Sig, D0bar Sig + registry.fill(HIST("hCorrel2DVsPtSignalRegionMCRecSigSig"), deltaPhi, deltaEta, ptD0, ptD0bar); + break; + default: //should not happen for MC reco + break; + } + } + + if ((pairEntry.mD0() > sidebandLeftInner->at(pTBinD0) && pairEntry.mD0() < sidebandLeftOuter->at(pTBinD0) && pairEntry.mD0bar() > sidebandLeftInner->at(pTBinD0bar) && pairEntry.mD0bar() < sidebandRightOuter->at(pTBinD0bar)) || + (pairEntry.mD0() > sidebandRightInner->at(pTBinD0) && pairEntry.mD0() < sidebandRightOuter->at(pTBinD0) && pairEntry.mD0bar() > sidebandLeftInner->at(pTBinD0bar) && pairEntry.mD0bar() < sidebandRightOuter->at(pTBinD0bar)) || + (pairEntry.mD0() > sidebandLeftInner->at(pTBinD0) && pairEntry.mD0() < sidebandRightOuter->at(pTBinD0) && pairEntry.mD0bar() > sidebandLeftInner->at(pTBinD0bar) && pairEntry.mD0bar() < sidebandLeftOuter->at(pTBinD0bar)) || + (pairEntry.mD0() > sidebandLeftInner->at(pTBinD0) && pairEntry.mD0() < sidebandRightOuter->at(pTBinD0) && pairEntry.mD0bar() > sidebandRightInner->at(pTBinD0bar) && pairEntry.mD0bar() < sidebandRightOuter->at(pTBinD0bar))) { + //in sideband region + registry.fill(HIST("hCorrel2DPtIntSidebandsMCRec"), deltaPhi, deltaEta); + registry.fill(HIST("hDeltaEtaPtIntSidebandsMCRec"), deltaEta); + registry.fill(HIST("hDeltaPhiPtIntSidebandsMCRec"), deltaPhi); + registry.fill(HIST("hDeltaPtDDbarSidebandsMCRec"), ptD0bar - ptD0); + registry.fill(HIST("hDeltaPtMaxMinSidebandsMCRec"), std::abs(ptD0bar - ptD0)); + switch (pairEntry.signalStatus()) { + case 0: //D0 Bkg, D0bar Bkg + registry.fill(HIST("hCorrel2DVsPtSidebandsMCRecBkgBkg"), deltaPhi, deltaEta, ptD0, ptD0bar); + break; + case 1: //D0 Bkg, D0bar Sig + registry.fill(HIST("hCorrel2DVsPtSidebandsMCRecBkgSig"), deltaPhi, deltaEta, ptD0, ptD0bar); + break; + case 2: //D0 Sig, D0bar Bkg + registry.fill(HIST("hCorrel2DVsPtSidebandsMCRecSigBkg"), deltaPhi, deltaEta, ptD0, ptD0bar); + break; + case 3: //D0 Sig, D0bar Sig + registry.fill(HIST("hCorrel2DVsPtSidebandsMCRecSigSig"), deltaPhi, deltaEta, ptD0, ptD0bar); + break; + default: //should not happen for MC reco + break; + } + } + } //end loop + } +}; + +/// D0-D0bar 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 TaskD0D0barCorrelationMCGen { + + HistogramRegistry registry{ + "registry", + {{"hDeltaEtaPtIntMCGen", "D0,D0bar particles - MC gen;#it{#eta}^{D0bar}-#it{#eta}^{D0};entries", {HistType::kTH1F, {{200, -10., 10.}}}}, + {"hDeltaPhiPtIntMCGen", "D0,D0bar particles - MC gen;#it{#varphi}^{D0bar}-#it{#varphi}^{D0};entries", {HistType::kTH1F, {{32, -o2::constants::math::PI / 2., 3. * o2::constants::math::PI / 2.}}}}, + {"hCorrel2DPtIntMCGen", "D0,D0bar particles - MC gen;#it{#varphi}^{D0bar}-#it{#varphi}^{D0};#it{#eta}^{D0bar}-#it{#eta}^{D0};entries", {HistType::kTH2F, {{32, -o2::constants::math::PI / 2., 3. * o2::constants::math::PI / 2.}, {200, -10., 10.}}}}, + {"hCorrel2DVsPtMCGen", "D0,D0bar particles - MC gen;#it{#varphi}^{D0bar}-#it{#varphi}^{D0};#it{#eta}^{D0bar}-#it{#eta}^{D0};#it{p}_{T}^{D0};#it{p}_{T}^{D0bar};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", "D0,D0bar particles - MC gen;#it{p}_{T}^{D0bar}-#it{p}_{T}^{D0};entries", {HistType::kTH1F, {{144, -36., 36.}}}}, + {"hDeltaPtMaxMinMCGen", "D0,D0bar particles - MC gen;#it{p}_{T}^{max}-#it{p}_{T}^{min};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::D0D0barPair const& pairEntries) + { + for (auto& pairEntry : pairEntries) { + //define variables for widely used quantities + double deltaPhi = pairEntry.deltaPhi(); + double deltaEta = pairEntry.deltaEta(); + double ptD0 = pairEntry.ptD0(); + double ptD0bar = pairEntry.ptD0bar(); + + //reject entries outside pT ranges of interest + double minPtAllowed = binsCorrelations->at(0); + double maxPtAllowed = binsCorrelations->at(binsCorrelations->size() - 1); + if (ptD0 < minPtAllowed || ptD0bar < minPtAllowed || ptD0 > maxPtAllowed || ptD0bar > maxPtAllowed) { + continue; + } + + registry.fill(HIST("hCorrel2DVsPtMCGen"), deltaPhi, deltaEta, ptD0, ptD0bar); + registry.fill(HIST("hCorrel2DPtIntMCGen"), deltaPhi, deltaEta); + registry.fill(HIST("hDeltaEtaPtIntMCGen"), deltaEta); + registry.fill(HIST("hDeltaPhiPtIntMCGen"), deltaPhi); + registry.fill(HIST("hDeltaPtDDbarMCGen"), ptD0bar - ptD0); + registry.fill(HIST("hDeltaPtMaxMinMCGen"), std::abs(ptD0bar - ptD0)); + } //end loop + } +}; + +/// checks phi resolution for standard definition and sec-vtx based definition +struct TaskD0D0barCorrelationCheckPhiResolution { + + HistogramRegistry registry{ + "registry", + {{"hMass", "D0,D0bar candidates;inv. mass (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH1F, {{120, 1.5848, 2.1848}}}}, + {"hEta", "D0,D0bar candidates;candidate #it{#eta};entries", {HistType::kTH1F, {{100, -5., 5.}}}}, + {"hPhiStdPhi", "D0,D0bar candidates;candidate #it{#varphi};#it{p}_{T}", {HistType::kTH2F, {{128, 0., 2. * o2::constants::math::PI}, {50, 0., 50.}}}}, + {"hPhiByVtxPhi", "D0,D0bar candidates;candidate #it{#varphi};#it{p}_{T}", {HistType::kTH2F, {{128, 0., 2. * o2::constants::math::PI}, {50, 0., 50.}}}}, + {"hPhiDifferenceTwoMethods", "D0,D0bar candidates;candidate #it{#Delta#varphi};#it{p}_{T}", {HistType::kTH2F, {{128, -o2::constants::math::PI, o2::constants::math::PI}, {50, 0., 50.}}}}, + {"hDifferenceGenPhiStdPhi", "D0,D0bar candidates;candidate #it{#varphi};#it{p}_{T}", {HistType::kTH2F, {{128, -o2::constants::math::PI, o2::constants::math::PI}, {50, 0., 50.}}}}, + {"hDifferenceGenPhiByVtxPhi", "D0,D0bar candidates;candidate #it{#varphi};#it{p}_{T}", {HistType::kTH2F, {{128, -o2::constants::math::PI, o2::constants::math::PI}, {50, 0., 50.}}}}, + {"hDeltaPhiPtIntStdPhi", "D0,D0bar candidates;#it{#varphi}^{D0bar}-#it{#varphi}^{D0};entries", {HistType::kTH1F, {{128, -o2::constants::math::PI / 2., 3. * o2::constants::math::PI / 2.}}}}, + {"hDeltaPhiPtIntByVtxPhi", "D0,D0bar candidates;#it{#varphi}^{D0bar}-#it{#varphi}^{D0};entries", {HistType::kTH1F, {{128, -o2::constants::math::PI / 2., 3. * o2::constants::math::PI / 2.}}}}, + {"hDeltaPhiVsPtStdPhi", "D0,D0bar candidates;#it{p}_{T}^{D0};#it{p}_{T}^{D0bar};#it{#varphi}^{D0bar}-#it{#varphi}^{D0};entries", {HistType::kTH3F, {{36, 0., 36.}, {36, 0., 36.}, {128, -o2::constants::math::PI / 2., 3. * o2::constants::math::PI / 2.}}}}, + {"hDeltaPhiVsPtByVtxPhi", "D0,D0bar candidates;#it{p}_{T}^{D0};#it{p}_{T}^{D0bar};#it{#varphi}^{D0bar}-#it{#varphi}^{D0};entries", {HistType::kTH3F, {{36, 0., 36.}, {36, 0., 36.}, {128, -o2::constants::math::PI / 2., 3. * o2::constants::math::PI / 2.}}}}}}; + + Configurable dSelectionFlagD0{"dSelectionFlagD0", 1, "Selection Flag for D0"}; + Configurable dSelectionFlagD0bar{"dSelectionFlagD0bar", 1, "Selection Flag for D0bar"}; + Configurable cutYCandMax{"cutYCandMax", -1., "max. cand. rapidity"}; + Configurable cutPtCandMin{"cutPtCandMin", -1., "min. cand. pT"}; + + Filter filterSelectCandidates = (aod::hf_selcandidate_d0::isSelD0 >= dSelectionFlagD0 || aod::hf_selcandidate_d0::isSelD0bar >= dSelectionFlagD0bar); + + void process(aod::Collision const& collision, soa::Filtered> const& candidates, aod::McParticles const& particlesMC, aod::BigTracksMC const& tracksMC) + { + for (auto& candidate1 : candidates) { + //check decay channel flag for candidate1 + if (!(candidate1.hfflag() & 1 << D0ToPiK)) { + continue; + } + if (cutYCandMax >= 0. && std::abs(YD0(candidate1)) > cutYCandMax) { + continue; + } + if (cutPtCandMin >= 0. && candidate1.pt() < cutPtCandMin) { + continue; + } + registry.fill(HIST("hMass"), InvMassD0(candidate1)); + registry.fill(HIST("hEta"), candidate1.eta()); + + //D-Dbar correlation dedicated section + //if it's a candidate D0, search for D0bar and evaluate correlations + if (candidate1.isSelD0() >= dSelectionFlagD0) { + double xPrimaryVertex = candidate1.index0_as().collision().posX(), yPrimaryVertex = candidate1.index0_as().collision().posY(); + double pt1 = candidate1.pt(), phi1Std = candidate1.phi(); + double phi1ByVtx = evaluatePhiByVertex(xPrimaryVertex, candidate1.xSecondaryVertex(), yPrimaryVertex, candidate1.ySecondaryVertex()); + registry.fill(HIST("hPhiStdPhi"), phi1Std, pt1); + registry.fill(HIST("hPhiByVtxPhi"), phi1ByVtx, pt1); + registry.fill(HIST("hPhiDifferenceTwoMethods"), getDeltaPhiForResolution(phi1ByVtx, phi1Std), pt1); + + //get corresponding gen-level D0, if exists, and evaluate gen-rec phi-difference with two approaches + if (std::abs(candidate1.flagMCMatchRec()) == 1 << D0ToPiK) { //ok to keep both D0 and D0bar + int indexGen = RecoDecay::getMother(particlesMC, candidate1.index0_as().mcParticle(), 421, true); //MC-gen corresponding index for MC-reco candidate + double phi1Gen = particlesMC.iteratorAt(indexGen).phi(); + registry.fill(HIST("hDifferenceGenPhiStdPhi"), getDeltaPhiForResolution(phi1Std, phi1Gen), pt1); + registry.fill(HIST("hDifferenceGenPhiByVtxPhi"), getDeltaPhiForResolution(phi1ByVtx, phi1Gen), pt1); + } + + for (auto& candidate2 : candidates) { + //check decay channel flag for candidate2 + if (!(candidate2.hfflag() & 1 << D0ToPiK)) { + continue; + } + if (candidate2.isSelD0bar() >= dSelectionFlagD0bar) { //accept only D0bar candidates + if (cutYCandMax >= 0. && std::abs(YD0(candidate2)) > cutYCandMax) { + continue; + } + if (cutPtCandMin >= 0. && candidate2.pt() < cutPtCandMin) { + continue; + } + //Excluding self-correlations (could happen in case of reflections) + if (candidate1.mRowIndex == candidate2.mRowIndex) { + continue; + } + double pt2 = candidate2.pt(), phi2Std = candidate2.phi(); + double phi2ByVtx = evaluatePhiByVertex(xPrimaryVertex, candidate2.xSecondaryVertex(), yPrimaryVertex, candidate2.ySecondaryVertex()); + registry.fill(HIST("hDeltaPhiPtIntStdPhi"), getDeltaPhi(phi2Std, phi1Std)); + registry.fill(HIST("hDeltaPhiPtIntByVtxPhi"), getDeltaPhi(phi2ByVtx, phi1ByVtx)); + registry.fill(HIST("hDeltaPhiVsPtStdPhi"), pt1, pt2, getDeltaPhi(phi2Std, phi1Std)); + registry.fill(HIST("hDeltaPhiVsPtByVtxPhi"), pt1, pt2, getDeltaPhi(phi2ByVtx, phi1ByVtx)); + } + } // end inner loop (Dbars) + } + } //end outer 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, TaskName{"task-d0d0bar-correlation-mc-gen"})); + } else if (doMCRec) { //MC-Rec analysis + workflow.push_back(adaptAnalysisTask(cfgc, TaskName{"task-d0d0bar-correlation-mc-rec"})); + workflow.push_back(adaptAnalysisTask(cfgc, TaskName{"task-d0d0bar-correlation-check-phi-resolution"})); + } else { //data analysis + workflow.push_back(adaptAnalysisTask(cfgc, TaskName{"task-d0d0bar-correlation"})); + } + return workflow; +}