diff --git a/ALICE3/DataModel/A3DecayFinderTables.h b/ALICE3/DataModel/A3DecayFinderTables.h index 05eb5ad9307..55229bbb5d4 100644 --- a/ALICE3/DataModel/A3DecayFinderTables.h +++ b/ALICE3/DataModel/A3DecayFinderTables.h @@ -20,6 +20,7 @@ // O2 includes #include "Framework/AnalysisDataModel.h" +#include "Common/Core/RecoDecay.h" enum a3selectionBit : uint32_t { kDCAxy = 0, kInnerTOFPion, @@ -59,6 +60,64 @@ DECLARE_SOA_TABLE(Alice3DecayMaps, "AOD", "ALICE3DECAYMAPS", using Alice3DecayMap = Alice3DecayMaps::iterator; +namespace a3D0meson +{ +DECLARE_SOA_INDEX_COLUMN(Collision, collision); //! +DECLARE_SOA_COLUMN(PxProng0, pxProng0, float); //! positive track +DECLARE_SOA_COLUMN(PyProng0, pyProng0, float); //! +DECLARE_SOA_COLUMN(PzProng0, pzProng0, float); //! +DECLARE_SOA_COLUMN(PxProng1, pxProng1, float); //! negative track +DECLARE_SOA_COLUMN(PyProng1, pyProng1, float); //! +DECLARE_SOA_COLUMN(PzProng1, pzProng1, float); //! +DECLARE_SOA_DYNAMIC_COLUMN(PtProng0, ptProng0, //! + [](float px, float py) -> float { return RecoDecay::pt(px, py); }); +DECLARE_SOA_DYNAMIC_COLUMN(PtProng1, ptProng1, //! + [](float px, float py) -> float { return RecoDecay::pt(px, py); }); +DECLARE_SOA_EXPRESSION_COLUMN(Px, px, //! + float, 1.f * aod::a3D0meson::pxProng0 + 1.f * aod::a3D0meson::pxProng1); +DECLARE_SOA_EXPRESSION_COLUMN(Py, py, //! + float, 1.f * aod::a3D0meson::pyProng0 + 1.f * aod::a3D0meson::pyProng1); +DECLARE_SOA_EXPRESSION_COLUMN(Pz, pz, //! + float, 1.f * aod::a3D0meson::pzProng0 + 1.f * aod::a3D0meson::pzProng1); +DECLARE_SOA_COLUMN(Pt, pt, float); //! +DECLARE_SOA_COLUMN(M, m, float); //! +DECLARE_SOA_DYNAMIC_COLUMN(E, e, //! + [](float px, float py, float pz, double m) -> float { return RecoDecay::e(px, py, pz, m); }); +DECLARE_SOA_COLUMN(Eta, eta, float); //! +DECLARE_SOA_COLUMN(Phi, phi, float); //! +DECLARE_SOA_COLUMN(Y, y, float); +} // namespace a3D0meson +DECLARE_SOA_TABLE(Alice3D0Meson, "AOD", "ALICE3D0MESON", //! + o2::soa::Index<>, + a3D0meson::CollisionId, + a3D0meson::PxProng0, a3D0meson::PyProng0, a3D0meson::PzProng0, // positive track + a3D0meson::PxProng1, a3D0meson::PyProng1, a3D0meson::PzProng1, // negative track + a3D0meson::PtProng0, + a3D0meson::PtProng1, + a3D0meson::Px, a3D0meson::Py, a3D0meson::Pz, + a3D0meson::Pt, + a3D0meson::M, + a3D0meson::E, + a3D0meson::Eta, + a3D0meson::Phi, + a3D0meson::Y); + +namespace a3D0Selection +{ +DECLARE_SOA_COLUMN(IsSelD0, isSelD0, int); //! +DECLARE_SOA_COLUMN(IsSelD0bar, isSelD0bar, int); //! +} // namespace a3D0Selection +DECLARE_SOA_TABLE(Alice3D0Sel, "AOD", "ALICE3D0SEL", //! + a3D0Selection::IsSelD0, + a3D0Selection::IsSelD0bar); + +namespace a3D0MCTruth +{ +DECLARE_SOA_COLUMN(McTruthInfo, mcTruthInfo, int); //! 0 for bkg, 1 for true D0, 2 for true D0bar +} // namespace a3D0MCTruth +DECLARE_SOA_TABLE(Alice3D0MCTruth, "AOD", "ALICE3D0MCTRUTH", //! + a3D0MCTruth::McTruthInfo); //! + } // namespace o2::aod #endif // ALICE3_DATAMODEL_A3DECAYFINDERTABLES_H_ diff --git a/ALICE3/TableProducer/CMakeLists.txt b/ALICE3/TableProducer/CMakeLists.txt index 9c3d710d358..8548ecd9897 100644 --- a/ALICE3/TableProducer/CMakeLists.txt +++ b/ALICE3/TableProducer/CMakeLists.txt @@ -45,3 +45,8 @@ o2physics_add_dpl_workflow(alice3-multicharm SOURCES alice3-multicharm.cxx PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2::DCAFitter COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(alice3-correlatorddbar + SOURCES alice3-correlatorDDbar.cxx + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2::DCAFitter + COMPONENT_NAME Analysis) diff --git a/ALICE3/TableProducer/alice3-correlatorDDbar.cxx b/ALICE3/TableProducer/alice3-correlatorDDbar.cxx new file mode 100644 index 00000000000..ad224442717 --- /dev/null +++ b/ALICE3/TableProducer/alice3-correlatorDDbar.cxx @@ -0,0 +1,315 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file alice3-correlatorDDbar.cxx +/// \brief D0-D0bar correlator task - data-like and MC-reco analysis performance on ALICE 3 detector. +/// +/// \author Fabio Colamaria , INFN Bari + +#include + +#include "CommonConstants/PhysicsConstants.h" +#include "Framework/AnalysisTask.h" +#include "Framework/HistogramRegistry.h" +#include "Framework/runDataProcessing.h" + +#include "Common/Core/TrackSelection.h" +#include "Common/DataModel/TrackSelectionTables.h" + +#include "ALICE3/DataModel/A3DecayFinderTables.h" + +#include "PWGHF/Core/HfHelper.h" +#include "PWGHF/HFC/DataModel/CorrelationTables.h" +// #include "PWGHF/DataModel/CandidateReconstructionTables.h" +// #include "PWGHF/DataModel/CandidateSelectionTables.h" + +using namespace o2; +using namespace o2::analysis; +using namespace o2::constants::physics; +using namespace o2::framework; +using namespace o2::framework::expressions; + +/// +/// Returns deltaPhi value in range [-pi/2., 3.*pi/2], typically used for correlation studies +/// +double getDeltaPhi(double phiD, double phiDbar) +{ + return RecoDecay::constrainAngle(phiDbar - phiD, -o2::constants::math::PIHalf); +} + +/// 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; + +const int npTBinsMassAndEfficiency = o2::analysis::hf_cuts_d0_to_pi_k::NBinsPt; +const double efficiencyDmesonDefault[npTBinsMassAndEfficiency] = {}; +auto efficiencyDmeson_v = std::vector{efficiencyDmesonDefault, efficiencyDmesonDefault + npTBinsMassAndEfficiency}; + +// histogram binning definition +const int massAxisBins = 120; +const double massAxisMin = 1.5848; +const double massAxisMax = 2.1848; +const int phiAxisBins = 32; +const double phiAxisMin = 0.; +const double phiAxisMax = o2::constants::math::TwoPI; +const int yAxisBins = 100; +const double yAxisMin = -5.; +const double yAxisMax = 5.; +const int ptDAxisBins = 180; +const double ptDAxisMin = 0.; +const double ptDAxisMax = 36.; + +struct alice3correlatorddbar { + SliceCache cache; + Preslice perCol = aod::a3D0meson::collisionId; + Produces entryD0D0barPair; + Produces entryD0D0barRecoInfo; + + Configurable selectionFlagD0{"selectionFlagD0", 1, "Selection Flag for D0"}; + Configurable selectionFlagD0bar{"selectionFlagD0bar", 1, "Selection Flag for D0bar"}; + Configurable applyEfficiency{"applyEfficiency", 1, "Flag for applying D-meson efficiency weights"}; + Configurable yCandMax{"yCandMax", 999., "max. cand. rapidity"}; + Configurable ptCandMin{"ptCandMin", -1., "min. cand. pT"}; + Configurable> binsPt{"binsPt", std::vector{o2::analysis::hf_cuts_d0_to_pi_k::vecBinsPt}, "pT bin limits for candidate mass plots and efficiency"}; + Configurable> efficiencyD{"efficiencyD", std::vector{efficiencyDmeson_v}, "Efficiency values for D0 meson"}; + + // HfHelper hfHelper; //not needed for now + + Partition> selectedCandidates = aod::a3D0meson::y > -yCandMax&& aod::a3D0meson::y ptCandMin && (aod::a3D0Selection::isSelD0 >= selectionFlagD0 || aod::a3D0Selection::isSelD0bar >= selectionFlagD0bar); + Partition> selectedCandidatesMC = aod::a3D0meson::y > -yCandMax&& aod::a3D0meson::y ptCandMin && (aod::a3D0Selection::isSelD0 >= selectionFlagD0 || aod::a3D0Selection::isSelD0bar >= selectionFlagD0bar); + + 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, {{ptDAxisBins, ptDAxisMin, ptDAxisMax}}}}, + {"hPtProng0", "D0,D0bar candidates;prong 0 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{ptDAxisBins, ptDAxisMin, ptDAxisMax}}}}, + {"hPtProng1", "D0,D0bar candidates;prong 1 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{ptDAxisBins, ptDAxisMin, ptDAxisMax}}}}, + {"hSelectionStatus", "D0,D0bar candidates;selection status;entries", {HistType::kTH1F, {{4, -0.5, 3.5}}}}, + {"hEta", "D0,D0bar candidates;candidate #it{#eta};entries", {HistType::kTH1F, {{yAxisBins, yAxisMin, yAxisMax}}}}, + {"hPhi", "D0,D0bar candidates;candidate #it{#varphi};entries", {HistType::kTH1F, {{phiAxisBins, phiAxisMin, phiAxisMax}}}}, + {"hY", "D0,D0bar candidates;candidate #it{y};entries", {HistType::kTH1F, {{yAxisBins, yAxisMin, yAxisMax}}}}, + {"hDDbarVsEtaCut", "D0,D0bar pairs vs #eta cut;#eta_{max};candidates #it{p}_{T} threshold (GeV/#it{c});entries", {HistType::kTH2F, {{static_cast(maxEtaCut / incrementEtaCut), 0., maxEtaCut}, {static_cast(ptThresholdForMaxEtaCut / incrementPtThreshold), 0., ptThresholdForMaxEtaCut}}}}, + {"hPtCandMCRec", "D0,D0bar candidates - MC reco;candidate #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{ptDAxisBins, ptDAxisMin, ptDAxisMax}}}}, + {"hPtProng0MCRec", "D0,D0bar candidates - MC reco;prong 0 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{ptDAxisBins, ptDAxisMin, ptDAxisMax}}}}, + {"hPtProng1MCRec", "D0,D0bar candidates - MC reco;prong 1 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{ptDAxisBins, ptDAxisMin, ptDAxisMax}}}}, + {"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, {{yAxisBins, yAxisMin, yAxisMax}}}}, + {"hPhiMCRec", "D0,D0bar candidates - MC reco;candidate #it{#varphi};entries", {HistType::kTH1F, {{phiAxisBins, phiAxisMin, phiAxisMax}}}}, + {"hYMCRec", "D0,D0bar candidates - MC reco;candidate #it{y};entries", {HistType::kTH1F, {{yAxisBins, yAxisMin, yAxisMax}}}}}}; + + void init(InitContext&) + { + auto vbins = (std::vector)binsPt; + registry.add("hMass", "D0,D0bar candidates;inv. mass (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{massAxisBins, massAxisMin, massAxisMax}, {vbins, "#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, {{massAxisBins, massAxisMin, massAxisMax}, {vbins, "#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, {{massAxisBins, massAxisMin, massAxisMax}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hMassD0MCRecSig", "D0 signal candidates - MC reco;inv. mass (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{massAxisBins, massAxisMin, massAxisMax}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hMassD0MCRecRefl", "D0 reflection candidates - MC reco;inv. mass (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{massAxisBins, massAxisMin, massAxisMax}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hMassD0MCRecBkg", "D0 background candidates - MC reco;inv. mass (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{massAxisBins, massAxisMin, massAxisMax}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hMassD0barMCRecSig", "D0bar signal candidates - MC reco;inv. mass D0bar only (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{massAxisBins, massAxisMin, massAxisMax}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hMassD0barMCRecRefl", "D0bar reflection candidates - MC reco;inv. mass D0bar only (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{massAxisBins, massAxisMin, massAxisMax}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hMassD0barMCRecBkg", "D0bar background candidates - MC reco;inv. mass D0bar only (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{massAxisBins, massAxisMin, massAxisMax}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + } + + /// D0-D0bar correlation pair builder - for real data and data-like analysis (i.e. reco-level w/o matching request via MC truth) + void processData(aod::Collision const& collision, + soa::Join const&) + { + auto selectedCandidatesGrouped = selectedCandidates->sliceByCached(aod::a3D0meson::collisionId, collision.globalIndex(), cache); + + for (const auto& candidate1 : selectedCandidatesGrouped) { // loop over reconstructed and selected D0 and D0bar (together, to fill mass plots first) + double efficiencyWeight = 1.; + if (applyEfficiency) { + efficiencyWeight = 1. / efficiencyD->at(o2::analysis::findBin(binsPt, candidate1.pt())); + } + + // fill invariant mass plots and generic info from all D0/D0bar candidates + if (candidate1.isSelD0() >= selectionFlagD0) { + registry.fill(HIST("hMass"), candidate1.m(), candidate1.pt(), efficiencyWeight); + registry.fill(HIST("hMassD0"), candidate1.m(), candidate1.pt(), efficiencyWeight); + } + if (candidate1.isSelD0bar() >= selectionFlagD0bar) { + registry.fill(HIST("hMass"), candidate1.m(), candidate1.pt(), efficiencyWeight); + registry.fill(HIST("hMassD0bar"), candidate1.m(), candidate1.pt(), efficiencyWeight); + } + 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"), candidate1.y()); + 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() < selectionFlagD0) { + continue; + } + for (const auto& candidate2 : selectedCandidatesGrouped) { + if (candidate2.isSelD0bar() < selectionFlagD0bar) { // keep only D0bar candidates passing the selection + continue; + } + // excluding trigger self-correlations (possible in case of both mass hypotheses accepted) + if (candidate1.mRowIndex == candidate2.mRowIndex) { // this by definition should never happen, since each candidate is either D0 or D0bar + continue; + } + if ((candidate1.pt() - candidate2.pt()) < 1e-5 && (candidate1.eta() - candidate2.eta()) < 1e-5 && (candidate1.phi() - candidate2.phi()) < 1e-5) { // revised, temporary condition to avoid self-correlations (the best would be check the daughterIDs, but we don't store them at the moment) + continue; + } + entryD0D0barPair(getDeltaPhi(candidate2.phi(), candidate1.phi()), + candidate2.eta() - candidate1.eta(), + candidate1.pt(), + candidate2.pt()); + entryD0D0barRecoInfo(candidate1.m(), // mD0 + candidate2.m(), // mD0bar + 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 + } + + PROCESS_SWITCH(alice3correlatorddbar, processData, "Process data", true); + + /// D0-D0bar correlation pair builder - for MC reco-level analysis (candidates matched to true signal only, but also the various bkg sources are studied) + void processMcRec(aod::Collision const& collision, + soa::Join const&) + { + auto selectedCandidatesGroupedMC = selectedCandidatesMC->sliceByCached(aod::a3D0meson::collisionId, collision.globalIndex(), cache); + + // MC reco level + bool flagD0Signal = false; + bool flagD0Reflection = false; + bool flagD0barSignal = false; + bool flagD0barReflection = false; + for (const auto& candidate1 : selectedCandidatesGroupedMC) { + double efficiencyWeight = 1.; + if (applyEfficiency) { + efficiencyWeight = 1. / efficiencyD->at(o2::analysis::findBin(binsPt, candidate1.pt())); + } + + if (candidate1.mcTruthInfo()) { // 1 or 2, i.e. true D0 or D0bar + // fill per-candidate distributions from D0/D0bar true candidates + 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"), candidate1.y()); + registry.fill(HIST("hSelectionStatusMCRec"), candidate1.isSelD0() + (candidate1.isSelD0bar() * 2)); + } + // fill invariant mass plots from D0/D0bar signal and background candidates + if (candidate1.isSelD0() >= selectionFlagD0) { // only reco as D0 + if (candidate1.mcTruthInfo() == 1) { // also matched as D0 + registry.fill(HIST("hMassD0MCRecSig"), candidate1.m(), candidate1.pt(), efficiencyWeight); // here m is univoque, since a given candidate passes the selection with only a single mass option + } else if (candidate1.mcTruthInfo() == 2) { + registry.fill(HIST("hMassD0MCRecRefl"), candidate1.m(), candidate1.pt(), efficiencyWeight); + } else { + registry.fill(HIST("hMassD0MCRecBkg"), candidate1.m(), candidate1.pt(), efficiencyWeight); + } + } + if (candidate1.isSelD0bar() >= selectionFlagD0bar) { // only reco as D0bar + if (candidate1.mcTruthInfo() == 2) { // also matched as D0bar + registry.fill(HIST("hMassD0barMCRecSig"), candidate1.m(), candidate1.pt(), efficiencyWeight); // here m is univoque, since a given candidate passes the selection with only a single mass option + } else if (candidate1.mcTruthInfo() == 1) { + registry.fill(HIST("hMassD0barMCRecRefl"), candidate1.m(), candidate1.pt(), efficiencyWeight); + } else { + registry.fill(HIST("hMassD0barMCRecBkg"), candidate1.m(), candidate1.pt(), efficiencyWeight); + } + } + + // 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.mcTruthInfo() == 1; // flagD0Signal 'true' if candidate1 matched to D0 (particle) + flagD0Reflection = candidate1.mcTruthInfo() == 2; // flagD0Reflection 'true' if candidate1, selected as D0 (particle), is matched to D0bar (antiparticle) + + for (const auto& candidate2 : selectedCandidatesGroupedMC) { + if (candidate2.isSelD0bar() < selectionFlagD0bar) { // discard candidates not selected as D0bar in inner loop + continue; + } + flagD0barSignal = candidate2.mcTruthInfo() == 2; // flagD0barSignal 'true' if candidate2 matched to D0bar (antiparticle) + flagD0barReflection = candidate2.mcTruthInfo() == 1; // flagD0barReflection 'true' if candidate2, selected as D0bar (antiparticle), is matched to D0 (particle) + + // Excluding trigger self-correlations (possible in case of both mass hypotheses of the same real particle, reconstructed as candidate1 for D0 and candidate2 for D0bar) + if (candidate1.mRowIndex == candidate2.mRowIndex) { // this by definition should never happen, since each candidate is either D0 or D0bar + continue; + } + if ((candidate1.pt() - candidate2.pt()) < 1e-5 && (candidate1.eta() - candidate2.eta()) < 1e-5 && (candidate1.phi() - candidate2.phi()) < 1e-5) { // revised, temporary condition to avoid self-correlations (the best would be check the daughterIDs, but we don't store them at the moment) + continue; + } + // choice of options (D0/D0bar signal/bkg) + int pairSignalStatus = 0; // 0 = bkg/bkg, 1 = bkg/ref, 2 = bkg/sig, 3 = ref/bkg, 4 = ref/ref, 5 = ref/sig, 6 = sig/bkg, 7 = sig/ref, 8 = sig/sig + if (flagD0Signal) { + pairSignalStatus += 6; + } + if (flagD0Reflection) { + pairSignalStatus += 3; + } + if (flagD0barSignal) { + pairSignalStatus += 2; + } + if (flagD0barReflection) { + pairSignalStatus += 1; + } + entryD0D0barPair(getDeltaPhi(candidate2.phi(), candidate1.phi()), + candidate2.eta() - candidate1.eta(), + candidate1.pt(), + candidate2.pt()); + entryD0D0barRecoInfo(candidate1.m(), // mD0 + candidate2.m(), // mD0bar + 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 + } + + PROCESS_SWITCH(alice3correlatorddbar, processMcRec, "Process MC Reco mode", false); +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{ + adaptAnalysisTask(cfgc)}; +} diff --git a/ALICE3/TableProducer/alice3-decayfinder.cxx b/ALICE3/TableProducer/alice3-decayfinder.cxx index ae0ea6b5e20..b0aeba5272d 100644 --- a/ALICE3/TableProducer/alice3-decayfinder.cxx +++ b/ALICE3/TableProducer/alice3-decayfinder.cxx @@ -49,6 +49,7 @@ using namespace o2; using namespace o2::framework; +using namespace o2::constants::physics; using namespace o2::framework::expressions; using std::array; @@ -68,6 +69,10 @@ using alice3tracks = soa::Join candidateD0meson; // contains D0 and D0bar selected candidates (separated, i.e. each row with a single mass hypothesis) + Produces selectionOutcome; // flags for isSelD0 and isSelD0bar + Produces mcTruthOutcome; // contains MC truth info (is true D0, true D0bar, or bkg) + // Operation and minimisation criteria Configurable magneticField{"magneticField", 20.0f, "Magnetic field (in kilogauss)"}; Configurable doDCAplotsD{"doDCAplotsD", true, "do daughter prong DCA plots for D mesons"}; @@ -160,23 +165,29 @@ struct alice3decayFinder { float mass; std::array posSV; std::array P; + std::array Pdaug; // positive track + std::array Ndaug; // negative track float pt; + float phi; float eta; + float y; float cosPA; float cosPAxy; float cosThetaStar; float normalizedDecayLength; + int mcTruth; // 0 = bkg, 1 = D0, 2 = D0bar } dmeson; struct { float dcaDau; float mass; float pt; + float phi; float eta; } lcbaryon; template - bool buildDecayCandidateTwoBody(TTrackType const& posTrackRow, TTrackType const& negTrackRow, float posMass, float negMass) + bool buildDecayCandidateTwoBody(TTrackType const& posTrackRow, TTrackType const& negTrackRow, float posMass, float negMass, aod::McParticles const& mcParticles) { o2::track::TrackParCov posTrack = getTrackParCov(posTrackRow); o2::track::TrackParCov negTrack = getTrackParCov(negTrackRow); @@ -201,11 +212,19 @@ struct alice3decayFinder { posTrack.getPxPyPzGlo(posP); negTrack.getPxPyPzGlo(negP); dmeson.dcaDau = TMath::Sqrt(fitter.getChi2AtPCACandidate()); - - // return mass + dmeson.Pdaug[0] = posP[0]; + dmeson.Pdaug[1] = posP[1]; + dmeson.Pdaug[2] = posP[2]; + dmeson.Ndaug[0] = negP[0]; + dmeson.Ndaug[1] = negP[1]; + dmeson.Ndaug[2] = negP[2]; + + // return mass and kinematic variables dmeson.mass = RecoDecay::m(array{array{posP[0], posP[1], posP[2]}, array{negP[0], negP[1], negP[2]}}, array{posMass, negMass}); dmeson.pt = std::hypot(posP[0] + negP[0], posP[1] + negP[1]); + dmeson.phi = RecoDecay::phi(array{posP[0] + negP[0], posP[1] + negP[1]}); dmeson.eta = RecoDecay::eta(array{posP[0] + negP[0], posP[1] + negP[1], posP[2] + negP[2]}); + dmeson.y = RecoDecay::y(std::array{posP[0] + negP[0], posP[1] + negP[1], posP[2] + negP[2]}, dmeson.mass); const auto posSV = fitter.getPCACandidate(); dmeson.posSV[0] = posSV[0]; dmeson.posSV[1] = posSV[1]; @@ -213,6 +232,21 @@ struct alice3decayFinder { o2::track::TrackParCov parentTrack = fitter.createParentTrackParCov(); parentTrack.getPxPyPzGlo(dmeson.P); dmeson.cosThetaStar = RecoDecay::cosThetaStar(std::array{std::array{posP[0], posP[1], posP[2]}, std::array{negP[0], negP[1], negP[2]}}, std::array{posMass, negMass}, dmeson.mass, 0); + + // MC truth check + int indexRec = -1; + int8_t sign = 0; + auto arrayDaughters = std::array{posTrackRow, negTrackRow}; + indexRec = RecoDecay::getMatchedMCRec(mcParticles, arrayDaughters, Pdg::kD0, std::array{+kPiPlus, -kKPlus}, true, &sign); + if (indexRec < 0) { + dmeson.mcTruth = 0; // bkg + } else { + if (sign > 0) { + dmeson.mcTruth = 1; // D0 + } else { + dmeson.mcTruth = 2; // D0bar + } + } return true; } @@ -253,6 +287,7 @@ struct alice3decayFinder { // return mass lcbaryon.mass = RecoDecay::m(array{array{P0[0], P0[1], P0[2]}, array{P1[0], P1[1], P1[2]}, array{P2[0], P2[1], P2[2]}}, array{p0mass, p1mass, p2mass}); lcbaryon.pt = std::hypot(P0[0] + P1[0] + P2[0], P0[1] + P1[1] + P2[1]); + lcbaryon.phi = RecoDecay::phi(array{P0[0] + P1[0] + P2[0], P0[1] + P1[1] + P2[1]}); lcbaryon.eta = RecoDecay::eta(array{P0[0] + P1[0] + P2[0], P0[1] + P1[1] + P2[1], P0[2] + P1[2] + P2[2]}); return true; } @@ -369,7 +404,7 @@ struct alice3decayFinder { } //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* - void processFindDmesons(aod::Collision const& collision, alice3tracks const&, aod::McParticles const&) + void processFindDmesons(aod::Collision const& collision, alice3tracks const&, aod::McParticles const& mcParticles) { // group with this collision auto tracksPiPlusFromDgrouped = tracksPiPlusFromD->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache); @@ -388,12 +423,12 @@ struct alice3decayFinder { histos.fill(HIST("h2dDCAxyVsPtKaMinusFromD"), track.pt(), track.dcaXY() * 1e+4); } - // D mesons + // D0 mesons for (auto const& posTrackRow : tracksPiPlusFromDgrouped) { for (auto const& negTrackRow : tracksKaMinusFromDgrouped) { if (mcSameMotherCheck && !checkSameMother(posTrackRow, negTrackRow)) continue; - if (!buildDecayCandidateTwoBody(posTrackRow, negTrackRow, o2::constants::physics::MassPionCharged, o2::constants::physics::MassKaonCharged)) + if (!buildDecayCandidateTwoBody(posTrackRow, negTrackRow, o2::constants::physics::MassPionCharged, o2::constants::physics::MassKaonCharged, mcParticles)) continue; dmeson.cosPA = RecoDecay::cpa(std::array{collision.posX(), collision.posY(), collision.posZ()}, std::array{dmeson.posSV[0], dmeson.posSV[1], dmeson.posSV[2]}, std::array{dmeson.P[0], dmeson.P[1], dmeson.P[2]}); @@ -433,14 +468,27 @@ struct alice3decayFinder { histos.fill(HIST("hDCADDaughters"), dmeson.dcaDau * 1e+4); histos.fill(HIST("hMassD"), dmeson.mass); histos.fill(HIST("h3dRecD"), dmeson.pt, dmeson.eta, dmeson.mass); + + // store D0 in output table + candidateD0meson(collision.globalIndex(), + dmeson.Pdaug[0], dmeson.Pdaug[1], dmeson.Pdaug[2], + dmeson.Ndaug[0], dmeson.Ndaug[1], dmeson.Ndaug[2], + dmeson.P[0], dmeson.P[1], dmeson.P[2], + dmeson.pt, + dmeson.mass, + dmeson.eta, + dmeson.phi, + dmeson.y); + selectionOutcome(1, 0); // isSelD0 true, isSelD0bar false + mcTruthOutcome(dmeson.mcTruth); } } - // D mesons + // D0bar mesons for (auto const& posTrackRow : tracksKaPlusFromDgrouped) { for (auto const& negTrackRow : tracksPiMinusFromDgrouped) { if (mcSameMotherCheck && !checkSameMother(posTrackRow, negTrackRow)) continue; - if (!buildDecayCandidateTwoBody(posTrackRow, negTrackRow, o2::constants::physics::MassKaonCharged, o2::constants::physics::MassPionCharged)) + if (!buildDecayCandidateTwoBody(posTrackRow, negTrackRow, o2::constants::physics::MassKaonCharged, o2::constants::physics::MassPionCharged, mcParticles)) continue; dmeson.cosPA = RecoDecay::cpa(std::array{collision.posX(), collision.posY(), collision.posZ()}, std::array{dmeson.posSV[0], dmeson.posSV[1], dmeson.posSV[2]}, std::array{dmeson.P[0], dmeson.P[1], dmeson.P[2]}); @@ -480,6 +528,19 @@ struct alice3decayFinder { histos.fill(HIST("hDCADbarDaughters"), dmeson.dcaDau * 1e+4); histos.fill(HIST("hMassDbar"), dmeson.mass); histos.fill(HIST("h3dRecDbar"), dmeson.pt, dmeson.eta, dmeson.mass); + + // store D0bar in output table + candidateD0meson(collision.globalIndex(), + dmeson.Pdaug[0], dmeson.Pdaug[1], dmeson.Pdaug[2], + dmeson.Ndaug[0], dmeson.Ndaug[1], dmeson.Ndaug[2], + dmeson.P[0], dmeson.P[1], dmeson.P[2], + dmeson.pt, + dmeson.mass, + dmeson.eta, + dmeson.phi, + dmeson.y); + selectionOutcome(0, 1); // isSelD0 true, isSelD0bar false + mcTruthOutcome(dmeson.mcTruth); } } } diff --git a/ALICE3/Tasks/CMakeLists.txt b/ALICE3/Tasks/CMakeLists.txt index 36a9d8d58e8..0b2e9972aa5 100644 --- a/ALICE3/Tasks/CMakeLists.txt +++ b/ALICE3/Tasks/CMakeLists.txt @@ -48,3 +48,8 @@ o2physics_add_dpl_workflow(alice3-dilepton SOURCES alice3-dilepton.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::FrameworkPhysicsSupport COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(alice3-taskcorrelationddbar + SOURCES alice3-taskcorrelationDDbar.cxx + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore + COMPONENT_NAME Analysis) \ No newline at end of file diff --git a/ALICE3/Tasks/alice3-taskcorrelationDDbar.cxx b/ALICE3/Tasks/alice3-taskcorrelationDDbar.cxx new file mode 100644 index 00000000000..66f641302ee --- /dev/null +++ b/ALICE3/Tasks/alice3-taskcorrelationDDbar.cxx @@ -0,0 +1,433 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file taskCorrelationDDbar.cxx +/// \brief D0-D0bar correlation analysis task - data-like and MC-reco analysis performance on ALICE 3 detector. +/// +/// \author Fabio Colamaria , INFN Bari + +#include + +#include "Framework/AnalysisTask.h" +#include "Framework/HistogramRegistry.h" +#include "Framework/runDataProcessing.h" + +#include "PWGHF/Utils/utilsAnalysis.h" +#include "PWGHF/HFC/DataModel/CorrelationTables.h" +// #include "PWGHF/DataModel/CandidateReconstructionTables.h" +// #include "PWGHF/DataModel/CandidateSelectionTables.h" + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; + +namespace o2::aod +{ +using DDbarPairFull = soa::Join; +} // namespace o2::aod + +/// +/// 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::PIHalf); +} + +/// +/// 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); +} + +// 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 "; + +// definition of vectors for standard ptbin and invariant mass configurables +const int npTBinsCorrelations = 8; +const double pTBinsCorrelations[npTBinsCorrelations + 1] = {0., 2., 4., 6., 8., 12., 16., 24., 99.}; +auto pTBinsCorrelations_v = std::vector{pTBinsCorrelations, pTBinsCorrelations + npTBinsCorrelations + 1}; +const double signalRegionInnerDefault[npTBinsCorrelations] = {1.810, 1.810, 1.810, 1.810, 1.810, 1.810, 1.810, 1.810}; +const double signalRegionOuterDefault[npTBinsCorrelations] = {1.922, 1.922, 1.922, 1.922, 1.922, 1.922, 1.922, 1.922}; +const double sidebandLeftInnerDefault[npTBinsCorrelations] = {1.642, 1.642, 1.642, 1.642, 1.642, 1.642, 1.642, 1.642}; +const double sidebandLeftOuterDefault[npTBinsCorrelations] = {1.754, 1.754, 1.754, 1.754, 1.754, 1.754, 1.754, 1.754}; +const double sidebandRightInnerDefault[npTBinsCorrelations] = {1.978, 1.978, 1.978, 1.978, 1.978, 1.978, 1.978, 1.978}; +const double sidebandRightOuterDefault[npTBinsCorrelations] = {2.090, 2.090, 2.090, 2.090, 2.090, 2.090, 2.090, 2.090}; +auto signalRegionInner_v = std::vector{signalRegionInnerDefault, signalRegionInnerDefault + npTBinsCorrelations}; +auto signalRegionOuter_v = std::vector{signalRegionOuterDefault, signalRegionOuterDefault + npTBinsCorrelations}; +auto sidebandLeftInner_v = std::vector{sidebandLeftInnerDefault, sidebandLeftInnerDefault + npTBinsCorrelations}; +auto sidebandLeftOuter_v = std::vector{sidebandLeftOuterDefault, sidebandLeftOuterDefault + npTBinsCorrelations}; +auto sidebandRightInner_v = std::vector{sidebandRightInnerDefault, sidebandRightInnerDefault + npTBinsCorrelations}; +auto sidebandRightOuter_v = std::vector{sidebandRightOuterDefault, sidebandRightOuterDefault + npTBinsCorrelations}; +const int npTBinsEfficiency = o2::analysis::hf_cuts_d0_to_pi_k::NBinsPt; +const double efficiencyDmesonDefault[npTBinsEfficiency] = {}; +auto efficiencyDmeson_v = std::vector{efficiencyDmesonDefault, efficiencyDmesonDefault + npTBinsEfficiency}; + +struct alice3taskcorrelationddbar { + Configurable applyEfficiency{"applyEfficiency", 1, "Flag for applying efficiency weights"}; + // pT ranges for correlation plots: the default values are those embedded in hf_cuts_d0_to_pi_k (i.e. the mass pT bins), but can be redefined via json files + Configurable> binsPtCorrelations{"binsPtCorrelations", std::vector{pTBinsCorrelations_v}, "pT bin limits for correlation plots"}; + // pT bins for effiencies: same as above + Configurable> binsPtEfficiency{"binsPtEfficiency", std::vector{o2::analysis::hf_cuts_d0_to_pi_k::vecBinsPt}, "pT bin limits for efficiency"}; + // signal and sideband region edges, to be defined via json file (initialised to empty) + Configurable> signalRegionInner{"signalRegionInner", std::vector{signalRegionInner_v}, "Inner values of signal region vs pT"}; + Configurable> signalRegionOuter{"signalRegionOuter", std::vector{signalRegionOuter_v}, "Outer values of signal region vs pT"}; + Configurable> sidebandLeftInner{"sidebandLeftInner", std::vector{sidebandLeftInner_v}, "Inner values of left sideband vs pT"}; + Configurable> sidebandLeftOuter{"sidebandLeftOuter", std::vector{sidebandLeftOuter_v}, "Outer values of left sideband vs pT"}; + Configurable> sidebandRightInner{"sidebandRightInner", std::vector{sidebandRightInner_v}, "Inner values of right sideband vs pT"}; + Configurable> sidebandRightOuter{"sidebandRightOuter", std::vector{sidebandRightOuter_v}, "Outer values of right sideband vs pT"}; + Configurable> efficiencyD{"efficiencyD", std::vector{efficiencyDmeson_v}, "Efficiency values for D meson specie under study"}; + + 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, {{64, -o2::constants::math::PIHalf, 3. * o2::constants::math::PIHalf}}}}, + {"hCorrel2DPtIntSignalRegion", stringDDbar + stringSignal + stringDeltaPhi + stringDeltaEta + "entries", {HistType::kTH2F, {{64, -o2::constants::math::PIHalf, 3. * o2::constants::math::PIHalf}, {200, -10., 10.}}}}, + {"hCorrel2DVsPtSignalRegion", stringDDbar + stringSignal + stringDeltaPhi + stringDeltaEta + stringPtD + stringPtDbar + "entries", {HistType::kTHnSparseD, {{64, -o2::constants::math::PIHalf, 3. * o2::constants::math::PIHalf}, {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, {{64, -o2::constants::math::PIHalf, 3. * o2::constants::math::PIHalf}}}}, + {"hCorrel2DPtIntSidebands", stringDDbar + stringSideband + stringDeltaPhi + stringDeltaEta + "entries", {HistType::kTH2F, {{64, -o2::constants::math::PIHalf, 3. * o2::constants::math::PIHalf}, {200, -10., 10.}}}}, + {"hCorrel2DVsPtSidebands", stringDDbar + stringSideband + stringDeltaPhi + stringDeltaEta + stringPtD + stringPtDbar + "entries", {HistType::kTHnSparseD, {{64, -o2::constants::math::PIHalf, 3. * o2::constants::math::PIHalf}, {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.}}}}, + {"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() + {"hMass2DCorrelationPairsMCRecBkgRef", stringDDbar + "2D BkgRef - 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() + {"hMass2DCorrelationPairsMCRecRefBkg", stringDDbar + "2D RefBkg - 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() + {"hMass2DCorrelationPairsMCRecRefRef", stringDDbar + "2D RefRef - 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() + {"hMass2DCorrelationPairsMCRecRefSig", stringDDbar + "2D RefSig - 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() + {"hMass2DCorrelationPairsMCRecSigRef", stringDDbar + "2D SigRef - 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() + {"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() + {"hDeltaEtaPtIntSignalRegionMCRec", stringMCReco + stringSignal + stringDeltaEta + "entries", {HistType::kTH1F, {{200, -10., 10.}}}}, + {"hDeltaPhiPtIntSignalRegionMCRec", stringMCReco + stringSignal + stringDeltaPhi + "entries", {HistType::kTH1F, {{64, -o2::constants::math::PIHalf, 3. * o2::constants::math::PIHalf}}}}, + {"hCorrel2DPtIntSignalRegionMCRec", stringMCReco + stringSignal + stringDeltaPhi + stringDeltaEta + "entries", {HistType::kTH2F, {{64, -o2::constants::math::PIHalf, 3. * o2::constants::math::PIHalf}, {200, -10., 10.}}}}, + {"hDeltaPtDDbarSignalRegionMCRec", stringMCReco + stringSignal + stringDeltaPt + "entries", {HistType::kTH1F, {{144, -36., 36.}}}}, + {"hDeltaPtMaxMinSignalRegionMCRec", stringMCReco + stringSignal + stringDeltaPtMaxMin + "entries", {HistType::kTH1F, {{72, 0., 36.}}}}, + {"hCorrel2DVsPtSignalRegionMCRecBkgBkg", stringMCReco + "BkgBkg" + stringSignal + stringDeltaPhi + stringDeltaEta + stringPtD + stringPtDbar + "entries", {HistType::kTHnSparseD, {{64, -o2::constants::math::PIHalf, 3. * o2::constants::math::PIHalf}, {120, -6., 6.}, {10, 0., 10.}, {10, 0., 10.}}}}, // note: axes 3 and 4 (the pT) are updated in the init() + {"hCorrel2DVsPtSignalRegionMCRecBkgRef", stringMCReco + "BkgRef" + stringSignal + stringDeltaPhi + stringDeltaEta + stringPtD + stringPtDbar + "entries", {HistType::kTHnSparseD, {{64, -o2::constants::math::PIHalf, 3. * o2::constants::math::PIHalf}, {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, {{64, -o2::constants::math::PIHalf, 3. * o2::constants::math::PIHalf}, {120, -6., 6.}, {10, 0., 10.}, {10, 0., 10.}}}}, // note: axes 3 and 4 (the pT) are updated in the init() + {"hCorrel2DVsPtSignalRegionMCRecRefBkg", stringMCReco + "RefBkg" + stringSignal + stringDeltaPhi + stringDeltaEta + stringPtD + stringPtDbar + "entries", {HistType::kTHnSparseD, {{64, -o2::constants::math::PIHalf, 3. * o2::constants::math::PIHalf}, {120, -6., 6.}, {10, 0., 10.}, {10, 0., 10.}}}}, // note: axes 3 and 4 (the pT) are updated in the init() + {"hCorrel2DVsPtSignalRegionMCRecRefRef", stringMCReco + "RefRef" + stringSignal + stringDeltaPhi + stringDeltaEta + stringPtD + stringPtDbar + "entries", {HistType::kTHnSparseD, {{64, -o2::constants::math::PIHalf, 3. * o2::constants::math::PIHalf}, {120, -6., 6.}, {10, 0., 10.}, {10, 0., 10.}}}}, // note: axes 3 and 4 (the pT) are updated in the init() + {"hCorrel2DVsPtSignalRegionMCRecRefSig", stringMCReco + "RefSig" + stringSignal + stringDeltaPhi + stringDeltaEta + stringPtD + stringPtDbar + "entries", {HistType::kTHnSparseD, {{64, -o2::constants::math::PIHalf, 3. * o2::constants::math::PIHalf}, {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, {{64, -o2::constants::math::PIHalf, 3. * o2::constants::math::PIHalf}, {120, -6., 6.}, {10, 0., 10.}, {10, 0., 10.}}}}, // note: axes 3 and 4 (the pT) are updated in the init() + {"hCorrel2DVsPtSignalRegionMCRecSigRef", stringMCReco + "SigRef" + stringSignal + stringDeltaPhi + stringDeltaEta + stringPtD + stringPtDbar + "entries", {HistType::kTHnSparseD, {{64, -o2::constants::math::PIHalf, 3. * o2::constants::math::PIHalf}, {120, -6., 6.}, {10, 0., 10.}, {10, 0., 10.}}}}, // note: axes 3 and 4 (the pT) are updated in the init() + {"hCorrel2DVsPtSignalRegionMCRecSigSig", stringMCReco + "SigSig" + stringSignal + stringDeltaPhi + stringDeltaEta + stringPtD + stringPtDbar + "entries", {HistType::kTHnSparseD, {{64, -o2::constants::math::PIHalf, 3. * o2::constants::math::PIHalf}, {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, {{64, -o2::constants::math::PIHalf, 3. * o2::constants::math::PIHalf}}}}, + {"hCorrel2DPtIntSidebandsMCRec", stringMCReco + stringSideband + stringDeltaPhi + stringDeltaEta + "entries", {HistType::kTH2F, {{64, -o2::constants::math::PIHalf, 3. * o2::constants::math::PIHalf}, {200, -10., 10.}}}}, + {"hDeltaPtDDbarSidebandsMCRec", stringMCReco + stringSideband + stringDeltaPt + "entries", {HistType::kTH1F, {{144, -36., 36.}}}}, + {"hDeltaPtMaxMinSidebandsMCRec", stringMCReco + stringSideband + stringDeltaPtMaxMin + "entries", {HistType::kTH1F, {{72, 0., 36.}}}}, + {"hCorrel2DVsPtSidebandsMCRecBkgBkg", stringMCReco + "BkgBkg" + stringSideband + stringDeltaPhi + stringDeltaEta + stringPtD + stringPtDbar + "entries", {HistType::kTHnSparseD, {{64, -o2::constants::math::PIHalf, 3. * o2::constants::math::PIHalf}, {120, -6., 6.}, {10, 0., 10.}, {10, 0., 10.}}}}, // note: axes 3 and 4 (the pT) are updated in the init() + {"hCorrel2DVsPtSidebandsMCRecBkgRef", stringMCReco + "BkgRef" + stringSideband + stringDeltaPhi + stringDeltaEta + stringPtD + stringPtDbar + "entries", {HistType::kTHnSparseD, {{64, -o2::constants::math::PIHalf, 3. * o2::constants::math::PIHalf}, {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, {{64, -o2::constants::math::PIHalf, 3. * o2::constants::math::PIHalf}, {120, -6., 6.}, {10, 0., 10.}, {10, 0., 10.}}}}, // note: axes 3 and 4 (the pT) are updated in the init() + {"hCorrel2DVsPtSidebandsMCRecRefBkg", stringMCReco + "RefBkg" + stringSideband + stringDeltaPhi + stringDeltaEta + stringPtD + stringPtDbar + "entries", {HistType::kTHnSparseD, {{64, -o2::constants::math::PIHalf, 3. * o2::constants::math::PIHalf}, {120, -6., 6.}, {10, 0., 10.}, {10, 0., 10.}}}}, // note: axes 3 and 4 (the pT) are updated in the init() + {"hCorrel2DVsPtSidebandsMCRecRefRef", stringMCReco + "RefRef" + stringSideband + stringDeltaPhi + stringDeltaEta + stringPtD + stringPtDbar + "entries", {HistType::kTHnSparseD, {{64, -o2::constants::math::PIHalf, 3. * o2::constants::math::PIHalf}, {120, -6., 6.}, {10, 0., 10.}, {10, 0., 10.}}}}, // note: axes 3 and 4 (the pT) are updated in the init() + {"hCorrel2DVsPtSidebandsMCRecRefSig", stringMCReco + "RefSig" + stringSideband + stringDeltaPhi + stringDeltaEta + stringPtD + stringPtDbar + "entries", {HistType::kTHnSparseD, {{64, -o2::constants::math::PIHalf, 3. * o2::constants::math::PIHalf}, {120, -6., 6.}, {10, 0., 10.}, {10, 0., 10.}}}}, // note: axes 3 and 4 (the pT) are updated in the init() + {"hCorrel2DVsPtSidebandsMCRecSigBkg", stringMCReco + "SigBkg" + stringSideband + stringDeltaPhi + stringDeltaEta + stringPtD + stringPtDbar + "entries", {HistType::kTHnSparseD, {{64, -o2::constants::math::PIHalf, 3. * o2::constants::math::PIHalf}, {120, -6., 6.}, {10, 0., 10.}, {10, 0., 10.}}}}, // note: axes 3 and 4 (the pT) are updated in the init() + {"hCorrel2DVsPtSidebandsMCRecSigRef", stringMCReco + "SigRef" + stringSideband + stringDeltaPhi + stringDeltaEta + stringPtD + stringPtDbar + "entries", {HistType::kTHnSparseD, {{64, -o2::constants::math::PIHalf, 3. * o2::constants::math::PIHalf}, {120, -6., 6.}, {10, 0., 10.}, {10, 0., 10.}}}}, // note: axes 3 and 4 (the pT) are updated in the init() + {"hCorrel2DVsPtSidebandsMCRecSigSig", stringMCReco + "SigSig" + stringSideband + stringDeltaPhi + stringDeltaEta + stringPtD + stringPtDbar + "entries", {HistType::kTHnSparseD, {{64, -o2::constants::math::PIHalf, 3. * o2::constants::math::PIHalf}, {120, -6., 6.}, {10, 0., 10.}, {10, 0., 10.}}}}}}; // note: axes 3 and 4 (the pT) are updated in the init() + + void init(InitContext&) + { + // redefinition of pT axes for THnSparse holding correlation entries + int nBinspTaxis = binsPtCorrelations->size() - 1; + const double* valuespTaxis = binsPtCorrelations->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); + registry.get(HIST("hMass2DCorrelationPairs"))->Sumw2(); + registry.get(HIST("hCorrel2DVsPtSignalRegion"))->Sumw2(); + registry.get(HIST("hCorrel2DVsPtSidebands"))->Sumw2(); + registry.get(HIST("hMass2DCorrelationPairsMCRecBkgBkg"))->GetAxis(i)->Set(nBinspTaxis, valuespTaxis); + registry.get(HIST("hMass2DCorrelationPairsMCRecBkgRef"))->GetAxis(i)->Set(nBinspTaxis, valuespTaxis); + registry.get(HIST("hMass2DCorrelationPairsMCRecBkgSig"))->GetAxis(i)->Set(nBinspTaxis, valuespTaxis); + registry.get(HIST("hMass2DCorrelationPairsMCRecRefBkg"))->GetAxis(i)->Set(nBinspTaxis, valuespTaxis); + registry.get(HIST("hMass2DCorrelationPairsMCRecRefRef"))->GetAxis(i)->Set(nBinspTaxis, valuespTaxis); + registry.get(HIST("hMass2DCorrelationPairsMCRecRefSig"))->GetAxis(i)->Set(nBinspTaxis, valuespTaxis); + registry.get(HIST("hMass2DCorrelationPairsMCRecSigBkg"))->GetAxis(i)->Set(nBinspTaxis, valuespTaxis); + registry.get(HIST("hMass2DCorrelationPairsMCRecSigRef"))->GetAxis(i)->Set(nBinspTaxis, valuespTaxis); + registry.get(HIST("hMass2DCorrelationPairsMCRecSigSig"))->GetAxis(i)->Set(nBinspTaxis, valuespTaxis); + registry.get(HIST("hCorrel2DVsPtSignalRegionMCRecBkgBkg"))->GetAxis(i)->Set(nBinspTaxis, valuespTaxis); + registry.get(HIST("hCorrel2DVsPtSignalRegionMCRecBkgRef"))->GetAxis(i)->Set(nBinspTaxis, valuespTaxis); + registry.get(HIST("hCorrel2DVsPtSignalRegionMCRecBkgSig"))->GetAxis(i)->Set(nBinspTaxis, valuespTaxis); + registry.get(HIST("hCorrel2DVsPtSignalRegionMCRecRefBkg"))->GetAxis(i)->Set(nBinspTaxis, valuespTaxis); + registry.get(HIST("hCorrel2DVsPtSignalRegionMCRecRefRef"))->GetAxis(i)->Set(nBinspTaxis, valuespTaxis); + registry.get(HIST("hCorrel2DVsPtSignalRegionMCRecRefSig"))->GetAxis(i)->Set(nBinspTaxis, valuespTaxis); + registry.get(HIST("hCorrel2DVsPtSignalRegionMCRecSigBkg"))->GetAxis(i)->Set(nBinspTaxis, valuespTaxis); + registry.get(HIST("hCorrel2DVsPtSignalRegionMCRecSigRef"))->GetAxis(i)->Set(nBinspTaxis, valuespTaxis); + registry.get(HIST("hCorrel2DVsPtSignalRegionMCRecSigSig"))->GetAxis(i)->Set(nBinspTaxis, valuespTaxis); + registry.get(HIST("hCorrel2DVsPtSidebandsMCRecBkgBkg"))->GetAxis(i)->Set(nBinspTaxis, valuespTaxis); + registry.get(HIST("hCorrel2DVsPtSidebandsMCRecBkgRef"))->GetAxis(i)->Set(nBinspTaxis, valuespTaxis); + registry.get(HIST("hCorrel2DVsPtSidebandsMCRecBkgSig"))->GetAxis(i)->Set(nBinspTaxis, valuespTaxis); + registry.get(HIST("hCorrel2DVsPtSidebandsMCRecRefBkg"))->GetAxis(i)->Set(nBinspTaxis, valuespTaxis); + registry.get(HIST("hCorrel2DVsPtSidebandsMCRecRefRef"))->GetAxis(i)->Set(nBinspTaxis, valuespTaxis); + registry.get(HIST("hCorrel2DVsPtSidebandsMCRecRefSig"))->GetAxis(i)->Set(nBinspTaxis, valuespTaxis); + registry.get(HIST("hCorrel2DVsPtSidebandsMCRecSigBkg"))->GetAxis(i)->Set(nBinspTaxis, valuespTaxis); + registry.get(HIST("hCorrel2DVsPtSidebandsMCRecSigRef"))->GetAxis(i)->Set(nBinspTaxis, valuespTaxis); + registry.get(HIST("hCorrel2DVsPtSidebandsMCRecSigSig"))->GetAxis(i)->Set(nBinspTaxis, valuespTaxis); + registry.get(HIST("hMass2DCorrelationPairsMCRecBkgBkg"))->Sumw2(); + registry.get(HIST("hMass2DCorrelationPairsMCRecBkgRef"))->Sumw2(); + registry.get(HIST("hMass2DCorrelationPairsMCRecBkgSig"))->Sumw2(); + registry.get(HIST("hMass2DCorrelationPairsMCRecRefBkg"))->Sumw2(); + registry.get(HIST("hMass2DCorrelationPairsMCRecRefRef"))->Sumw2(); + registry.get(HIST("hMass2DCorrelationPairsMCRecRefSig"))->Sumw2(); + registry.get(HIST("hMass2DCorrelationPairsMCRecSigBkg"))->Sumw2(); + registry.get(HIST("hMass2DCorrelationPairsMCRecSigRef"))->Sumw2(); + registry.get(HIST("hMass2DCorrelationPairsMCRecSigSig"))->Sumw2(); + registry.get(HIST("hCorrel2DVsPtSignalRegionMCRecBkgBkg"))->Sumw2(); + registry.get(HIST("hCorrel2DVsPtSignalRegionMCRecBkgRef"))->Sumw2(); + registry.get(HIST("hCorrel2DVsPtSignalRegionMCRecBkgSig"))->Sumw2(); + registry.get(HIST("hCorrel2DVsPtSignalRegionMCRecRefBkg"))->Sumw2(); + registry.get(HIST("hCorrel2DVsPtSignalRegionMCRecRefRef"))->Sumw2(); + registry.get(HIST("hCorrel2DVsPtSignalRegionMCRecRefSig"))->Sumw2(); + registry.get(HIST("hCorrel2DVsPtSignalRegionMCRecSigBkg"))->Sumw2(); + registry.get(HIST("hCorrel2DVsPtSignalRegionMCRecSigRef"))->Sumw2(); + registry.get(HIST("hCorrel2DVsPtSignalRegionMCRecSigSig"))->Sumw2(); + registry.get(HIST("hCorrel2DVsPtSidebandsMCRecBkgBkg"))->Sumw2(); + registry.get(HIST("hCorrel2DVsPtSidebandsMCRecBkgRef"))->Sumw2(); + registry.get(HIST("hCorrel2DVsPtSidebandsMCRecBkgSig"))->Sumw2(); + registry.get(HIST("hCorrel2DVsPtSidebandsMCRecRefBkg"))->Sumw2(); + registry.get(HIST("hCorrel2DVsPtSidebandsMCRecRefRef"))->Sumw2(); + registry.get(HIST("hCorrel2DVsPtSidebandsMCRecRefSig"))->Sumw2(); + registry.get(HIST("hCorrel2DVsPtSidebandsMCRecSigBkg"))->Sumw2(); + registry.get(HIST("hCorrel2DVsPtSidebandsMCRecSigRef"))->Sumw2(); + registry.get(HIST("hCorrel2DVsPtSidebandsMCRecSigSig"))->Sumw2(); + } + } + + /// 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 + void processData(aod::DDbarPairFull const& pairEntries) + { + for (const 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(); + + int pTBinD = o2::analysis::findBin(binsPtCorrelations, ptD); + int pTBinDbar = o2::analysis::findBin(binsPtCorrelations, ptDbar); + + double efficiencyWeight = 1.; + if (applyEfficiency) { + efficiencyWeight = 1. / (efficiencyD->at(o2::analysis::findBin(binsPtEfficiency, ptD)) * efficiencyD->at(o2::analysis::findBin(binsPtEfficiency, ptDbar))); + } + + // fill 2D invariant mass plots + registry.fill(HIST("hMass2DCorrelationPairs"), massD, massDbar, ptD, ptDbar, efficiencyWeight); + + // reject entries outside pT ranges of interest + 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, efficiencyWeight); + registry.fill(HIST("hCorrel2DPtIntSignalRegion"), deltaPhi, deltaEta, efficiencyWeight); + registry.fill(HIST("hDeltaEtaPtIntSignalRegion"), deltaEta, efficiencyWeight); + registry.fill(HIST("hDeltaPhiPtIntSignalRegion"), deltaPhi, efficiencyWeight); + registry.fill(HIST("hDeltaPtDDbarSignalRegion"), ptDbar - ptD, efficiencyWeight); + registry.fill(HIST("hDeltaPtMaxMinSignalRegion"), std::abs(ptDbar - ptD), efficiencyWeight); + } + + 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, efficiencyWeight); + registry.fill(HIST("hCorrel2DPtIntSidebands"), deltaPhi, deltaEta, efficiencyWeight); + registry.fill(HIST("hDeltaEtaPtIntSidebands"), deltaEta, efficiencyWeight); + registry.fill(HIST("hDeltaPhiPtIntSidebands"), deltaPhi, efficiencyWeight); + registry.fill(HIST("hDeltaPtDDbarSidebands"), ptDbar - ptD, efficiencyWeight); + registry.fill(HIST("hDeltaPtMaxMinSidebands"), std::abs(ptDbar - ptD), efficiencyWeight); + } + } // end loop + } + + PROCESS_SWITCH(alice3taskcorrelationddbar, processData, "Process data", true); + + /// 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 + void processMcRec(aod::DDbarPairFull const& pairEntries) + { + for (const 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(); + + int pTBinD = o2::analysis::findBin(binsPtCorrelations, ptD); + int pTBinDbar = o2::analysis::findBin(binsPtCorrelations, ptDbar); + + double efficiencyWeight = 1.; + if (applyEfficiency) { + efficiencyWeight = 1. / (efficiencyD->at(o2::analysis::findBin(binsPtEfficiency, ptD)) * efficiencyD->at(o2::analysis::findBin(binsPtEfficiency, ptDbar))); + } + + // fill 2D invariant mass plots + switch (pairEntry.signalStatus()) { + case 0: // D Bkg, Dbar Bkg + registry.fill(HIST("hMass2DCorrelationPairsMCRecBkgBkg"), massD, massDbar, ptD, ptDbar, efficiencyWeight); + break; + case 1: // D Bkg, Dbar Ref + registry.fill(HIST("hMass2DCorrelationPairsMCRecBkgRef"), massD, massDbar, ptD, ptDbar, efficiencyWeight); + break; + case 2: // D Bkg, Dbar Sig + registry.fill(HIST("hMass2DCorrelationPairsMCRecBkgSig"), massD, massDbar, ptD, ptDbar, efficiencyWeight); + break; + case 3: // D Ref, Dbar Bkg + registry.fill(HIST("hMass2DCorrelationPairsMCRecRefBkg"), massD, massDbar, ptD, ptDbar, efficiencyWeight); + break; + case 4: // D Ref, Dbar Ref + registry.fill(HIST("hMass2DCorrelationPairsMCRecRefRef"), massD, massDbar, ptD, ptDbar, efficiencyWeight); + break; + case 5: // D Ref, Dbar Sig + registry.fill(HIST("hMass2DCorrelationPairsMCRecRefSig"), massD, massDbar, ptD, ptDbar, efficiencyWeight); + break; + case 6: // D Sig, Dbar Bkg + registry.fill(HIST("hMass2DCorrelationPairsMCRecSigBkg"), massD, massDbar, ptD, ptDbar, efficiencyWeight); + break; + case 7: // D Sig, Dbar Ref + registry.fill(HIST("hMass2DCorrelationPairsMCRecSigRef"), massD, massDbar, ptD, ptDbar, efficiencyWeight); + break; + case 8: // D Sig, Dbar Sig + registry.fill(HIST("hMass2DCorrelationPairsMCRecSigSig"), massD, massDbar, ptD, ptDbar, efficiencyWeight); + break; + default: // should not happen for MC reco + break; + } + + // reject entries outside pT ranges of interest + 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("hCorrel2DPtIntSignalRegionMCRec"), deltaPhi, deltaEta, efficiencyWeight); + registry.fill(HIST("hDeltaEtaPtIntSignalRegionMCRec"), deltaEta, efficiencyWeight); + registry.fill(HIST("hDeltaPhiPtIntSignalRegionMCRec"), deltaPhi, efficiencyWeight); + registry.fill(HIST("hDeltaPtDDbarSignalRegionMCRec"), ptDbar - ptD, efficiencyWeight); + registry.fill(HIST("hDeltaPtMaxMinSignalRegionMCRec"), std::abs(ptDbar - ptD), efficiencyWeight); + switch (pairEntry.signalStatus()) { + case 0: // D Bkg, Dbar Bkg + registry.fill(HIST("hCorrel2DVsPtSignalRegionMCRecBkgBkg"), deltaPhi, deltaEta, ptD, ptDbar, efficiencyWeight); + break; + case 1: // D Bkg, Dbar Ref + registry.fill(HIST("hCorrel2DVsPtSignalRegionMCRecBkgRef"), deltaPhi, deltaEta, ptD, ptDbar, efficiencyWeight); + break; + case 2: // D Bkg, Dbar Sig + registry.fill(HIST("hCorrel2DVsPtSignalRegionMCRecBkgSig"), deltaPhi, deltaEta, ptD, ptDbar, efficiencyWeight); + break; + case 3: // D Ref, Dbar Bkg + registry.fill(HIST("hCorrel2DVsPtSignalRegionMCRecRefBkg"), deltaPhi, deltaEta, ptD, ptDbar, efficiencyWeight); + break; + case 4: // D Ref, Dbar Ref + registry.fill(HIST("hCorrel2DVsPtSignalRegionMCRecRefRef"), deltaPhi, deltaEta, ptD, ptDbar, efficiencyWeight); + break; + case 5: // D Ref, Dbar Sig + registry.fill(HIST("hCorrel2DVsPtSignalRegionMCRecRefSig"), deltaPhi, deltaEta, ptD, ptDbar, efficiencyWeight); + break; + case 6: // D Sig, Dbar Bkg + registry.fill(HIST("hCorrel2DVsPtSignalRegionMCRecSigBkg"), deltaPhi, deltaEta, ptD, ptDbar, efficiencyWeight); + break; + case 7: // D Sig, Dbar Ref + registry.fill(HIST("hCorrel2DVsPtSignalRegionMCRecSigRef"), deltaPhi, deltaEta, ptD, ptDbar, efficiencyWeight); + break; + case 8: // D Sig, Dbar Sig + registry.fill(HIST("hCorrel2DVsPtSignalRegionMCRecSigSig"), deltaPhi, deltaEta, ptD, ptDbar, efficiencyWeight); + 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, efficiencyWeight); + registry.fill(HIST("hDeltaEtaPtIntSidebandsMCRec"), deltaEta, efficiencyWeight); + registry.fill(HIST("hDeltaPhiPtIntSidebandsMCRec"), deltaPhi, efficiencyWeight); + registry.fill(HIST("hDeltaPtDDbarSidebandsMCRec"), ptDbar - ptD, efficiencyWeight); + registry.fill(HIST("hDeltaPtMaxMinSidebandsMCRec"), std::abs(ptDbar - ptD), efficiencyWeight); + switch (pairEntry.signalStatus()) { + case 0: // D Bkg, Dbar Bkg + registry.fill(HIST("hCorrel2DVsPtSidebandsMCRecBkgBkg"), deltaPhi, deltaEta, ptD, ptDbar, efficiencyWeight); + break; + case 1: // D Bkg, Dbar Ref + registry.fill(HIST("hCorrel2DVsPtSidebandsMCRecBkgRef"), deltaPhi, deltaEta, ptD, ptDbar, efficiencyWeight); + break; + case 2: // D Bkg, Dbar Sig + registry.fill(HIST("hCorrel2DVsPtSidebandsMCRecBkgSig"), deltaPhi, deltaEta, ptD, ptDbar, efficiencyWeight); + break; + case 3: // D Ref, Dbar Bkg + registry.fill(HIST("hCorrel2DVsPtSidebandsMCRecRefBkg"), deltaPhi, deltaEta, ptD, ptDbar, efficiencyWeight); + break; + case 4: // D Ref, Dbar Ref + registry.fill(HIST("hCorrel2DVsPtSidebandsMCRecRefRef"), deltaPhi, deltaEta, ptD, ptDbar, efficiencyWeight); + break; + case 5: // D Ref, Dbar Sig + registry.fill(HIST("hCorrel2DVsPtSidebandsMCRecRefSig"), deltaPhi, deltaEta, ptD, ptDbar, efficiencyWeight); + break; + case 6: // D Sig, Dbar Bkg + registry.fill(HIST("hCorrel2DVsPtSidebandsMCRecSigBkg"), deltaPhi, deltaEta, ptD, ptDbar, efficiencyWeight); + break; + case 7: // D Sig, Dbar Ref + registry.fill(HIST("hCorrel2DVsPtSidebandsMCRecSigRef"), deltaPhi, deltaEta, ptD, ptDbar, efficiencyWeight); + break; + case 8: // D Sig, Dbar Sig + registry.fill(HIST("hCorrel2DVsPtSidebandsMCRecSigSig"), deltaPhi, deltaEta, ptD, ptDbar, efficiencyWeight); + break; + default: // should not happen for MC reco + break; + } + } + } // end loop + } + + PROCESS_SWITCH(alice3taskcorrelationddbar, processMcRec, "Process MC Reco mode", false); +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{adaptAnalysisTask(cfgc)}; +}