diff --git a/PWGHF/Tasks/HFCorrelatorD0D0bar.cxx b/PWGHF/Tasks/HFCorrelatorD0D0bar.cxx index cc4f24e543f..78ac9815211 100644 --- a/PWGHF/Tasks/HFCorrelatorD0D0bar.cxx +++ b/PWGHF/Tasks/HFCorrelatorD0D0bar.cxx @@ -29,18 +29,6 @@ using namespace o2::aod::hf_correlation_ddbar; using namespace o2::analysis::hf_cuts_d0_topik; using namespace o2::constants::math; -void customize(std::vector& workflowOptions) -{ - ConfigParamSpec optionDoLikeSign{"doLikeSign", VariantType::Bool, false, {"Run Like-Sign analysis."}}; - ConfigParamSpec optionDoMCccbar{"doMCccbar", VariantType::Bool, false, {"Run MC-Gen dedicated tasks."}}; - ConfigParamSpec optionDoMCGen{"doMCGen", VariantType::Bool, false, {"Run MC-Gen dedicated tasks."}}; - ConfigParamSpec optionDoMCRec{"doMCRec", VariantType::Bool, true, {"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" /// @@ -62,26 +50,55 @@ const int npTBinsMassAndEfficiency = o2::analysis::hf_cuts_d0_topik::npTBins; 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 = 2. * o2::constants::math::PI; +const int yAxisBins = 100; +const double yAxisMin = -5.; +const double yAxisMax = 5.; +const int ptDAxisBins = 180; +const double ptDAxisMin = 0.; +const double ptDAxisMax = 36.; + using MCParticlesPlus = soa::Join; -/// D0-D0bar correlation pair builder - for real data and data-like analysis (i.e. reco-level w/o matching request via MC truth) struct HfCorrelatorD0D0bar { Produces entryD0D0barPair; Produces entryD0D0barRecoInfo; HistogramRegistry registry{ "registry", - //NOTE: use hMassD0 for trigger normalisation (S*0.955), and hMass2DCorrelationPairs (in final task) for 2D-sideband-subtraction purposes - {{"hPtCand", "D0,D0bar candidates;candidate #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, - {"hPtProng0", "D0,D0bar candidates;prong 0 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, - {"hPtProng1", "D0,D0bar candidates;prong 1 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, + // 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, {{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.}}}}, + {"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}}}}, {"hMultiplicityPreSelection", "multiplicity prior to selection;multiplicity;entries", {HistType::kTH1F, {{10000, 0., 10000.}}}}, {"hMultiplicity", "multiplicity;multiplicity;entries", {HistType::kTH1F, {{10000, 0., 10000.}}}}, - {"hDDbarVsEtaCut", "D0,D0bar pairs vs #eta cut;#eta_{max};entries", {HistType::kTH2F, {{(int)(maxEtaCut / incrementEtaCut), 0., maxEtaCut}, {(int)(ptThresholdForMaxEtaCut / incrementPtThreshold), 0., ptThresholdForMaxEtaCut}}}}}}; + {"hDDbarVsEtaCut", "D0,D0bar pairs vs #eta cut;#eta_{max};candidates #it{p}_{T} threshold (GeV/#it{c});entries", {HistType::kTH2F, {{(int)(maxEtaCut / incrementEtaCut), 0., maxEtaCut}, {(int)(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}}}}, + {"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, {{ptDAxisBins, ptDAxisMin, ptDAxisMax}}}}, + {"hEtaMCGen", "D0,D0bar particles - MC gen;particle #it{#eta};entries", {HistType::kTH1F, {{yAxisBins, yAxisMin, yAxisMax}}}}, + {"hPhiMCGen", "D0,D0bar particles - MC gen;particle #it{#varphi};entries", {HistType::kTH1F, {{phiAxisBins, phiAxisMin, phiAxisMax}}}}, + {"hYMCGen", "D0,D0bar candidates - MC gen;candidate #it{y};entries", {HistType::kTH1F, {{yAxisBins, yAxisMin, yAxisMax}}}}, + {"hCountD0D0barPerEvent", "D0,D0bar particles - MC gen;Number per event;entries", {HistType::kTH1F, {{20, 0., 20.}}}}, + {"hDDbarVsDaughterEtaCut", "D0,D0bar pairs vs #eta cut on D daughters;#eta_{max};candidates #it{p}_{T} threshold (GeV/#it{c});entries", {HistType::kTH2F, {{(int)(maxEtaCut / incrementEtaCut), 0., maxEtaCut}, {(int)(ptThresholdForMaxEtaCut / incrementPtThreshold), 0., ptThresholdForMaxEtaCut}}}}, + {"hCountCCbarPerEvent", "c,cbar particles - MC gen;Number per event;entries", {HistType::kTH1F, {{20, 0., 20.}}}}, + {"hCountCCbarPerEventBeforeEtaCut", "c,cbar particles - MC gen;Number per event pre #eta cut;entries", {HistType::kTH1F, {{20, 0., 20.}}}}}}; Configurable selectionFlagD0{"selectionFlagD0", 1, "Selection Flag for D0"}; Configurable selectionFlagD0bar{"selectionFlagD0bar", 1, "Selection Flag for D0bar"}; @@ -95,14 +112,24 @@ struct HfCorrelatorD0D0bar { 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})"}}}); + auto vbins = (std::vector)bins; + 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})"}}}); + registry.add("hCountD0triggersMCGen", "D0 trigger particles - MC gen;;N of trigger D0", {HistType::kTH2F, {{1, -0.5, 0.5}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hCountCtriggersMCGen", "c trigger particles - MC gen;;N of trigger c quark", {HistType::kTH2F, {{1, -0.5, 0.5}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); } Filter filterSelectCandidates = (aod::hf_selcandidate_d0::isSelD0 >= selectionFlagD0 || aod::hf_selcandidate_d0::isSelD0bar >= selectionFlagD0bar); - void process(aod::Collision const& collision, soa::Join& tracks, soa::Filtered> const& candidates) + /// 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& tracks, soa::Filtered> const& candidates) { int nTracks = 0; if (collision.numContrib() > 1) { @@ -129,7 +156,7 @@ struct HfCorrelatorD0D0bar { if (cutPtCandMin >= 0. && candidate1.pt() < cutPtCandMin) { continue; } - //check decay channel flag for candidate1 + // check decay channel flag for candidate1 if (!(candidate1.hfflag() & 1 << DecayType::D0ToPiK)) { continue; } @@ -139,7 +166,7 @@ struct HfCorrelatorD0D0bar { efficiencyWeight = 1. / efficiencyDmeson->at(o2::analysis::findBin(bins, candidate1.pt())); } - //fill invariant mass plots and generic info from all D0/D0bar candidates + // fill invariant mass plots and generic info from all D0/D0bar candidates if (candidate1.isSelD0() >= selectionFlagD0) { registry.fill(HIST("hMass"), InvMassD0(candidate1), candidate1.pt(), efficiencyWeight); registry.fill(HIST("hMassD0"), InvMassD0(candidate1), candidate1.pt(), efficiencyWeight); @@ -156,26 +183,26 @@ struct HfCorrelatorD0D0bar { registry.fill(HIST("hY"), YD0(candidate1)); registry.fill(HIST("hSelectionStatus"), candidate1.isSelD0bar() + (candidate1.isSelD0() * 2)); - //D-Dbar correlation dedicated section - //if the candidate is a D0, search for D0bar and evaluate correlations + // D-Dbar correlation dedicated section + // if the candidate is a D0, search for D0bar and evaluate correlations if (candidate1.isSelD0() < selectionFlagD0) { continue; } for (auto& candidate2 : candidates) { - if (!(candidate2.hfflag() & 1 << DecayType::D0ToPiK)) { //check decay channel flag for candidate2 + if (!(candidate2.hfflag() & 1 << DecayType::D0ToPiK)) { // check decay channel flag for candidate2 continue; } - if (candidate2.isSelD0bar() < selectionFlagD0bar) { //keep only D0bar candidates passing the selection + if (candidate2.isSelD0bar() < selectionFlagD0bar) { // keep only D0bar candidates passing the selection continue; } - //kinematic selection on D0bar candidates + // kinematic selection on D0bar candidates if (cutYCandMax >= 0. && std::abs(YD0(candidate2)) > cutYCandMax) { continue; } if (cutPtCandMin >= 0. && candidate2.pt() < cutPtCandMin) { continue; } - //excluding trigger self-correlations (possible in case of both mass hypotheses accepted) + // excluding trigger self-correlations (possible in case of both mass hypotheses accepted) if (candidate1.mRowIndex == candidate2.mRowIndex) { continue; } @@ -188,68 +215,28 @@ struct HfCorrelatorD0D0bar { 0); double etaCut = 0.; double ptCut = 0.; - do { //fill pairs vs etaCut plot + do { // fill pairs vs etaCut plot ptCut = 0.; etaCut += incrementEtaCut; - do { //fill pairs vs etaCut plot + 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) + // note: candidates selected as both D0 and D0bar are used, and considered in both situation (but not auto-correlated): reflections could play a relevant role. + // another, more restrictive, option, could be to consider only candidates selected with a single option (D0 xor D0bar) } // end inner loop (Dbars) - } //end outer loop - } -}; - -/// D0-D0bar correlation pair builder - for MC reco-level analysis (candidates matched to true signal only, but also the various bkg sources are studied) -struct HfCorrelatorD0D0barMcRec { - - Produces entryD0D0barPair; - Produces entryD0D0barRecoInfo; - - HistogramRegistry registry{ - "registry", - //NOTE: use hMassD0 for trigger normalisation (S*0.955), and hMass2DCorrelationPairs (in final task) for 2D-sideband-subtraction purposes - {{"hPtCandMCRec", "D0,D0bar candidates - MC reco;candidate #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, - {"hPtProng0MCRec", "D0,D0bar candidates - MC reco;prong 0 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, - {"hPtProng1MCRec", "D0,D0bar candidates - MC reco;prong 1 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, - {"hSelectionStatusMCRec", "D0,D0bar candidates - MC reco;selection status;entries", {HistType::kTH1F, {{4, -0.5, 3.5}}}}, - {"hEtaMCRec", "D0,D0bar candidates - MC reco;candidate #it{#eta};entries", {HistType::kTH1F, {{100, -5., 5.}}}}, - {"hPhiMCRec", "D0,D0bar candidates - MC reco;candidate #it{#varphi};entries", {HistType::kTH1F, {{32, 0., 2. * o2::constants::math::PI}}}}, - {"hYMCRec", "D0,D0bar candidates - MC reco;candidate #it{#y};entries", {HistType::kTH1F, {{100, -5., 5.}}}}, - {"hMultiplicityPreSelection", "multiplicity prior to selection;multiplicity;entries", {HistType::kTH1F, {{10000, 0., 10000.}}}}, - {"hMultiplicity", "multiplicity;multiplicity;entries", {HistType::kTH1F, {{10000, 0., 10000.}}}}, - {"hDDbarVsEtaCut", "D0,D0bar pairs vs #eta cut;#eta_{max};entries", {HistType::kTH2F, {{(int)(maxEtaCut / incrementEtaCut), 0., maxEtaCut}, {(int)(ptThresholdForMaxEtaCut / incrementPtThreshold), 0., ptThresholdForMaxEtaCut}}}}}}; - - Configurable selectionFlagD0{"selectionFlagD0", 1, "Selection Flag for D0"}; - Configurable selectionFlagD0bar{"selectionFlagD0bar", 1, "Selection Flag for D0bar"}; - Configurable cutYCandMax{"cutYCandMax", -1., "max. cand. rapidity"}; - Configurable cutPtCandMin{"cutPtCandMin", -1., "min. cand. pT"}; - Configurable> bins{"ptBinsForMassAndEfficiency", std::vector{o2::analysis::hf_cuts_d0_topik::pTBins_v}, "pT bin limits for candidate mass plots and efficiency"}; - Configurable> efficiencyDmeson{"efficiencyDmeson", std::vector{efficiencyDmeson_v}, "Efficiency values for D0 meson"}; - Configurable flagApplyEfficiency{"efficiencyFlagD", 1, "Flag for applying D-meson efficiency weights"}; - Configurable multMin{"multMin", 0., "minimum multiplicity accepted"}; - Configurable multMax{"multMax", 10000., "maximum multiplicity accepted"}; - - void init(o2::framework::InitContext&) - { - registry.add("hMassD0MCRecSig", "D0 signal 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("hMassD0MCRecRefl", "D0 reflection 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("hMassD0MCRecBkg", "D0 background 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("hMassD0barMCRecSig", "D0bar signal candidates - MC reco;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})"}}}); - registry.add("hMassD0barMCRecRefl", "D0bar reflection candidates - MC reco;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})"}}}); - registry.add("hMassD0barMCRecBkg", "D0bar background candidates - MC reco;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})"}}}); + } // end outer loop } - Filter filterSelectCandidates = (aod::hf_selcandidate_d0::isSelD0 >= selectionFlagD0 || aod::hf_selcandidate_d0::isSelD0bar >= selectionFlagD0bar); + PROCESS_SWITCH(HfCorrelatorD0D0bar, processData, "Process data", false); - void process(aod::Collision const& collision, soa::Join& tracks, soa::Filtered> const& candidates) + /// 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& tracks, soa::Filtered> const& candidates) { int nTracks = 0; if (collision.numContrib() > 1) { @@ -269,13 +256,13 @@ struct HfCorrelatorD0D0barMcRec { } registry.fill(HIST("hMultiplicity"), nTracks); - //MC reco level + // MC reco level bool flagD0Signal = false; bool flagD0Reflection = false; bool flagD0barSignal = false; bool flagD0barReflection = false; for (auto& candidate1 : candidates) { - //check decay channel flag for candidate1 + // check decay channel flag for candidate1 if (!(candidate1.hfflag() & 1 << DecayType::D0ToPiK)) { continue; } @@ -292,7 +279,7 @@ struct HfCorrelatorD0D0barMcRec { } if (std::abs(candidate1.flagMCMatchRec()) == 1 << DecayType::D0ToPiK) { - //fill per-candidate distributions from D0/D0bar true candidates + // 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()); @@ -301,9 +288,9 @@ struct HfCorrelatorD0D0barMcRec { registry.fill(HIST("hYMCRec"), YD0(candidate1)); registry.fill(HIST("hSelectionStatusMCRec"), candidate1.isSelD0bar() + (candidate1.isSelD0() * 2)); } - //fill invariant mass plots from D0/D0bar signal and background candidates - if (candidate1.isSelD0() >= selectionFlagD0) { //only reco as D0 - if (candidate1.flagMCMatchRec() == 1 << DecayType::D0ToPiK) { //also matched as D0 + // fill invariant mass plots from D0/D0bar signal and background candidates + if (candidate1.isSelD0() >= selectionFlagD0) { // only reco as D0 + if (candidate1.flagMCMatchRec() == 1 << DecayType::D0ToPiK) { // also matched as D0 registry.fill(HIST("hMassD0MCRecSig"), InvMassD0(candidate1), candidate1.pt(), efficiencyWeight); } else if (candidate1.flagMCMatchRec() == -(1 << DecayType::D0ToPiK)) { registry.fill(HIST("hMassD0MCRecRefl"), InvMassD0(candidate1), candidate1.pt(), efficiencyWeight); @@ -311,8 +298,8 @@ struct HfCorrelatorD0D0barMcRec { registry.fill(HIST("hMassD0MCRecBkg"), InvMassD0(candidate1), candidate1.pt(), efficiencyWeight); } } - if (candidate1.isSelD0bar() >= selectionFlagD0bar) { //only reco as D0bar - if (candidate1.flagMCMatchRec() == -(1 << DecayType::D0ToPiK)) { //also matched as D0bar + if (candidate1.isSelD0bar() >= selectionFlagD0bar) { // only reco as D0bar + if (candidate1.flagMCMatchRec() == -(1 << DecayType::D0ToPiK)) { // also matched as D0bar registry.fill(HIST("hMassD0barMCRecSig"), InvMassD0bar(candidate1), candidate1.pt(), efficiencyWeight); } else if (candidate1.flagMCMatchRec() == 1 << DecayType::D0ToPiK) { registry.fill(HIST("hMassD0barMCRecRefl"), InvMassD0bar(candidate1), candidate1.pt(), efficiencyWeight); @@ -321,34 +308,34 @@ struct HfCorrelatorD0D0barMcRec { } } - //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 + // D-Dbar correlation dedicated section + // if the candidate is selected ad D0, search for D0bar and evaluate correlations + if (candidate1.isSelD0() < selectionFlagD0) { // discard candidates not selected as D0 in outer loop continue; } - flagD0Signal = candidate1.flagMCMatchRec() == 1 << DecayType::D0ToPiK; //flagD0Signal 'true' if candidate1 matched to D0 (particle) - flagD0Reflection = candidate1.flagMCMatchRec() == -(1 << DecayType::D0ToPiK); //flagD0Reflection 'true' if candidate1, selected as D0 (particle), is matched to D0bar (antiparticle) + flagD0Signal = candidate1.flagMCMatchRec() == 1 << DecayType::D0ToPiK; // flagD0Signal 'true' if candidate1 matched to D0 (particle) + flagD0Reflection = candidate1.flagMCMatchRec() == -(1 << DecayType::D0ToPiK); // flagD0Reflection 'true' if candidate1, selected as D0 (particle), is matched to D0bar (antiparticle) for (auto& candidate2 : candidates) { - if (!(candidate2.hfflag() & 1 << DecayType::D0ToPiK)) { //check decay channel flag for candidate2 + if (!(candidate2.hfflag() & 1 << DecayType::D0ToPiK)) { // check decay channel flag for candidate2 continue; } - if (candidate2.isSelD0bar() < selectionFlagD0bar) { //discard candidates not selected as D0bar in inner loop + if (candidate2.isSelD0bar() < selectionFlagD0bar) { // discard candidates not selected as D0bar in inner loop continue; } - flagD0barSignal = candidate2.flagMCMatchRec() == -(1 << DecayType::D0ToPiK); //flagD0barSignal 'true' if candidate2 matched to D0bar (antiparticle) - flagD0barReflection = candidate2.flagMCMatchRec() == 1 << DecayType::D0ToPiK; //flagD0barReflection 'true' if candidate2, selected as D0bar (antiparticle), is matched to D0 (particle) + flagD0barSignal = candidate2.flagMCMatchRec() == -(1 << DecayType::D0ToPiK); // flagD0barSignal 'true' if candidate2 matched to D0bar (antiparticle) + flagD0barReflection = candidate2.flagMCMatchRec() == 1 << DecayType::D0ToPiK; // flagD0barReflection 'true' if candidate2, selected as D0bar (antiparticle), is matched to D0 (particle) if (cutYCandMax >= 0. && std::abs(YD0(candidate2)) > cutYCandMax) { continue; } if (cutPtCandMin >= 0. && candidate2.pt() < cutPtCandMin) { continue; } - //Excluding trigger self-correlations (possible in case of both mass hypotheses accepted) + // 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/ref, 2 = bkg/sig, 3 = ref/bkg, 4 = ref/ref, 5 = ref/sig, 6 = sig/bkg, 7 = sig/ref, 8 = sig/sig + // 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; } @@ -370,10 +357,10 @@ struct HfCorrelatorD0D0barMcRec { pairSignalStatus); double etaCut = 0.; double ptCut = 0.; - do { //fill pairs vs etaCut plot + do { // fill pairs vs etaCut plot ptCut = 0.; etaCut += incrementEtaCut; - do { //fill pairs vs etaCut plot + 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); } @@ -382,42 +369,19 @@ struct HfCorrelatorD0D0barMcRec { } while (etaCut < maxEtaCut - epsilon); } // end inner loop (Dbars) - } //end outer loop + } // end outer loop } -}; -/// D0-D0bar correlation pair builder - for MC gen-level analysis (no filter/selection, only true signal) -struct HfCorrelatorD0D0barMcGen { + PROCESS_SWITCH(HfCorrelatorD0D0bar, processMcRec, "Process MC Reco mode", true); - Produces entryD0D0barPair; - - HistogramRegistry registry{ - "registry", - {{"hMCEvtCount", "Event counter - MC gen;;entries", {HistType::kTH1F, {{1, -0.5, 0.5}}}}, - {"hPtCandMCGen", "D0,D0bar particles - MC gen;particle #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, - {"hEtaMCGen", "D0,D0bar particles - MC gen;particle #it{#eta};entries", {HistType::kTH1F, {{100, -5., 5.}}}}, - {"hPhiMCGen", "D0,D0bar particles - MC gen;particle #it{#varphi};entries", {HistType::kTH1F, {{32, 0., 2. * o2::constants::math::PI}}}}, - {"hYMCGen", "D0,D0bar candidates - MC gen;candidate #it{#y};entries", {HistType::kTH1F, {{100, -5., 5.}}}}, - {"hcountD0D0barPerEvent", "D0,D0bar particles - MC gen;Number per event;entries", {HistType::kTH1F, {{20, 0., 20.}}}}, - {"hDDbarVsEtaCut", "D0,D0bar pairs vs #eta cut of D mesons;#eta_{max};entries", {HistType::kTH2F, {{(int)(maxEtaCut / incrementEtaCut), 0., maxEtaCut}, {(int)(ptThresholdForMaxEtaCut / incrementPtThreshold), 0., ptThresholdForMaxEtaCut}}}}, - {"hDDbarVsDaughterEtaCut", "D0,D0bar pairs vs #eta cut on D daughters;#eta_{max};entries", {HistType::kTH2F, {{(int)(maxEtaCut / incrementEtaCut), 0., maxEtaCut}, {(int)(ptThresholdForMaxEtaCut / incrementPtThreshold), 0., ptThresholdForMaxEtaCut}}}}}}; - - Configurable cutYCandMax{"cutYCandMax", -1., "max. cand. rapidity"}; - Configurable cutPtCandMin{"cutPtCandMin", -1., "min. cand. pT"}; - Configurable> bins{"ptBinsForMass", std::vector{o2::analysis::hf_cuts_d0_topik::pTBins_v}, "pT bin limits for trigger counters"}; - - void init(o2::framework::InitContext&) - { - registry.add("hcountD0triggersMCGen", "D0 trigger particles - MC gen;;N of trigger D0", {HistType::kTH2F, {{1, -0.5, 0.5}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); - } - - void process(aod::McCollision const& mccollision, MCParticlesPlus const& particlesMC) + /// D0-D0bar correlation pair builder - for MC gen-level analysis (no filter/selection, only true signal) + void processMcGen(aod::McCollision const& mccollision, MCParticlesPlus const& particlesMC) { int counterD0D0bar = 0; registry.fill(HIST("hMCEvtCount"), 0); - //MC gen level + // 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! + // check if the particle is D0 or D0bar (for general plot filling and selection, so both cases are fine) - NOTE: decay channel is not probed! if (std::abs(particle1.pdgCode()) != pdg::Code::kD0) { continue; } @@ -434,14 +398,14 @@ struct HfCorrelatorD0D0barMcGen { registry.fill(HIST("hYMCGen"), yD); counterD0D0bar++; - //D-Dbar correlation dedicated section - //if it's a D0 particle, search for D0bar and evaluate correlations - if (particle1.pdgCode() != pdg::Code::kD0) { //just checking the particle PDG, not the decay channel (differently from Reco: you have a BR factor btw such levels!) + // D-Dbar correlation dedicated section + // if it's a D0 particle, search for D0bar and evaluate correlations + if (particle1.pdgCode() != pdg::Code::kD0) { // just checking the particle PDG, not the decay channel (differently from Reco: you have a BR factor btw such levels!) continue; } - registry.fill(HIST("hcountD0triggersMCGen"), 0, particle1.pt()); //to count trigger D0 (for normalisation) + registry.fill(HIST("hCountD0triggersMCGen"), 0, particle1.pt()); // to count trigger D0 (for normalisation) for (auto& particle2 : particlesMC) { - if (particle2.pdgCode() != pdg::Code::kD0bar) { //check that inner particle is D0bar + if (particle2.pdgCode() != pdg::Code::kD0bar) { // check that inner particle is D0bar continue; } if (cutYCandMax >= 0. && std::abs(RecoDecay::Y(array{particle2.px(), particle2.py(), particle2.pz()}, RecoDecay::getMassPDG(particle2.pdgCode()))) > cutYCandMax) { @@ -454,10 +418,13 @@ struct HfCorrelatorD0D0barMcGen { particle2.eta() - particle1.eta(), particle1.pt(), particle2.pt()); + entryD0D0barRecoInfo(1.864, + 1.864, + 8); // dummy information double etaCut = 0.; double ptCut = 0.; - //fill pairs vs etaCut plot + // fill pairs vs etaCut plot bool rightDecayChannels = false; if ((std::abs(particle1.flagMCMatchGen()) == 1 << DecayType::D0ToPiK) && (std::abs(particle2.flagMCMatchGen()) == 1 << DecayType::D0ToPiK)) { rightDecayChannels = true; @@ -465,11 +432,11 @@ struct HfCorrelatorD0D0barMcGen { do { ptCut = 0.; etaCut += incrementEtaCut; - do { //fill pairs vs etaCut plot - if (std::abs(particle1.eta()) < etaCut && std::abs(particle2.eta()) < etaCut && particle1.pt() > ptCut && particle2.pt() > ptCut) { //fill with D and Dbar acceptance checks + do { // fill pairs vs etaCut plot + if (std::abs(particle1.eta()) < etaCut && std::abs(particle2.eta()) < etaCut && particle1.pt() > ptCut && particle2.pt() > ptCut) { // fill with D and Dbar acceptance checks registry.fill(HIST("hDDbarVsEtaCut"), etaCut - epsilon, ptCut + epsilon); } - if (rightDecayChannels) { //fill with D and Dbar daughter particls acceptance checks + if (rightDecayChannels) { // fill with D and Dbar daughter particls acceptance checks double etaCandidate1Daughter1 = particle1.daughter0_as().eta(); double etaCandidate1Daughter2 = particle1.daughter1_as().eta(); double etaCandidate2Daughter1 = particle2.daughter0_as().eta(); @@ -483,325 +450,30 @@ struct HfCorrelatorD0D0barMcGen { ptCut += incrementPtThreshold; } while (ptCut < ptThresholdForMaxEtaCut - epsilon); } while (etaCut < maxEtaCut - epsilon); - } //end inner loop - } //end outer loop - registry.fill(HIST("hcountD0D0barPerEvent"), counterD0D0bar); - } -}; - -/// D0-D0bar correlation pair builder - LIKE SIGN - for real data and data-like analysis (i.e. reco-level w/o matching request via MC truth) -/// NOTE: At the moment, both dPhi-symmetrical correlation pairs (part1-part2 and part2-part1) are filled, -/// since we bin in pT and selecting as trigger the largest pT particle would bias the distributions w.r.t. the ULS case. -struct HfCorrelatorD0D0barLs { - - Produces entryD0D0barPair; - Produces entryD0D0barRecoInfo; - - HistogramRegistry registry{ - "registry", - //NOTE: use hMassD0 for trigger normalisation (S*0.955), and hMass2DCorrelationPairs (in final task) for 2D-sideband-subtraction purposes - {{"hPtCand", "D0,D0bar candidates;candidate #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, - {"hPtProng0", "D0,D0bar candidates;prong 0 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, - {"hPtProng1", "D0,D0bar candidates;prong 1 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, - {"hSelectionStatus", "D0,D0bar candidates;selection status;entries", {HistType::kTH1F, {{4, -0.5, 3.5}}}}, - {"hEta", "D0,D0bar candidates;candidate #it{#eta};entries", {HistType::kTH1F, {{100, -5., 5.}}}}, - {"hPhi", "D0,D0bar candidates;candidate #it{#varphi};entries", {HistType::kTH1F, {{32, 0., 2. * o2::constants::math::PI}}}}, - {"hY", "D0,D0bar candidates;candidate #it{#y};entries", {HistType::kTH1F, {{100, -5., 5.}}}}}}; - - Configurable selectionFlagD0{"selectionFlagD0", 1, "Selection Flag for D0"}; - Configurable selectionFlagD0bar{"selectionFlagD0bar", 1, "Selection Flag for D0bar"}; - Configurable cutYCandMax{"cutYCandMax", -1., "max. cand. rapidity"}; - Configurable cutPtCandMin{"cutPtCandMin", -1., "min. cand. pT"}; - Configurable> bins{"ptBinsForMass", std::vector{o2::analysis::hf_cuts_d0_topik::pTBins_v}, "pT bin limits for candidate mass plots"}; - - void init(o2::framework::InitContext&) - { - registry.add("hMass", "D0,D0bar candidates;inv. mass (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{120, 1.5848, 2.1848}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); - registry.add("hMassD0", "D0,D0bar candidates;inv. mass D0 only (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{120, 1.5848, 2.1848}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); - registry.add("hMassD0bar", "D0,D0bar candidates;inv. mass D0bar only (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{120, 1.5848, 2.1848}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); - } - - Filter filterSelectCandidates = (aod::hf_selcandidate_d0::isSelD0 >= selectionFlagD0 || aod::hf_selcandidate_d0::isSelD0bar >= selectionFlagD0bar); - - void process(aod::Collision const& collision, soa::Filtered> const& candidates) - { - for (auto& candidate1 : candidates) { - //check decay channel flag for candidate1 - if (!(candidate1.hfflag() & 1 << DecayType::D0ToPiK)) { - continue; - } - if (cutYCandMax >= 0. && std::abs(YD0(candidate1)) > cutYCandMax) { - continue; - } - if (cutPtCandMin >= 0. && candidate1.pt() < cutPtCandMin) { - continue; - } - //fill invariant mass plots and generic info from all D0/D0bar candidates - if (candidate1.isSelD0() >= selectionFlagD0) { - registry.fill(HIST("hMass"), InvMassD0(candidate1), candidate1.pt()); - registry.fill(HIST("hMassD0"), InvMassD0(candidate1), candidate1.pt()); - } - if (candidate1.isSelD0bar() >= selectionFlagD0bar) { - registry.fill(HIST("hMass"), InvMassD0bar(candidate1), candidate1.pt()); - registry.fill(HIST("hMassD0bar"), InvMassD0bar(candidate1), candidate1.pt()); - } - registry.fill(HIST("hPtCand"), candidate1.pt()); - registry.fill(HIST("hPtProng0"), candidate1.ptProng0()); - registry.fill(HIST("hPtProng1"), candidate1.ptProng1()); - registry.fill(HIST("hEta"), candidate1.eta()); - registry.fill(HIST("hPhi"), candidate1.phi()); - registry.fill(HIST("hY"), YD0(candidate1)); - registry.fill(HIST("hSelectionStatus"), candidate1.isSelD0bar() + (candidate1.isSelD0() * 2)); - - //D-Dbar correlation dedicated section - //For like-sign, first loop on both D0 and D0bars. First candidate is for sure a D0 and D0bars (checked before, so don't re-check anything on it) - for (auto& candidate2 : candidates) { - //check decay channel flag for candidate2 - if (!(candidate2.hfflag() & 1 << DecayType::D0ToPiK)) { - continue; - } - //for the associated, has to have smaller pT, and pass D0sel if trigger passes D0sel, or D0barsel if trigger passes D0barsel - if ((candidate1.isSelD0() >= selectionFlagD0 && candidate2.isSelD0() >= selectionFlagD0) || (candidate1.isSelD0bar() >= selectionFlagD0bar && candidate2.isSelD0bar() >= selectionFlagD0bar)) { - if (cutYCandMax >= 0. && std::abs(YD0(candidate2)) > cutYCandMax) { - continue; - } - if (cutPtCandMin >= 0. && candidate2.pt() < cutPtCandMin) { - continue; - } - //Excluding self-correlations - if (candidate1.mRowIndex == candidate2.mRowIndex) { - continue; - } - entryD0D0barPair(getDeltaPhi(candidate2.phi(), candidate1.phi()), - candidate2.eta() - candidate1.eta(), - candidate1.pt(), - candidate2.pt()); - entryD0D0barRecoInfo(InvMassD0(candidate1), - InvMassD0bar(candidate2), - 0); - } - //note: candidates selected as both D0 and D0bar are used, and considered in both situation (but not auto-correlated): reflections could play a relevant role. - //another, more restrictive, option, could be to consider only candidates selected with a single option (D0 xor D0bar) - } // end inner loop (Dbars) - } //end outer loop - } -}; - -/// D0-D0bar correlation pair builder - LIKE SIGN - for MC reco analysis (data-like but matching to true DO and D0bar) -/// NOTE: At the moment, both dPhi-symmetrical correlation pairs (part1-part2 and part2-part1) are filled, -/// since we bin in pT and selecting as trigger the largest pT particle would bias the distributions w.r.t. the ULS case. -struct HfCorrelatorD0D0barMcRecLs { - - Produces entryD0D0barPair; - Produces entryD0D0barRecoInfo; - - HistogramRegistry registry{ - "registry", - //NOTE: use hMassD0 for trigger normalisation (S*0.955), and hMass2DCorrelationPairs (in final task) for 2D-sideband-subtraction purposes - {{"hPtCandMCRec", "D0,D0bar candidates - MC reco;candidate #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, - {"hPtProng0MCRec", "D0,D0bar candidates - MC reco;prong 0 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, - {"hPtProng1MCRec", "D0,D0bar candidates - MC reco;prong 1 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, - {"hSelectionStatusMCRec", "D0,D0bar candidates - MC reco;selection status;entries", {HistType::kTH1F, {{4, -0.5, 3.5}}}}, - {"hEtaMCRec", "D0,D0bar candidates - MC reco;candidate #it{#eta};entries", {HistType::kTH1F, {{100, -5., 5.}}}}, - {"hPhiMCRec", "D0,D0bar candidates - MC reco;candidate #it{#varphi};entries", {HistType::kTH1F, {{32, 0., 2. * o2::constants::math::PI}}}}, - {"hYMCRec", "D0,D0bar candidates - MC reco;candidate #it{#y};entries", {HistType::kTH1F, {{100, -5., 5.}}}}}}; - - Configurable selectionFlagD0{"selectionFlagD0", 1, "Selection Flag for D0"}; - Configurable selectionFlagD0bar{"selectionFlagD0bar", 1, "Selection Flag for D0bar"}; - Configurable cutYCandMax{"cutYCandMax", -1., "max. cand. rapidity"}; - Configurable cutPtCandMin{"cutPtCandMin", -1., "min. cand. pT"}; - Configurable> bins{"ptBinsForMass", std::vector{o2::analysis::hf_cuts_d0_topik::pTBins_v}, "pT bin limits for candidate mass plots"}; - - void init(o2::framework::InitContext&) - { - registry.add("hMassD0MCRec", "D0,D0bar candidates - MC reco;inv. mass (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{120, 1.5848, 2.1848}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); - registry.add("hMassD0barMCRec", "D0,D0bar candidates - MC reco;inv. mass D0 only (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{120, 1.5848, 2.1848}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); - } - - Filter filterSelectCandidates = (aod::hf_selcandidate_d0::isSelD0 >= selectionFlagD0 || aod::hf_selcandidate_d0::isSelD0bar >= selectionFlagD0bar); - - void process(aod::Collision const& collision, soa::Filtered> const& candidates) - { - //MC reco level - for (auto& candidate1 : candidates) { - //check decay channel flag for candidate1 - if (!(candidate1.hfflag() & 1 << DecayType::D0ToPiK)) { - continue; - } - if (cutYCandMax >= 0. && std::abs(YD0(candidate1)) > cutYCandMax) { - continue; - } - if (cutPtCandMin >= 0. && candidate1.pt() < cutPtCandMin) { - continue; - } - if (std::abs(candidate1.flagMCMatchRec()) == 1 << DecayType::D0ToPiK) { - //fill invariant mass plots and generic info from all D0/D0bar candidates - if (candidate1.isSelD0() >= selectionFlagD0 && candidate1.flagMCMatchRec() == DecayType::D0ToPiK) { //only reco and matched as D0 - registry.fill(HIST("hMassD0MCRec"), InvMassD0(candidate1)); - } - if (candidate1.isSelD0bar() >= selectionFlagD0bar && candidate1.flagMCMatchRec() == DecayType::D0ToPiK) { //only reco and matched as D0bar - registry.fill(HIST("hMassD0barMCRec"), InvMassD0bar(candidate1)); - } - registry.fill(HIST("hPtCandMCRec"), candidate1.pt()); - registry.fill(HIST("hPtProng0MCRec"), candidate1.ptProng0()); - registry.fill(HIST("hPtProng1MCRec"), candidate1.ptProng1()); - registry.fill(HIST("hEtaMCRec"), candidate1.eta()); - registry.fill(HIST("hPhiMCRec"), candidate1.phi()); - registry.fill(HIST("hYMCRec"), YD0(candidate1)); - registry.fill(HIST("hSelectionStatusMCRec"), candidate1.isSelD0bar() + (candidate1.isSelD0() * 2)); - - //D-Dbar correlation dedicated section - //For like-sign, first loop on both D0 and D0bars. First candidate is for sure a D0 and D0bars (looping on filtered) and was already matched, so don't re-check anything on it) - for (auto& candidate2 : candidates) { - //check decay channel flag for candidate2 - if (!(candidate2.hfflag() & 1 << DecayType::D0ToPiK)) { - continue; - } - bool conditionLSForD0 = (candidate1.isSelD0() >= selectionFlagD0bar && candidate1.flagMCMatchRec() == 1 << DecayType::D0ToPiK) && (candidate2.isSelD0() >= selectionFlagD0bar && candidate2.flagMCMatchRec() == 1 << DecayType::D0ToPiK); - bool conditionLSForD0bar = (candidate1.isSelD0bar() >= selectionFlagD0bar && candidate1.flagMCMatchRec() == -(1 << DecayType::D0ToPiK)) && (candidate2.isSelD0bar() >= selectionFlagD0bar && candidate2.flagMCMatchRec() == -(1 << DecayType::D0ToPiK)); - if (conditionLSForD0 || conditionLSForD0bar) { //LS pair (of D0 or of D0bar) + pt2= 0. && std::abs(YD0(candidate2)) > cutYCandMax) { - continue; - } - if (cutPtCandMin >= 0. && candidate2.pt() < cutPtCandMin) { - continue; - } - //Excluding self-correlations - if (candidate1.mRowIndex == candidate2.mRowIndex) { - continue; - } - entryD0D0barPair(getDeltaPhi(candidate2.phi(), candidate1.phi()), - candidate2.eta() - candidate1.eta(), - candidate1.pt(), - candidate2.pt()); - entryD0D0barRecoInfo(InvMassD0(candidate1), - InvMassD0bar(candidate2), - 0); //for LS studies we set a dummy 0 for pairSignalStatus (there are no more the usual 4 possible combinations) - - } //end inner if (MC match) - } // end inner loop (Dbars) - } //end outer if (MC match) - } //end outer loop + } // end inner loop + } // end outer loop + registry.fill(HIST("hCountD0D0barPerEvent"), counterD0D0bar); } -}; - -/// D0-D0bar correlation pair builder - for MC gen-level analysis, like sign particles -/// NOTE: At the moment, both dPhi-symmetrical correlation pairs (part1-part2 and part2-part1) are filled, -/// since we bin in pT and selecting as trigger the largest pT particle would bias the distributions w.r.t. the ULS case. -struct HfCorrelatorD0D0barMcGenLs { - Produces entryD0D0barPair; - - HistogramRegistry registry{ - "registry", - {{"hMCEvtCount", "Event counter - MC gen;;entries", {HistType::kTH1F, {{1, -0.5, 0.5}}}}, - {"hPtCandMCGen", "D0,D0bar particles - MC gen;particle #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, - {"hEtaMCGen", "D0,D0bar particles - MC gen;particle #it{#eta};entries", {HistType::kTH1F, {{100, -5., 5.}}}}, - {"hPhiMCGen", "D0,D0bar particles - MC gen;particle #it{#varphi};entries", {HistType::kTH1F, {{32, 0., 2. * o2::constants::math::PI}}}}, - {"hYMCGen", "D0,D0bar candidates - MC gen;candidate #it{#y};entries", {HistType::kTH1F, {{100, -5., 5.}}}}, - {"hcountD0D0barPerEvent", "D0,D0bar particles - MC gen;Number per event;entries", {HistType::kTH1F, {{20, 0., 20.}}}}}}; - - Configurable cutYCandMax{"cutYCandMax", -1., "max. cand. rapidity"}; - Configurable cutPtCandMin{"cutPtCandMin", -1., "min. cand. pT"}; - Configurable> bins{"ptBinsForMass", std::vector{o2::analysis::hf_cuts_d0_topik::pTBins_v}, "pT bin limits for trigger counters"}; - - void init(o2::framework::InitContext&) - { - registry.add("hcountD0triggersMCGen", "D0 trigger particles - MC gen;;N of trigger D0", {HistType::kTH2F, {{1, -0.5, 0.5}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); - } + PROCESS_SWITCH(HfCorrelatorD0D0bar, processMcGen, "Process MC Gen mode", false); - void process(aod::McCollision const& mccollision, MCParticlesPlus const& particlesMC) + /// c-cbar correlator table builder - for MC gen-level analysis + void processccbar(aod::McCollision const& mccollision, MCParticlesPlus const& particlesMC) { - int counterD0D0bar = 0; registry.fill(HIST("hMCEvtCount"), 0); - //MC gen level - for (auto& particle1 : particlesMC) { - //check if the particle is D0 or D0bar (both can be trigger) - NOTE: decay channel is not probed! - if (std::abs(particle1.pdgCode()) != pdg::Code::kD0) { - continue; - } - double yD = RecoDecay::Y(array{particle1.px(), particle1.py(), particle1.pz()}, RecoDecay::getMassPDG(particle1.pdgCode())); - if (cutYCandMax >= 0. && std::abs(yD) > cutYCandMax) { - continue; - } - if (cutPtCandMin >= 0. && particle1.pt() < cutPtCandMin) { - continue; - } - - registry.fill(HIST("hPtCandMCGen"), particle1.pt()); - registry.fill(HIST("hEtaMCGen"), particle1.eta()); - registry.fill(HIST("hPhiMCGen"), particle1.phi()); - registry.fill(HIST("hYMCGen"), yD); - counterD0D0bar++; - //D-Dbar correlation dedicated section - //if it's D0, search for D0bar and evaluate correlations. - registry.fill(HIST("hcountD0triggersMCGen"), 0, particle1.pt()); //to count trigger D0 (normalisation) - for (auto& particle2 : particlesMC) { - if (std::abs(particle2.pdgCode()) != pdg::Code::kD0) { //check that associated is a D0/D0bar (both are fine) - continue; - } - if (cutYCandMax >= 0. && std::abs(RecoDecay::Y(array{particle2.px(), particle2.py(), particle2.pz()}, RecoDecay::getMassPDG(particle2.pdgCode()))) > cutYCandMax) { - continue; - } - if (cutPtCandMin >= 0. && particle2.pt() < cutPtCandMin) { - continue; - } - if (particle2.pdgCode() == particle1.pdgCode()) { //like-sign condition (both 421 or both -421) and pT_Trig>pT_assoc - //Excluding self-correlations - if (particle1.mRowIndex == particle2.mRowIndex) { - continue; - } - entryD0D0barPair(getDeltaPhi(particle2.phi(), particle1.phi()), - particle2.eta() - particle1.eta(), - particle1.pt(), - particle2.pt()); - } - } // end inner loop (Dbars) - } //end outer loop - registry.fill(HIST("hcountD0D0barPerEvent"), counterD0D0bar); - } -}; - -/// c-cbar correlator table builder - for MC gen-level analysis -struct HfCorrelatorCCbarMcGen { - - Produces entryccbarPair; + int counterCCbar = 0, counterCCbarBeforeEtasel = 0; - HistogramRegistry registry{ - "registry", - {{"hMCEvtCount", "Event counter - MC gen;;entries", {HistType::kTH1F, {{1, -0.5, 0.5}}}}, - {"hPtCandMCGen", "c,cbar particles - MC gen;particle #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, - {"hEtaMCGen", "c,cbar particles - MC gen;particle #it{#eta};entries", {HistType::kTH1F, {{100, -5., 5.}}}}, - {"hYMCGen", "c,cbar candidates - MC gen;candidate #it{#y};entries", {HistType::kTH1F, {{100, -5., 5.}}}}, - {"hPhiMCGen", "c,cbar particles - MC gen;particle #it{#varphi};entries", {HistType::kTH1F, {{32, 0., 2. * o2::constants::math::PI}}}}, - {"hcountCCbarPerEvent", "c,cbar particles - MC gen;Number per event;entries", {HistType::kTH1F, {{20, 0., 20.}}}}, - {"hcountCCbarPerEventPreEtaCut", "c,cbar particles - MC gen;Number per event pre #eta cut;entries", {HistType::kTH1F, {{20, 0., 20.}}}}}}; - - Configurable cutYCandMax{"cutYCandMax", -1., "max. cand. rapidity"}; - Configurable cutPtCandMin{"cutPtCandMin", -1., "min. cand. pT"}; - Configurable> bins{"ptBinsForMass", std::vector{o2::analysis::hf_cuts_d0_topik::pTBins_v}, "pT bin limits for trigger counters"}; - - void init(o2::framework::InitContext&) - { - registry.add("hcountCtriggersMCGen", "c trigger particles - MC gen;;N of trigger c quark", {HistType::kTH2F, {{1, -0.5, 0.5}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); - } - - void process(aod::McCollision const& mccollision, MCParticlesPlus const& particlesMC) - { - registry.fill(HIST("hMCEvtCount"), 0); - int counterccbar = 0, counterccbarPreEtasel = 0; - - //loop over particles at MC gen level + // loop over particles at MC gen level for (auto& particle1 : particlesMC) { - if (std::abs(particle1.pdgCode()) != PDG_t::kCharm) { //search c or cbar particles + if (std::abs(particle1.pdgCode()) != PDG_t::kCharm) { // search c or cbar particles continue; } int partMothPDG = particle1.mother0_as().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 + // 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) + counterCCbarBeforeEtasel++; // count c or cbar (before kinematic selection) double yC = RecoDecay::Y(array{particle1.px(), particle1.py(), particle1.pz()}, RecoDecay::getMassPDG(particle1.pdgCode())); if (cutYCandMax >= 0. && std::abs(yC) > cutYCandMax) { continue; @@ -813,17 +485,17 @@ struct HfCorrelatorCCbarMcGen { registry.fill(HIST("hEtaMCGen"), particle1.eta()); registry.fill(HIST("hPhiMCGen"), particle1.phi()); registry.fill(HIST("hYMCGen"), yC); - counterccbar++; //count if c or cbar don't come from themselves during fragmentation (after kinematic selection) + 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. + // c-cbar correlation dedicated section + // if it's c, search for cbar and evaluate correlations. if (particle1.pdgCode() != PDG_t::kCharm) { continue; } - registry.fill(HIST("hcountCtriggersMCGen"), 0, particle1.pt()); //to count trigger c quark (for normalisation) + registry.fill(HIST("hCountCtriggersMCGen"), 0, particle1.pt()); // to count trigger c quark (for normalisation) for (auto& particle2 : particlesMC) { - if (particle2.pdgCode() != PDG_t::kCharmBar) { //check that inner particle is a cbar + if (particle2.pdgCode() != PDG_t::kCharmBar) { // check that inner particle is a cbar continue; } if (cutYCandMax >= 0. && std::abs(RecoDecay::Y(array{particle2.px(), particle2.py(), particle2.pz()}, RecoDecay::getMassPDG(particle2.pdgCode()))) > cutYCandMax) { @@ -832,136 +504,27 @@ struct HfCorrelatorCCbarMcGen { if (cutPtCandMin >= 0. && particle2.pt() < cutPtCandMin) { continue; } - //check whether mothers of quark cbar (from associated loop) are still '-4' particles - in that case the cbar quark comes from its own fragmentation, skip it + // 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 (particle2.mother0_as().pdgCode() == PDG_t::kCharmBar) { continue; } - entryccbarPair(getDeltaPhi(particle2.phi(), particle1.phi()), - particle2.eta() - particle1.eta(), - particle1.pt(), - particle2.pt()); - } // end inner loop - } //end outer loop - registry.fill(HIST("hcountCCbarPerEvent"), counterccbar); - registry.fill(HIST("hcountCCbarPerEventPreEtaCut"), counterccbarPreEtasel); - } -}; - -/// c-cbar correlator table builder - for MC gen-level analysis - Like Sign -struct HfCorrelatorCCbarMcGenLs { - - Produces entryccbarPair; - - HistogramRegistry registry{ - "registry", - {{"hMCEvtCount", "Event counter - MC gen;;entries", {HistType::kTH1F, {{1, -0.5, 0.5}}}}, - {"hPtCandMCGen", "c,cbar particles - MC gen;particle #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, - {"hEtaMCGen", "c,cbar particles - MC gen;particle #it{#eta};entries", {HistType::kTH1F, {{100, -5., 5.}}}}, - {"hYMCGen", "c,cbar candidates - MC gen;candidate #it{#y};entries", {HistType::kTH1F, {{100, -5., 5.}}}}, - {"hPhiMCGen", "c,cbar particles - MC gen;particle #it{#varphi};entries", {HistType::kTH1F, {{32, 0., 2. * o2::constants::math::PI}}}}, - {"hcountCCbarPerEvent", "c,cbar particles - MC gen;Number per event;entries", {HistType::kTH1F, {{20, 0., 20.}}}}, - {"hcountCCbarPerEventPreEtaCut", "c,cbar particles - MC gen;Number per event pre #eta cut;entries", {HistType::kTH1F, {{20, 0., 20.}}}}}}; - - Configurable cutYCandMax{"cutYCandMax", -1., "max. cand. rapidity"}; - Configurable cutPtCandMin{"cutPtCandMin", -1., "min. cand. pT"}; - Configurable> bins{"ptBinsForMass", std::vector{o2::analysis::hf_cuts_d0_topik::pTBins_v}, "pT bin limits for trigger counters"}; - - void init(o2::framework::InitContext&) - { - registry.add("hcountCtriggersMCGen", "c trigger particles - MC gen;;N of trigger c quark", {HistType::kTH2F, {{1, -0.5, 0.5}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); - } - - void process(aod::McCollision const& mccollision, MCParticlesPlus const& particlesMC) - { - registry.fill(HIST("hMCEvtCount"), 0); - int counterccbar = 0, counterccbarPreEtasel = 0; - - //loop over particles at MC gen level - for (auto& particle1 : particlesMC) { - if (std::abs(particle1.pdgCode()) != PDG_t::kCharm) { //search c or cbar particles - continue; - } - int partMothPDG = particle1.mother0_as().pdgCode(); - //check whether mothers of quark c/cbar are still '4'/'-4' particles - in that case the c/cbar quark comes from its own fragmentation, skip it - if (partMothPDG == particle1.pdgCode()) { - continue; - } - counterccbarPreEtasel++; //count c or cbar (before kinematic selection) - double yC = RecoDecay::Y(array{particle1.px(), particle1.py(), particle1.pz()}, RecoDecay::getMassPDG(particle1.pdgCode())); - if (cutYCandMax >= 0. && std::abs(yC) > cutYCandMax) { - continue; - } - if (cutPtCandMin >= 0. && particle1.pt() < cutPtCandMin) { - continue; - } - registry.fill(HIST("hPtCandMCGen"), particle1.pt()); - registry.fill(HIST("hEtaMCGen"), particle1.eta()); - registry.fill(HIST("hPhiMCGen"), particle1.phi()); - registry.fill(HIST("hYMCGen"), yC); - counterccbar++; //count if c or cbar don't come from themselves during fragmentation (after kinematic selection) - - //c-cbar correlation dedicated section - registry.fill(HIST("hcountCtriggersMCGen"), 0, particle1.pt()); //to count trigger c quark (for normalisation) - - for (auto& particle2 : particlesMC) { - if (std::abs(particle2.pdgCode()) != PDG_t::kCharm) { //search c or cbar for associated particles - continue; - } - if (cutYCandMax >= 0. && std::abs(RecoDecay::Y(array{particle2.px(), particle2.py(), particle2.pz()}, RecoDecay::getMassPDG(particle2.pdgCode()))) > cutYCandMax) { - continue; - } - if (cutPtCandMin >= 0. && particle2.pt() < cutPtCandMin) { - continue; - } - if (particle2.pdgCode() == particle1.pdgCode()) { - //check whether mothers of quark cbar (from associated loop) are still '-4' particles - in that case the cbar quark comes from its own fragmentation, skip it - if (particle2.mother0_as().pdgCode() == particle2.pdgCode()) { - continue; - } - //Excluding self-correlations - if (particle1.mRowIndex == particle2.mRowIndex) { - continue; - } - entryccbarPair(getDeltaPhi(particle2.phi(), particle1.phi()), + 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); + entryD0D0barRecoInfo(1.864, + 1.864, + 8); // dummy information + } // end inner loop + } // end outer loop + registry.fill(HIST("hCountCCbarPerEvent"), counterCCbar); + registry.fill(HIST("hCountCCbarPerEventBeforeEtaCut"), counterCCbarBeforeEtasel); } + + PROCESS_SWITCH(HfCorrelatorD0D0bar, processccbar, "Process ccbar pairs", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { - WorkflowSpec workflow{}; - const bool doMCccbar = cfgc.options().get("doMCccbar"); - const bool doMCGen = cfgc.options().get("doMCGen"); - const bool doMCRec = cfgc.options().get("doMCRec"); - const bool doLikeSign = cfgc.options().get("doLikeSign"); - if (!doLikeSign) { //unlike-sign analyses - if (doMCGen) { //MC-Gen analysis - workflow.push_back(adaptAnalysisTask(cfgc)); - } else if (doMCRec) { //MC-Reco analysis - workflow.push_back(adaptAnalysisTask(cfgc)); - } else if (doMCccbar) { //MC-Reco analysis - workflow.push_back(adaptAnalysisTask(cfgc)); - } else { //data analysis - workflow.push_back(adaptAnalysisTask(cfgc)); - } - } else { //like-sign analyses - if (doMCGen) { //MC-Gen analysis - workflow.push_back(adaptAnalysisTask(cfgc)); - } else if (doMCRec) { //MC-Reco analysis - workflow.push_back(adaptAnalysisTask(cfgc)); - } else if (doMCccbar) { //MC-Reco analysis - workflow.push_back(adaptAnalysisTask(cfgc)); - } else { //data analysis - workflow.push_back(adaptAnalysisTask(cfgc)); - } - } - - return workflow; + return WorkflowSpec{adaptAnalysisTask(cfgc)}; } diff --git a/PWGHF/Tasks/HFCorrelatorD0D0barBarrelFullPID.cxx b/PWGHF/Tasks/HFCorrelatorD0D0barBarrelFullPID.cxx index 31e323dce8d..e216dcc692d 100644 --- a/PWGHF/Tasks/HFCorrelatorD0D0barBarrelFullPID.cxx +++ b/PWGHF/Tasks/HFCorrelatorD0D0barBarrelFullPID.cxx @@ -9,7 +9,7 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. -/// \file HFCorrelatorD0D0barBarrelFullPID.cxx +/// \file HFCorrelatorD0D0barBarrelFullPid.cxx /// \brief Temporary D0-D0bar correlator task with full barrel PID implementation - data-like, MC-reco and MC-kine analyses. For ULS and LS pairs /// /// \author Fabio Colamaria , INFN Bari @@ -29,18 +29,6 @@ using namespace o2::aod::hf_correlation_ddbar; using namespace o2::analysis::hf_cuts_d0_topik; using namespace o2::constants::math; -void customize(std::vector& workflowOptions) -{ - ConfigParamSpec optionDoLikeSign{"doLikeSign", VariantType::Bool, false, {"Run Like-Sign analysis."}}; - ConfigParamSpec optionDoMCccbar{"doMCccbar", VariantType::Bool, false, {"Run MC-Gen dedicated tasks."}}; - ConfigParamSpec optionDoMCGen{"doMCGen", VariantType::Bool, false, {"Run MC-Gen dedicated tasks."}}; - ConfigParamSpec optionDoMCRec{"doMCRec", VariantType::Bool, true, {"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" /// @@ -62,26 +50,55 @@ const int npTBinsMassAndEfficiency = o2::analysis::hf_cuts_d0_topik::npTBins; 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 = 2. * o2::constants::math::PI; +const int yAxisBins = 100; +const double yAxisMin = -5.; +const double yAxisMax = 5.; +const int ptDAxisBins = 180; +const double ptDAxisMin = 0.; +const double ptDAxisMax = 36.; + using MCParticlesPlus = soa::Join; -/// D0-D0bar correlation pair builder - for real data and data-like analysis (i.e. reco-level w/o matching request via MC truth) -struct HFCorrelatorD0D0barBarrelFullPID { +struct HfCorrelatorD0D0barBarrelFullPid { 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.}}}}, + // 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, {{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.}}}}, + {"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}}}}, {"hMultiplicityPreSelection", "multiplicity prior to selection;multiplicity;entries", {HistType::kTH1F, {{10000, 0., 10000.}}}}, {"hMultiplicity", "multiplicity;multiplicity;entries", {HistType::kTH1F, {{10000, 0., 10000.}}}}, - {"hDDbarVsEtaCut", "D0,D0bar pairs vs #eta cut;#eta_{max};entries", {HistType::kTH2F, {{(int)(maxEtaCut / incrementEtaCut), 0., maxEtaCut}, {(int)(ptThresholdForMaxEtaCut / incrementPtThreshold), 0., ptThresholdForMaxEtaCut}}}}}}; + {"hDDbarVsEtaCut", "D0,D0bar pairs vs #eta cut;#eta_{max};candidates #it{p}_{T} threshold (GeV/#it{c});entries", {HistType::kTH2F, {{(int)(maxEtaCut / incrementEtaCut), 0., maxEtaCut}, {(int)(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}}}}, + {"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, {{ptDAxisBins, ptDAxisMin, ptDAxisMax}}}}, + {"hEtaMCGen", "D0,D0bar particles - MC gen;particle #it{#eta};entries", {HistType::kTH1F, {{yAxisBins, yAxisMin, yAxisMax}}}}, + {"hPhiMCGen", "D0,D0bar particles - MC gen;particle #it{#varphi};entries", {HistType::kTH1F, {{phiAxisBins, phiAxisMin, phiAxisMax}}}}, + {"hYMCGen", "D0,D0bar candidates - MC gen;candidate #it{y};entries", {HistType::kTH1F, {{yAxisBins, yAxisMin, yAxisMax}}}}, + {"hCountD0D0barPerEvent", "D0,D0bar particles - MC gen;Number per event;entries", {HistType::kTH1F, {{20, 0., 20.}}}}, + {"hDDbarVsDaughterEtaCut", "D0,D0bar pairs vs #eta cut on D daughters;#eta_{max};candidates #it{p}_{T} threshold (GeV/#it{c});entries", {HistType::kTH2F, {{(int)(maxEtaCut / incrementEtaCut), 0., maxEtaCut}, {(int)(ptThresholdForMaxEtaCut / incrementPtThreshold), 0., ptThresholdForMaxEtaCut}}}}, + {"hCountCCbarPerEvent", "c,cbar particles - MC gen;Number per event;entries", {HistType::kTH1F, {{20, 0., 20.}}}}, + {"hCountCCbarPerEventBeforeEtaCut", "c,cbar particles - MC gen;Number per event pre #eta cut;entries", {HistType::kTH1F, {{20, 0., 20.}}}}}}; Configurable selectionFlagD0{"selectionFlagD0", 1, "Selection Flag for D0"}; Configurable selectionFlagD0bar{"selectionFlagD0bar", 1, "Selection Flag for D0bar"}; @@ -95,14 +112,24 @@ struct HFCorrelatorD0D0barBarrelFullPID { 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})"}}}); + auto vbins = (std::vector)bins; + 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})"}}}); + registry.add("hCountD0triggersMCGen", "D0 trigger particles - MC gen;;N of trigger D0", {HistType::kTH2F, {{1, -0.5, 0.5}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hCountCtriggersMCGen", "c trigger particles - MC gen;;N of trigger c quark", {HistType::kTH2F, {{1, -0.5, 0.5}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); } Filter filterSelectCandidates = (aod::hf_selcandidate_d0_ALICE3_Barrel::isSelD0TOFplusRICHPID >= selectionFlagD0 || aod::hf_selcandidate_d0_ALICE3_Barrel::isSelD0barTOFplusRICHPID >= selectionFlagD0bar); - void process(aod::Collision const& collision, soa::Join& tracks, soa::Filtered> const& candidates) + /// 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& tracks, soa::Filtered> const& candidates) { int nTracks = 0; if (collision.numContrib() > 1) { @@ -129,7 +156,7 @@ struct HFCorrelatorD0D0barBarrelFullPID { if (cutPtCandMin >= 0. && candidate1.pt() < cutPtCandMin) { continue; } - //check decay channel flag for candidate1 + // check decay channel flag for candidate1 if (!(candidate1.hfflag() & 1 << DecayType::D0ToPiK)) { continue; } @@ -139,7 +166,7 @@ struct HFCorrelatorD0D0barBarrelFullPID { efficiencyWeight = 1. / efficiencyDmeson->at(o2::analysis::findBin(bins, candidate1.pt())); } - //fill invariant mass plots and generic info from all D0/D0bar candidates + // fill invariant mass plots and generic info from all D0/D0bar candidates if (candidate1.isSelD0TOFplusRICHPID() >= selectionFlagD0) { registry.fill(HIST("hMass"), InvMassD0(candidate1), candidate1.pt(), efficiencyWeight); registry.fill(HIST("hMassD0"), InvMassD0(candidate1), candidate1.pt(), efficiencyWeight); @@ -156,26 +183,26 @@ struct HFCorrelatorD0D0barBarrelFullPID { registry.fill(HIST("hY"), YD0(candidate1)); registry.fill(HIST("hSelectionStatus"), candidate1.isSelD0barTOFplusRICHPID() + (candidate1.isSelD0TOFplusRICHPID() * 2)); - //D-Dbar correlation dedicated section - //if the candidate is a D0, search for D0bar and evaluate correlations + // D-Dbar correlation dedicated section + // if the candidate is a D0, search for D0bar and evaluate correlations if (candidate1.isSelD0TOFplusRICHPID() < selectionFlagD0) { continue; } for (auto& candidate2 : candidates) { - if (!(candidate2.hfflag() & 1 << DecayType::D0ToPiK)) { //check decay channel flag for candidate2 + if (!(candidate2.hfflag() & 1 << DecayType::D0ToPiK)) { // check decay channel flag for candidate2 continue; } - if (candidate2.isSelD0barTOFplusRICHPID() < selectionFlagD0bar) { //keep only D0bar candidates passing the selection + if (candidate2.isSelD0barTOFplusRICHPID() < selectionFlagD0bar) { // keep only D0bar candidates passing the selection continue; } - //kinematic selection on D0bar candidates + // kinematic selection on D0bar candidates if (cutYCandMax >= 0. && std::abs(YD0(candidate2)) > cutYCandMax) { continue; } if (cutPtCandMin >= 0. && candidate2.pt() < cutPtCandMin) { continue; } - //excluding trigger self-correlations (possible in case of both mass hypotheses accepted) + // excluding trigger self-correlations (possible in case of both mass hypotheses accepted) if (candidate1.mRowIndex == candidate2.mRowIndex) { continue; } @@ -188,68 +215,28 @@ struct HFCorrelatorD0D0barBarrelFullPID { 0); double etaCut = 0.; double ptCut = 0.; - do { //fill pairs vs etaCut plot + do { // fill pairs vs etaCut plot ptCut = 0.; etaCut += incrementEtaCut; - do { //fill pairs vs etaCut plot + 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) + // 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 HFCorrelatorD0D0barBarrelFullPIDMcRec { - - Produces entryD0D0barPair; - Produces entryD0D0barRecoInfo; - - HistogramRegistry registry{ - "registry", - //NOTE: use hMassD0 for trigger normalisation (S*0.955), and hMass2DCorrelationPairs (in final task) for 2D-sideband-subtraction purposes - {{"hPtCandMCRec", "D0,D0bar candidates - MC reco;candidate #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, - {"hPtProng0MCRec", "D0,D0bar candidates - MC reco;prong 0 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, - {"hPtProng1MCRec", "D0,D0bar candidates - MC reco;prong 1 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, - {"hSelectionStatusMCRec", "D0,D0bar candidates - MC reco;selection status;entries", {HistType::kTH1F, {{4, -0.5, 3.5}}}}, - {"hEtaMCRec", "D0,D0bar candidates - MC reco;candidate #it{#eta};entries", {HistType::kTH1F, {{100, -5., 5.}}}}, - {"hPhiMCRec", "D0,D0bar candidates - MC reco;candidate #it{#varphi};entries", {HistType::kTH1F, {{32, 0., 2. * o2::constants::math::PI}}}}, - {"hYMCRec", "D0,D0bar candidates - MC reco;candidate #it{#y};entries", {HistType::kTH1F, {{100, -5., 5.}}}}, - {"hMultiplicityPreSelection", "multiplicity prior to selection;multiplicity;entries", {HistType::kTH1F, {{10000, 0., 10000.}}}}, - {"hMultiplicity", "multiplicity;multiplicity;entries", {HistType::kTH1F, {{10000, 0., 10000.}}}}, - {"hDDbarVsEtaCut", "D0,D0bar pairs vs #eta cut;#eta_{max};entries", {HistType::kTH2F, {{(int)(maxEtaCut / incrementEtaCut), 0., maxEtaCut}, {(int)(ptThresholdForMaxEtaCut / incrementPtThreshold), 0., ptThresholdForMaxEtaCut}}}}}}; - - Configurable selectionFlagD0{"selectionFlagD0", 1, "Selection Flag for D0"}; - Configurable selectionFlagD0bar{"selectionFlagD0bar", 1, "Selection Flag for D0bar"}; - Configurable cutYCandMax{"cutYCandMax", -1., "max. cand. rapidity"}; - Configurable cutPtCandMin{"cutPtCandMin", -1., "min. cand. pT"}; - Configurable> bins{"ptBinsForMassAndEfficiency", std::vector{o2::analysis::hf_cuts_d0_topik::pTBins_v}, "pT bin limits for candidate mass plots and efficiency"}; - Configurable> efficiencyDmeson{"efficiencyDmeson", std::vector{efficiencyDmeson_v}, "Efficiency values for D0 meson"}; - Configurable flagApplyEfficiency{"efficiencyFlagD", 1, "Flag for applying D-meson efficiency weights"}; - Configurable multMin{"multMin", 0., "minimum multiplicity accepted"}; - Configurable multMax{"multMax", 10000., "maximum multiplicity accepted"}; - - void init(o2::framework::InitContext&) - { - registry.add("hMassD0MCRecSig", "D0 signal 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("hMassD0MCRecRefl", "D0 reflection 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("hMassD0MCRecBkg", "D0 background 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("hMassD0barMCRecSig", "D0bar signal candidates - MC reco;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})"}}}); - registry.add("hMassD0barMCRecRefl", "D0bar reflection candidates - MC reco;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})"}}}); - registry.add("hMassD0barMCRecBkg", "D0bar background candidates - MC reco;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})"}}}); + } // end outer loop } - Filter filterSelectCandidates = (aod::hf_selcandidate_d0_ALICE3_Barrel::isSelD0TOFplusRICHPID >= selectionFlagD0 || aod::hf_selcandidate_d0_ALICE3_Barrel::isSelD0barTOFplusRICHPID >= selectionFlagD0bar); + PROCESS_SWITCH(HfCorrelatorD0D0barBarrelFullPid, processData, "Process data", false); - void process(aod::Collision const& collision, soa::Join& tracks, soa::Filtered> const& candidates) + /// 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& tracks, soa::Filtered> const& candidates) { int nTracks = 0; if (collision.numContrib() > 1) { @@ -269,13 +256,13 @@ struct HFCorrelatorD0D0barBarrelFullPIDMcRec { } registry.fill(HIST("hMultiplicity"), nTracks); - //MC reco level + // MC reco level bool flagD0Signal = false; bool flagD0Reflection = false; bool flagD0barSignal = false; bool flagD0barReflection = false; for (auto& candidate1 : candidates) { - //check decay channel flag for candidate1 + // check decay channel flag for candidate1 if (!(candidate1.hfflag() & 1 << DecayType::D0ToPiK)) { continue; } @@ -292,7 +279,7 @@ struct HFCorrelatorD0D0barBarrelFullPIDMcRec { } if (std::abs(candidate1.flagMCMatchRec()) == 1 << DecayType::D0ToPiK) { - //fill per-candidate distributions from D0/D0bar true candidates + // 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()); @@ -301,9 +288,9 @@ struct HFCorrelatorD0D0barBarrelFullPIDMcRec { registry.fill(HIST("hYMCRec"), YD0(candidate1)); registry.fill(HIST("hSelectionStatusMCRec"), candidate1.isSelD0barTOFplusRICHPID() + (candidate1.isSelD0TOFplusRICHPID() * 2)); } - //fill invariant mass plots from D0/D0bar signal and background candidates - if (candidate1.isSelD0TOFplusRICHPID() >= selectionFlagD0) { //only reco as D0 - if (candidate1.flagMCMatchRec() == 1 << DecayType::D0ToPiK) { //also matched as D0 + // fill invariant mass plots from D0/D0bar signal and background candidates + if (candidate1.isSelD0TOFplusRICHPID() >= selectionFlagD0) { // only reco as D0 + if (candidate1.flagMCMatchRec() == 1 << DecayType::D0ToPiK) { // also matched as D0 registry.fill(HIST("hMassD0MCRecSig"), InvMassD0(candidate1), candidate1.pt(), efficiencyWeight); } else if (candidate1.flagMCMatchRec() == -(1 << DecayType::D0ToPiK)) { registry.fill(HIST("hMassD0MCRecRefl"), InvMassD0(candidate1), candidate1.pt(), efficiencyWeight); @@ -311,8 +298,8 @@ struct HFCorrelatorD0D0barBarrelFullPIDMcRec { registry.fill(HIST("hMassD0MCRecBkg"), InvMassD0(candidate1), candidate1.pt(), efficiencyWeight); } } - if (candidate1.isSelD0barTOFplusRICHPID() >= selectionFlagD0bar) { //only reco as D0bar - if (candidate1.flagMCMatchRec() == -(1 << DecayType::D0ToPiK)) { //also matched as D0bar + if (candidate1.isSelD0barTOFplusRICHPID() >= selectionFlagD0bar) { // only reco as D0bar + if (candidate1.flagMCMatchRec() == -(1 << DecayType::D0ToPiK)) { // also matched as D0bar registry.fill(HIST("hMassD0barMCRecSig"), InvMassD0bar(candidate1), candidate1.pt(), efficiencyWeight); } else if (candidate1.flagMCMatchRec() == 1 << DecayType::D0ToPiK) { registry.fill(HIST("hMassD0barMCRecRefl"), InvMassD0bar(candidate1), candidate1.pt(), efficiencyWeight); @@ -321,34 +308,34 @@ struct HFCorrelatorD0D0barBarrelFullPIDMcRec { } } - //D-Dbar correlation dedicated section - //if the candidate is selected ad D0, search for D0bar and evaluate correlations - if (candidate1.isSelD0TOFplusRICHPID() < selectionFlagD0) { //discard candidates not selected as D0 in outer loop + // D-Dbar correlation dedicated section + // if the candidate is selected ad D0, search for D0bar and evaluate correlations + if (candidate1.isSelD0TOFplusRICHPID() < selectionFlagD0) { // discard candidates not selected as D0 in outer loop continue; } - flagD0Signal = candidate1.flagMCMatchRec() == 1 << DecayType::D0ToPiK; //flagD0Signal 'true' if candidate1 matched to D0 (particle) - flagD0Reflection = candidate1.flagMCMatchRec() == -(1 << DecayType::D0ToPiK); //flagD0Reflection 'true' if candidate1, selected as D0 (particle), is matched to D0bar (antiparticle) + flagD0Signal = candidate1.flagMCMatchRec() == 1 << DecayType::D0ToPiK; // flagD0Signal 'true' if candidate1 matched to D0 (particle) + flagD0Reflection = candidate1.flagMCMatchRec() == -(1 << DecayType::D0ToPiK); // flagD0Reflection 'true' if candidate1, selected as D0 (particle), is matched to D0bar (antiparticle) for (auto& candidate2 : candidates) { - if (!(candidate2.hfflag() & 1 << DecayType::D0ToPiK)) { //check decay channel flag for candidate2 + if (!(candidate2.hfflag() & 1 << DecayType::D0ToPiK)) { // check decay channel flag for candidate2 continue; } - if (candidate2.isSelD0barTOFplusRICHPID() < selectionFlagD0bar) { //discard candidates not selected as D0bar in inner loop + if (candidate2.isSelD0barTOFplusRICHPID() < selectionFlagD0bar) { // discard candidates not selected as D0bar in inner loop continue; } - flagD0barSignal = candidate2.flagMCMatchRec() == -(1 << DecayType::D0ToPiK); //flagD0barSignal 'true' if candidate2 matched to D0bar (antiparticle) - flagD0barReflection = candidate2.flagMCMatchRec() == 1 << DecayType::D0ToPiK; //flagD0barReflection 'true' if candidate2, selected as D0bar (antiparticle), is matched to D0 (particle) + flagD0barSignal = candidate2.flagMCMatchRec() == -(1 << DecayType::D0ToPiK); // flagD0barSignal 'true' if candidate2 matched to D0bar (antiparticle) + flagD0barReflection = candidate2.flagMCMatchRec() == 1 << DecayType::D0ToPiK; // flagD0barReflection 'true' if candidate2, selected as D0bar (antiparticle), is matched to D0 (particle) if (cutYCandMax >= 0. && std::abs(YD0(candidate2)) > cutYCandMax) { continue; } if (cutPtCandMin >= 0. && candidate2.pt() < cutPtCandMin) { continue; } - //Excluding trigger self-correlations (possible in case of both mass hypotheses accepted) + // 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/ref, 2 = bkg/sig, 3 = ref/bkg, 4 = ref/ref, 5 = ref/sig, 6 = sig/bkg, 7 = sig/ref, 8 = sig/sig + // 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; } @@ -370,10 +357,10 @@ struct HFCorrelatorD0D0barBarrelFullPIDMcRec { pairSignalStatus); double etaCut = 0.; double ptCut = 0.; - do { //fill pairs vs etaCut plot + do { // fill pairs vs etaCut plot ptCut = 0.; etaCut += incrementEtaCut; - do { //fill pairs vs etaCut plot + 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); } @@ -382,42 +369,19 @@ struct HFCorrelatorD0D0barBarrelFullPIDMcRec { } while (etaCut < maxEtaCut - epsilon); } // end inner loop (Dbars) - } //end outer loop + } // end outer loop } -}; -/// D0-D0bar correlation pair builder - for MC gen-level analysis (no filter/selection, only true signal) -struct HFCorrelatorD0D0barBarrelFullPIDMcGen { + PROCESS_SWITCH(HfCorrelatorD0D0barBarrelFullPid, processMcRec, "Process MC Reco mode", true); - Produces entryD0D0barPair; - - HistogramRegistry registry{ - "registry", - {{"hMCEvtCount", "Event counter - MC gen;;entries", {HistType::kTH1F, {{1, -0.5, 0.5}}}}, - {"hPtCandMCGen", "D0,D0bar particles - MC gen;particle #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, - {"hEtaMCGen", "D0,D0bar particles - MC gen;particle #it{#eta};entries", {HistType::kTH1F, {{100, -5., 5.}}}}, - {"hPhiMCGen", "D0,D0bar particles - MC gen;particle #it{#varphi};entries", {HistType::kTH1F, {{32, 0., 2. * o2::constants::math::PI}}}}, - {"hYMCGen", "D0,D0bar candidates - MC gen;candidate #it{#y};entries", {HistType::kTH1F, {{100, -5., 5.}}}}, - {"hcountD0D0barPerEvent", "D0,D0bar particles - MC gen;Number per event;entries", {HistType::kTH1F, {{20, 0., 20.}}}}, - {"hDDbarVsEtaCut", "D0,D0bar pairs vs #eta cut of D mesons;#eta_{max};entries", {HistType::kTH2F, {{(int)(maxEtaCut / incrementEtaCut), 0., maxEtaCut}, {(int)(ptThresholdForMaxEtaCut / incrementPtThreshold), 0., ptThresholdForMaxEtaCut}}}}, - {"hDDbarVsDaughterEtaCut", "D0,D0bar pairs vs #eta cut on D daughters;#eta_{max};entries", {HistType::kTH2F, {{(int)(maxEtaCut / incrementEtaCut), 0., maxEtaCut}, {(int)(ptThresholdForMaxEtaCut / incrementPtThreshold), 0., ptThresholdForMaxEtaCut}}}}}}; - - Configurable cutYCandMax{"cutYCandMax", -1., "max. cand. rapidity"}; - Configurable cutPtCandMin{"cutPtCandMin", -1., "min. cand. pT"}; - Configurable> bins{"ptBinsForMass", std::vector{o2::analysis::hf_cuts_d0_topik::pTBins_v}, "pT bin limits for trigger counters"}; - - void init(o2::framework::InitContext&) - { - registry.add("hcountD0triggersMCGen", "D0 trigger particles - MC gen;;N of trigger D0", {HistType::kTH2F, {{1, -0.5, 0.5}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); - } - - void process(aod::McCollision const& mccollision, MCParticlesPlus const& particlesMC) + /// D0-D0bar correlation pair builder - for MC gen-level analysis (no filter/selection, only true signal) + void processMcGen(aod::McCollision const& mccollision, MCParticlesPlus const& particlesMC) { int counterD0D0bar = 0; registry.fill(HIST("hMCEvtCount"), 0); - //MC gen level + // 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! + // check if the particle is D0 or D0bar (for general plot filling and selection, so both cases are fine) - NOTE: decay channel is not probed! if (std::abs(particle1.pdgCode()) != pdg::Code::kD0) { continue; } @@ -434,14 +398,14 @@ struct HFCorrelatorD0D0barBarrelFullPIDMcGen { registry.fill(HIST("hYMCGen"), yD); counterD0D0bar++; - //D-Dbar correlation dedicated section - //if it's a D0 particle, search for D0bar and evaluate correlations - if (particle1.pdgCode() != pdg::Code::kD0) { //just checking the particle PDG, not the decay channel (differently from Reco: you have a BR factor btw such levels!) + // D-Dbar correlation dedicated section + // if it's a D0 particle, search for D0bar and evaluate correlations + if (particle1.pdgCode() != pdg::Code::kD0) { // just checking the particle PDG, not the decay channel (differently from Reco: you have a BR factor btw such levels!) continue; } - registry.fill(HIST("hcountD0triggersMCGen"), 0, particle1.pt()); //to count trigger D0 (for normalisation) + registry.fill(HIST("hCountD0triggersMCGen"), 0, particle1.pt()); // to count trigger D0 (for normalisation) for (auto& particle2 : particlesMC) { - if (particle2.pdgCode() != pdg::Code::kD0bar) { //check that inner particle is D0bar + if (particle2.pdgCode() != pdg::Code::kD0bar) { // check that inner particle is D0bar continue; } if (cutYCandMax >= 0. && std::abs(RecoDecay::Y(array{particle2.px(), particle2.py(), particle2.pz()}, RecoDecay::getMassPDG(particle2.pdgCode()))) > cutYCandMax) { @@ -454,10 +418,13 @@ struct HFCorrelatorD0D0barBarrelFullPIDMcGen { particle2.eta() - particle1.eta(), particle1.pt(), particle2.pt()); + entryD0D0barRecoInfo(1.864, + 1.864, + 8); // dummy information double etaCut = 0.; double ptCut = 0.; - //fill pairs vs etaCut plot + // fill pairs vs etaCut plot bool rightDecayChannels = false; if ((std::abs(particle1.flagMCMatchGen()) == 1 << DecayType::D0ToPiK) && (std::abs(particle2.flagMCMatchGen()) == 1 << DecayType::D0ToPiK)) { rightDecayChannels = true; @@ -465,11 +432,11 @@ struct HFCorrelatorD0D0barBarrelFullPIDMcGen { do { ptCut = 0.; etaCut += incrementEtaCut; - do { //fill pairs vs etaCut plot - if (std::abs(particle1.eta()) < etaCut && std::abs(particle2.eta()) < etaCut && particle1.pt() > ptCut && particle2.pt() > ptCut) { //fill with D and Dbar acceptance checks + do { // fill pairs vs etaCut plot + if (std::abs(particle1.eta()) < etaCut && std::abs(particle2.eta()) < etaCut && particle1.pt() > ptCut && particle2.pt() > ptCut) { // fill with D and Dbar acceptance checks registry.fill(HIST("hDDbarVsEtaCut"), etaCut - epsilon, ptCut + epsilon); } - if (rightDecayChannels) { //fill with D and Dbar daughter particls acceptance checks + if (rightDecayChannels) { // fill with D and Dbar daughter particls acceptance checks double etaCandidate1Daughter1 = particle1.daughter0_as().eta(); double etaCandidate1Daughter2 = particle1.daughter1_as().eta(); double etaCandidate2Daughter1 = particle2.daughter0_as().eta(); @@ -483,325 +450,30 @@ struct HFCorrelatorD0D0barBarrelFullPIDMcGen { ptCut += incrementPtThreshold; } while (ptCut < ptThresholdForMaxEtaCut - epsilon); } while (etaCut < maxEtaCut - epsilon); - } //end inner loop - } //end outer loop - registry.fill(HIST("hcountD0D0barPerEvent"), counterD0D0bar); - } -}; - -/// D0-D0bar correlation pair builder - LIKE SIGN - for real data and data-like analysis (i.e. reco-level w/o matching request via MC truth) -/// NOTE: At the moment, both dPhi-symmetrical correlation pairs (part1-part2 and part2-part1) are filled, -/// since we bin in pT and selecting as trigger the largest pT particle would bias the distributions w.r.t. the ULS case. -struct HFCorrelatorD0D0barBarrelFullPIDLs { - - Produces entryD0D0barPair; - Produces entryD0D0barRecoInfo; - - HistogramRegistry registry{ - "registry", - //NOTE: use hMassD0 for trigger normalisation (S*0.955), and hMass2DCorrelationPairs (in final task) for 2D-sideband-subtraction purposes - {{"hPtCand", "D0,D0bar candidates;candidate #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, - {"hPtProng0", "D0,D0bar candidates;prong 0 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, - {"hPtProng1", "D0,D0bar candidates;prong 1 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, - {"hSelectionStatus", "D0,D0bar candidates;selection status;entries", {HistType::kTH1F, {{4, -0.5, 3.5}}}}, - {"hEta", "D0,D0bar candidates;candidate #it{#eta};entries", {HistType::kTH1F, {{100, -5., 5.}}}}, - {"hPhi", "D0,D0bar candidates;candidate #it{#varphi};entries", {HistType::kTH1F, {{32, 0., 2. * o2::constants::math::PI}}}}, - {"hY", "D0,D0bar candidates;candidate #it{#y};entries", {HistType::kTH1F, {{100, -5., 5.}}}}}}; - - Configurable selectionFlagD0{"selectionFlagD0", 1, "Selection Flag for D0"}; - Configurable selectionFlagD0bar{"selectionFlagD0bar", 1, "Selection Flag for D0bar"}; - Configurable cutYCandMax{"cutYCandMax", -1., "max. cand. rapidity"}; - Configurable cutPtCandMin{"cutPtCandMin", -1., "min. cand. pT"}; - Configurable> bins{"ptBinsForMass", std::vector{o2::analysis::hf_cuts_d0_topik::pTBins_v}, "pT bin limits for candidate mass plots"}; - - void init(o2::framework::InitContext&) - { - registry.add("hMass", "D0,D0bar candidates;inv. mass (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{120, 1.5848, 2.1848}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); - registry.add("hMassD0", "D0,D0bar candidates;inv. mass D0 only (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{120, 1.5848, 2.1848}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); - registry.add("hMassD0bar", "D0,D0bar candidates;inv. mass D0bar only (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{120, 1.5848, 2.1848}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); - } - - Filter filterSelectCandidates = (aod::hf_selcandidate_d0_ALICE3_Barrel::isSelD0TOFplusRICHPID >= selectionFlagD0 || aod::hf_selcandidate_d0_ALICE3_Barrel::isSelD0barTOFplusRICHPID >= selectionFlagD0bar); - - void process(aod::Collision const& collision, soa::Filtered> const& candidates) - { - for (auto& candidate1 : candidates) { - //check decay channel flag for candidate1 - if (!(candidate1.hfflag() & 1 << DecayType::D0ToPiK)) { - continue; - } - if (cutYCandMax >= 0. && std::abs(YD0(candidate1)) > cutYCandMax) { - continue; - } - if (cutPtCandMin >= 0. && candidate1.pt() < cutPtCandMin) { - continue; - } - //fill invariant mass plots and generic info from all D0/D0bar candidates - if (candidate1.isSelD0TOFplusRICHPID() >= selectionFlagD0) { - registry.fill(HIST("hMass"), InvMassD0(candidate1), candidate1.pt()); - registry.fill(HIST("hMassD0"), InvMassD0(candidate1), candidate1.pt()); - } - if (candidate1.isSelD0barTOFplusRICHPID() >= selectionFlagD0bar) { - registry.fill(HIST("hMass"), InvMassD0bar(candidate1), candidate1.pt()); - registry.fill(HIST("hMassD0bar"), InvMassD0bar(candidate1), candidate1.pt()); - } - registry.fill(HIST("hPtCand"), candidate1.pt()); - registry.fill(HIST("hPtProng0"), candidate1.ptProng0()); - registry.fill(HIST("hPtProng1"), candidate1.ptProng1()); - registry.fill(HIST("hEta"), candidate1.eta()); - registry.fill(HIST("hPhi"), candidate1.phi()); - registry.fill(HIST("hY"), YD0(candidate1)); - registry.fill(HIST("hSelectionStatus"), candidate1.isSelD0barTOFplusRICHPID() + (candidate1.isSelD0TOFplusRICHPID() * 2)); - - //D-Dbar correlation dedicated section - //For like-sign, first loop on both D0 and D0bars. First candidate is for sure a D0 and D0bars (checked before, so don't re-check anything on it) - for (auto& candidate2 : candidates) { - //check decay channel flag for candidate2 - if (!(candidate2.hfflag() & 1 << DecayType::D0ToPiK)) { - continue; - } - //for the associated, has to have smaller pT, and pass D0sel if trigger passes D0sel, or D0barsel if trigger passes D0barsel - if ((candidate1.isSelD0TOFplusRICHPID() >= selectionFlagD0 && candidate2.isSelD0TOFplusRICHPID() >= selectionFlagD0) || (candidate1.isSelD0barTOFplusRICHPID() >= selectionFlagD0bar && candidate2.isSelD0barTOFplusRICHPID() >= selectionFlagD0bar)) { - if (cutYCandMax >= 0. && std::abs(YD0(candidate2)) > cutYCandMax) { - continue; - } - if (cutPtCandMin >= 0. && candidate2.pt() < cutPtCandMin) { - continue; - } - //Excluding self-correlations - if (candidate1.mRowIndex == candidate2.mRowIndex) { - continue; - } - entryD0D0barPair(getDeltaPhi(candidate2.phi(), candidate1.phi()), - candidate2.eta() - candidate1.eta(), - candidate1.pt(), - candidate2.pt()); - entryD0D0barRecoInfo(InvMassD0(candidate1), - InvMassD0bar(candidate2), - 0); - } - //note: candidates selected as both D0 and D0bar are used, and considered in both situation (but not auto-correlated): reflections could play a relevant role. - //another, more restrictive, option, could be to consider only candidates selected with a single option (D0 xor D0bar) - } // end inner loop (Dbars) - } //end outer loop - } -}; - -/// D0-D0bar correlation pair builder - LIKE SIGN - for MC reco analysis (data-like but matching to true DO and D0bar) -/// NOTE: At the moment, both dPhi-symmetrical correlation pairs (part1-part2 and part2-part1) are filled, -/// since we bin in pT and selecting as trigger the largest pT particle would bias the distributions w.r.t. the ULS case. -struct HFCorrelatorD0D0barBarrelFullPIDMcRecLs { - - Produces entryD0D0barPair; - Produces entryD0D0barRecoInfo; - - HistogramRegistry registry{ - "registry", - //NOTE: use hMassD0 for trigger normalisation (S*0.955), and hMass2DCorrelationPairs (in final task) for 2D-sideband-subtraction purposes - {{"hPtCandMCRec", "D0,D0bar candidates - MC reco;candidate #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, - {"hPtProng0MCRec", "D0,D0bar candidates - MC reco;prong 0 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, - {"hPtProng1MCRec", "D0,D0bar candidates - MC reco;prong 1 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, - {"hSelectionStatusMCRec", "D0,D0bar candidates - MC reco;selection status;entries", {HistType::kTH1F, {{4, -0.5, 3.5}}}}, - {"hEtaMCRec", "D0,D0bar candidates - MC reco;candidate #it{#eta};entries", {HistType::kTH1F, {{100, -5., 5.}}}}, - {"hPhiMCRec", "D0,D0bar candidates - MC reco;candidate #it{#varphi};entries", {HistType::kTH1F, {{32, 0., 2. * o2::constants::math::PI}}}}, - {"hYMCRec", "D0,D0bar candidates - MC reco;candidate #it{#y};entries", {HistType::kTH1F, {{100, -5., 5.}}}}}}; - - Configurable selectionFlagD0{"selectionFlagD0", 1, "Selection Flag for D0"}; - Configurable selectionFlagD0bar{"selectionFlagD0bar", 1, "Selection Flag for D0bar"}; - Configurable cutYCandMax{"cutYCandMax", -1., "max. cand. rapidity"}; - Configurable cutPtCandMin{"cutPtCandMin", -1., "min. cand. pT"}; - Configurable> bins{"ptBinsForMass", std::vector{o2::analysis::hf_cuts_d0_topik::pTBins_v}, "pT bin limits for candidate mass plots"}; - - void init(o2::framework::InitContext&) - { - registry.add("hMassD0MCRec", "D0,D0bar candidates - MC reco;inv. mass (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{120, 1.5848, 2.1848}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); - registry.add("hMassD0barMCRec", "D0,D0bar candidates - MC reco;inv. mass D0 only (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{120, 1.5848, 2.1848}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); - } - - Filter filterSelectCandidates = (aod::hf_selcandidate_d0_ALICE3_Barrel::isSelD0TOFplusRICHPID >= selectionFlagD0 || aod::hf_selcandidate_d0_ALICE3_Barrel::isSelD0barTOFplusRICHPID >= selectionFlagD0bar); - - void process(aod::Collision const& collision, soa::Filtered> const& candidates) - { - //MC reco level - for (auto& candidate1 : candidates) { - //check decay channel flag for candidate1 - if (!(candidate1.hfflag() & 1 << DecayType::D0ToPiK)) { - continue; - } - if (cutYCandMax >= 0. && std::abs(YD0(candidate1)) > cutYCandMax) { - continue; - } - if (cutPtCandMin >= 0. && candidate1.pt() < cutPtCandMin) { - continue; - } - if (std::abs(candidate1.flagMCMatchRec()) == 1 << DecayType::D0ToPiK) { - //fill invariant mass plots and generic info from all D0/D0bar candidates - if (candidate1.isSelD0TOFplusRICHPID() >= selectionFlagD0 && candidate1.flagMCMatchRec() == DecayType::D0ToPiK) { //only reco and matched as D0 - registry.fill(HIST("hMassD0MCRec"), InvMassD0(candidate1)); - } - if (candidate1.isSelD0barTOFplusRICHPID() >= selectionFlagD0bar && candidate1.flagMCMatchRec() == DecayType::D0ToPiK) { //only reco and matched as D0bar - registry.fill(HIST("hMassD0barMCRec"), InvMassD0bar(candidate1)); - } - registry.fill(HIST("hPtCandMCRec"), candidate1.pt()); - registry.fill(HIST("hPtProng0MCRec"), candidate1.ptProng0()); - registry.fill(HIST("hPtProng1MCRec"), candidate1.ptProng1()); - registry.fill(HIST("hEtaMCRec"), candidate1.eta()); - registry.fill(HIST("hPhiMCRec"), candidate1.phi()); - registry.fill(HIST("hYMCRec"), YD0(candidate1)); - registry.fill(HIST("hSelectionStatusMCRec"), candidate1.isSelD0barTOFplusRICHPID() + (candidate1.isSelD0TOFplusRICHPID() * 2)); - - //D-Dbar correlation dedicated section - //For like-sign, first loop on both D0 and D0bars. First candidate is for sure a D0 and D0bars (looping on filtered) and was already matched, so don't re-check anything on it) - for (auto& candidate2 : candidates) { - //check decay channel flag for candidate2 - if (!(candidate2.hfflag() & 1 << DecayType::D0ToPiK)) { - continue; - } - bool conditionLSForD0 = (candidate1.isSelD0TOFplusRICHPID() >= selectionFlagD0bar && candidate1.flagMCMatchRec() == 1 << DecayType::D0ToPiK) && (candidate2.isSelD0TOFplusRICHPID() >= selectionFlagD0bar && candidate2.flagMCMatchRec() == 1 << DecayType::D0ToPiK); - bool conditionLSForD0bar = (candidate1.isSelD0barTOFplusRICHPID() >= selectionFlagD0bar && candidate1.flagMCMatchRec() == -(1 << DecayType::D0ToPiK)) && (candidate2.isSelD0barTOFplusRICHPID() >= selectionFlagD0bar && candidate2.flagMCMatchRec() == -(1 << DecayType::D0ToPiK)); - if (conditionLSForD0 || conditionLSForD0bar) { //LS pair (of D0 or of D0bar) + pt2= 0. && std::abs(YD0(candidate2)) > cutYCandMax) { - continue; - } - if (cutPtCandMin >= 0. && candidate2.pt() < cutPtCandMin) { - continue; - } - //Excluding self-correlations - if (candidate1.mRowIndex == candidate2.mRowIndex) { - continue; - } - entryD0D0barPair(getDeltaPhi(candidate2.phi(), candidate1.phi()), - candidate2.eta() - candidate1.eta(), - candidate1.pt(), - candidate2.pt()); - entryD0D0barRecoInfo(InvMassD0(candidate1), - InvMassD0bar(candidate2), - 0); //for LS studies we set a dummy 0 for pairSignalStatus (there are no more the usual 4 possible combinations) - - } //end inner if (MC match) - } // end inner loop (Dbars) - } //end outer if (MC match) - } //end outer loop + } // end inner loop + } // end outer loop + registry.fill(HIST("hCountD0D0barPerEvent"), counterD0D0bar); } -}; - -/// D0-D0bar correlation pair builder - for MC gen-level analysis, like sign particles -/// NOTE: At the moment, both dPhi-symmetrical correlation pairs (part1-part2 and part2-part1) are filled, -/// since we bin in pT and selecting as trigger the largest pT particle would bias the distributions w.r.t. the ULS case. -struct HFCorrelatorD0D0barBarrelFullPIDMcGenLs { - 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})"}}}); - } + PROCESS_SWITCH(HfCorrelatorD0D0barBarrelFullPid, processMcGen, "Process MC Gen mode", false); - void process(aod::McCollision const& mccollision, MCParticlesPlus const& particlesMC) + /// c-cbar correlator table builder - for MC gen-level analysis + void processccbar(aod::McCollision const& mccollision, MCParticlesPlus const& particlesMC) { - int counterD0D0bar = 0; registry.fill(HIST("hMCEvtCount"), 0); - //MC gen level - for (auto& particle1 : particlesMC) { - //check if the particle is D0 or D0bar (both can be trigger) - NOTE: decay channel is not probed! - if (std::abs(particle1.pdgCode()) != pdg::Code::kD0) { - continue; - } - double yD = RecoDecay::Y(array{particle1.px(), particle1.py(), particle1.pz()}, RecoDecay::getMassPDG(particle1.pdgCode())); - if (cutYCandMax >= 0. && std::abs(yD) > cutYCandMax) { - continue; - } - if (cutPtCandMin >= 0. && particle1.pt() < cutPtCandMin) { - continue; - } - - registry.fill(HIST("hPtCandMCGen"), particle1.pt()); - registry.fill(HIST("hEtaMCGen"), particle1.eta()); - registry.fill(HIST("hPhiMCGen"), particle1.phi()); - registry.fill(HIST("hYMCGen"), yD); - counterD0D0bar++; - //D-Dbar correlation dedicated section - //if it's D0, search for D0bar and evaluate correlations. - registry.fill(HIST("hcountD0triggersMCGen"), 0, particle1.pt()); //to count trigger D0 (normalisation) - for (auto& particle2 : particlesMC) { - if (std::abs(particle2.pdgCode()) != pdg::Code::kD0) { //check that associated is a D0/D0bar (both are fine) - continue; - } - if (cutYCandMax >= 0. && std::abs(RecoDecay::Y(array{particle2.px(), particle2.py(), particle2.pz()}, RecoDecay::getMassPDG(particle2.pdgCode()))) > cutYCandMax) { - continue; - } - if (cutPtCandMin >= 0. && particle2.pt() < cutPtCandMin) { - continue; - } - if (particle2.pdgCode() == particle1.pdgCode()) { //like-sign condition (both 421 or both -421) and pT_Trig>pT_assoc - //Excluding self-correlations - if (particle1.mRowIndex == particle2.mRowIndex) { - continue; - } - entryD0D0barPair(getDeltaPhi(particle2.phi(), particle1.phi()), - particle2.eta() - particle1.eta(), - particle1.pt(), - particle2.pt()); - } - } // end inner loop (Dbars) - } //end outer loop - registry.fill(HIST("hcountD0D0barPerEvent"), counterD0D0bar); - } -}; - -/// c-cbar correlator table builder - for MC gen-level analysis -struct HfCorrelatorCCbarBarrelFullPIDMCGen { - - Produces entryccbarPair; + int counterCCbar = 0, counterCCbarBeforeEtasel = 0; - HistogramRegistry registry{ - "registry", - {{"hMCEvtCount", "Event counter - MC gen;;entries", {HistType::kTH1F, {{1, -0.5, 0.5}}}}, - {"hPtCandMCGen", "c,cbar particles - MC gen;particle #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, - {"hEtaMCGen", "c,cbar particles - MC gen;particle #it{#eta};entries", {HistType::kTH1F, {{100, -5., 5.}}}}, - {"hYMCGen", "c,cbar candidates - MC gen;candidate #it{#y};entries", {HistType::kTH1F, {{100, -5., 5.}}}}, - {"hPhiMCGen", "c,cbar particles - MC gen;particle #it{#varphi};entries", {HistType::kTH1F, {{32, 0., 2. * o2::constants::math::PI}}}}, - {"hcountCCbarPerEvent", "c,cbar particles - MC gen;Number per event;entries", {HistType::kTH1F, {{20, 0., 20.}}}}, - {"hcountCCbarPerEventPreEtaCut", "c,cbar particles - MC gen;Number per event pre #eta cut;entries", {HistType::kTH1F, {{20, 0., 20.}}}}}}; - - Configurable cutYCandMax{"cutYCandMax", -1., "max. cand. rapidity"}; - Configurable cutPtCandMin{"cutPtCandMin", -1., "min. cand. pT"}; - Configurable> bins{"ptBinsForMass", std::vector{o2::analysis::hf_cuts_d0_topik::pTBins_v}, "pT bin limits for trigger counters"}; - - void init(o2::framework::InitContext&) - { - registry.add("hcountCtriggersMCGen", "c trigger particles - MC gen;;N of trigger c quark", {HistType::kTH2F, {{1, -0.5, 0.5}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); - } - - void process(aod::McCollision const& mccollision, MCParticlesPlus const& particlesMC) - { - registry.fill(HIST("hMCEvtCount"), 0); - int counterccbar = 0, counterccbarPreEtasel = 0; - - //loop over particles at MC gen level + // loop over particles at MC gen level for (auto& particle1 : particlesMC) { - if (std::abs(particle1.pdgCode()) != PDG_t::kCharm) { //search c or cbar particles + if (std::abs(particle1.pdgCode()) != PDG_t::kCharm) { // search c or cbar particles continue; } int partMothPDG = particle1.mother0_as().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 + // 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) + counterCCbarBeforeEtasel++; // count c or cbar (before kinematic selection) double yC = RecoDecay::Y(array{particle1.px(), particle1.py(), particle1.pz()}, RecoDecay::getMassPDG(particle1.pdgCode())); if (cutYCandMax >= 0. && std::abs(yC) > cutYCandMax) { continue; @@ -813,17 +485,17 @@ struct HfCorrelatorCCbarBarrelFullPIDMCGen { registry.fill(HIST("hEtaMCGen"), particle1.eta()); registry.fill(HIST("hPhiMCGen"), particle1.phi()); registry.fill(HIST("hYMCGen"), yC); - counterccbar++; //count if c or cbar don't come from themselves during fragmentation (after kinematic selection) + 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. + // c-cbar correlation dedicated section + // if it's c, search for cbar and evaluate correlations. if (particle1.pdgCode() != PDG_t::kCharm) { continue; } - registry.fill(HIST("hcountCtriggersMCGen"), 0, particle1.pt()); //to count trigger c quark (for normalisation) + registry.fill(HIST("hCountCtriggersMCGen"), 0, particle1.pt()); // to count trigger c quark (for normalisation) for (auto& particle2 : particlesMC) { - if (particle2.pdgCode() != PDG_t::kCharmBar) { //check that inner particle is a cbar + if (particle2.pdgCode() != PDG_t::kCharmBar) { // check that inner particle is a cbar continue; } if (cutYCandMax >= 0. && std::abs(RecoDecay::Y(array{particle2.px(), particle2.py(), particle2.pz()}, RecoDecay::getMassPDG(particle2.pdgCode()))) > cutYCandMax) { @@ -832,136 +504,27 @@ struct HfCorrelatorCCbarBarrelFullPIDMCGen { if (cutPtCandMin >= 0. && particle2.pt() < cutPtCandMin) { continue; } - //check whether mothers of quark cbar (from associated loop) are still '-4' particles - in that case the cbar quark comes from its own fragmentation, skip it + // 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 (particle2.mother0_as().pdgCode() == PDG_t::kCharmBar) { continue; } - entryccbarPair(getDeltaPhi(particle2.phi(), particle1.phi()), - particle2.eta() - particle1.eta(), - particle1.pt(), - particle2.pt()); - } // end inner loop - } //end outer loop - registry.fill(HIST("hcountCCbarPerEvent"), counterccbar); - registry.fill(HIST("hcountCCbarPerEventPreEtaCut"), counterccbarPreEtasel); - } -}; - -/// c-cbar correlator table builder - for MC gen-level analysis - Like Sign -struct HfCorrelatorCCbarBarrelFullPIDMCGenLs { - - Produces entryccbarPair; - - HistogramRegistry registry{ - "registry", - {{"hMCEvtCount", "Event counter - MC gen;;entries", {HistType::kTH1F, {{1, -0.5, 0.5}}}}, - {"hPtCandMCGen", "c,cbar particles - MC gen;particle #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, - {"hEtaMCGen", "c,cbar particles - MC gen;particle #it{#eta};entries", {HistType::kTH1F, {{100, -5., 5.}}}}, - {"hYMCGen", "c,cbar candidates - MC gen;candidate #it{#y};entries", {HistType::kTH1F, {{100, -5., 5.}}}}, - {"hPhiMCGen", "c,cbar particles - MC gen;particle #it{#varphi};entries", {HistType::kTH1F, {{32, 0., 2. * o2::constants::math::PI}}}}, - {"hcountCCbarPerEvent", "c,cbar particles - MC gen;Number per event;entries", {HistType::kTH1F, {{20, 0., 20.}}}}, - {"hcountCCbarPerEventPreEtaCut", "c,cbar particles - MC gen;Number per event pre #eta cut;entries", {HistType::kTH1F, {{20, 0., 20.}}}}}}; - - Configurable cutYCandMax{"cutYCandMax", -1., "max. cand. rapidity"}; - Configurable cutPtCandMin{"cutPtCandMin", -1., "min. cand. pT"}; - Configurable> bins{"ptBinsForMass", std::vector{o2::analysis::hf_cuts_d0_topik::pTBins_v}, "pT bin limits for trigger counters"}; - - void init(o2::framework::InitContext&) - { - registry.add("hcountCtriggersMCGen", "c trigger particles - MC gen;;N of trigger c quark", {HistType::kTH2F, {{1, -0.5, 0.5}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); - } - - void process(aod::McCollision const& mccollision, MCParticlesPlus const& particlesMC) - { - registry.fill(HIST("hMCEvtCount"), 0); - int counterccbar = 0, counterccbarPreEtasel = 0; - - //loop over particles at MC gen level - for (auto& particle1 : particlesMC) { - if (std::abs(particle1.pdgCode()) != PDG_t::kCharm) { //search c or cbar particles - continue; - } - int partMothPDG = particle1.mother0_as().pdgCode(); - //check whether mothers of quark c/cbar are still '4'/'-4' particles - in that case the c/cbar quark comes from its own fragmentation, skip it - if (partMothPDG == particle1.pdgCode()) { - continue; - } - counterccbarPreEtasel++; //count c or cbar (before kinematic selection) - double yC = RecoDecay::Y(array{particle1.px(), particle1.py(), particle1.pz()}, RecoDecay::getMassPDG(particle1.pdgCode())); - if (cutYCandMax >= 0. && std::abs(yC) > cutYCandMax) { - continue; - } - if (cutPtCandMin >= 0. && particle1.pt() < cutPtCandMin) { - continue; - } - registry.fill(HIST("hPtCandMCGen"), particle1.pt()); - registry.fill(HIST("hEtaMCGen"), particle1.eta()); - registry.fill(HIST("hPhiMCGen"), particle1.phi()); - registry.fill(HIST("hYMCGen"), yC); - counterccbar++; //count if c or cbar don't come from themselves during fragmentation (after kinematic selection) - - //c-cbar correlation dedicated section - registry.fill(HIST("hcountCtriggersMCGen"), 0, particle1.pt()); //to count trigger c quark (for normalisation) - - for (auto& particle2 : particlesMC) { - if (std::abs(particle2.pdgCode()) != PDG_t::kCharm) { //search c or cbar for associated particles - continue; - } - if (cutYCandMax >= 0. && std::abs(RecoDecay::Y(array{particle2.px(), particle2.py(), particle2.pz()}, RecoDecay::getMassPDG(particle2.pdgCode()))) > cutYCandMax) { - continue; - } - if (cutPtCandMin >= 0. && particle2.pt() < cutPtCandMin) { - continue; - } - if (particle2.pdgCode() == particle1.pdgCode()) { - //check whether mothers of quark cbar (from associated loop) are still '-4' particles - in that case the cbar quark comes from its own fragmentation, skip it - if (particle2.mother0_as().pdgCode() == particle2.pdgCode()) { - continue; - } - //Excluding self-correlations - if (particle1.mRowIndex == particle2.mRowIndex) { - continue; - } - entryccbarPair(getDeltaPhi(particle2.phi(), particle1.phi()), + 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); + entryD0D0barRecoInfo(1.864, + 1.864, + 8); // dummy information + } // end inner loop + } // end outer loop + registry.fill(HIST("hCountCCbarPerEvent"), counterCCbar); + registry.fill(HIST("hCountCCbarPerEventBeforeEtaCut"), counterCCbarBeforeEtasel); } + + PROCESS_SWITCH(HfCorrelatorD0D0barBarrelFullPid, processccbar, "Process ccbar pairs", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { - WorkflowSpec workflow{}; - const bool doMCccbar = cfgc.options().get("doMCccbar"); - const bool doMCGen = cfgc.options().get("doMCGen"); - const bool doMCRec = cfgc.options().get("doMCRec"); - const bool doLikeSign = cfgc.options().get("doLikeSign"); - if (!doLikeSign) { //unlike-sign analyses - if (doMCGen) { //MC-Gen analysis - workflow.push_back(adaptAnalysisTask(cfgc)); - } else if (doMCRec) { //MC-Reco analysis - workflow.push_back(adaptAnalysisTask(cfgc)); - } else if (doMCccbar) { //MC-Reco analysis - workflow.push_back(adaptAnalysisTask(cfgc)); - } else { //data analysis - workflow.push_back(adaptAnalysisTask(cfgc)); - } - } else { //like-sign analyses - if (doMCGen) { //MC-Gen analysis - workflow.push_back(adaptAnalysisTask(cfgc)); - } else if (doMCRec) { //MC-Reco analysis - workflow.push_back(adaptAnalysisTask(cfgc)); - } else if (doMCccbar) { //MC-Reco analysis - workflow.push_back(adaptAnalysisTask(cfgc)); - } else { //data analysis - workflow.push_back(adaptAnalysisTask(cfgc)); - } - } - - return workflow; + return WorkflowSpec{adaptAnalysisTask(cfgc)}; } diff --git a/PWGHF/Tasks/HFCorrelatorDplusDminus.cxx b/PWGHF/Tasks/HFCorrelatorDplusDminus.cxx index 7401e3318fc..dcb4c87c7cc 100644 --- a/PWGHF/Tasks/HFCorrelatorDplusDminus.cxx +++ b/PWGHF/Tasks/HFCorrelatorDplusDminus.cxx @@ -9,7 +9,7 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. -/// \file HFCorrelatorDplusDminus.cxx +/// \file HfCorrelatorDplusDminus.cxx /// \brief Dplus-Dminus correlator task - data-like, MC-reco and MC-kine analyses. For ULS and LS pairs /// /// \author Fabio Colamaria , INFN Bari @@ -29,18 +29,6 @@ using namespace o2::aod::hf_correlation_ddbar; using namespace o2::analysis::hf_cuts_dplus_topikpi; using namespace o2::constants::math; -void customize(std::vector& workflowOptions) -{ - ConfigParamSpec optionDoLikeSign{"doLikeSign", VariantType::Bool, false, {"Run Like-Sign analysis."}}; - ConfigParamSpec optionDoMCccbar{"doMCccbar", VariantType::Bool, false, {"Run MC-Gen dedicated tasks."}}; - ConfigParamSpec optionDoMCGen{"doMCGen", VariantType::Bool, false, {"Run MC-Gen dedicated tasks."}}; - ConfigParamSpec optionDoMCRec{"doMCRec", VariantType::Bool, true, {"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" /// @@ -62,25 +50,58 @@ const int npTBinsMassAndEfficiency = o2::analysis::hf_cuts_dplus_topikpi::npTBin const double efficiencyDmesonDefault[npTBinsMassAndEfficiency] = {}; auto efficiencyDmeson_v = std::vector{efficiencyDmesonDefault, efficiencyDmesonDefault + npTBinsMassAndEfficiency}; -/// Dplus-Dminus correlation pair builder - for real data and data-like analysis (i.e. reco-level w/o matching request via MC truth) +// 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 = 2. * o2::constants::math::PI; +const int yAxisBins = 100; +const double yAxisMin = -5.; +const double yAxisMax = 5.; +const int ptDAxisBins = 180; +const double ptDAxisMin = 0.; +const double ptDAxisMax = 36.; + +using MCParticlesPlus2Prong = soa::Join; +using MCParticlesPlus3Prong = soa::Join; + struct HfCorrelatorDplusDminus { Produces entryDplusDminusPair; Produces entryDplusDminusRecoInfo; HistogramRegistry registry{ "registry", - //NOTE: use hMassDplus for trigger normalisation (S*0.955), and hMass2DCorrelationPairs (in final task) for 2D-sideband-subtraction purposes - {{"hPtCand", "Dplus,Dminus candidates;candidate #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, - {"hPtProng0", "Dplus,Dminus candidates;prong 0 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, - {"hPtProng1", "Dplus,Dminus candidates;prong 1 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, - {"hPtProng2", "Dplus,Dminus candidates;prong 2 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, + // NOTE: use hMassDplus for trigger normalisation (S*0.955), and hMass2DCorrelationPairs (in final task) for 2D-sideband-subtraction purposes + {{"hPtCand", "Dplus,Dminus candidates;candidate #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{ptDAxisBins, ptDAxisMin, ptDAxisMax}}}}, + {"hPtProng0", "Dplus,Dminus candidates;prong 0 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{ptDAxisBins, ptDAxisMin, ptDAxisMax}}}}, + {"hPtProng1", "Dplus,Dminus candidates;prong 1 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{ptDAxisBins, ptDAxisMin, ptDAxisMax}}}}, + {"hPtProng2", "Dplus,Dminus candidates;prong 2 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{ptDAxisBins, ptDAxisMin, ptDAxisMax}}}}, {"hSelectionStatus", "Dplus,Dminus candidates;selection status;entries", {HistType::kTH1F, {{2, -0.5, 1.5}}}}, - {"hEta", "Dplus,Dminus candidates;candidate #it{#eta};entries", {HistType::kTH1F, {{100, -5., 5.}}}}, - {"hPhi", "Dplus,Dminus candidates;candidate #it{#varphi};entries", {HistType::kTH1F, {{32, 0., 2. * o2::constants::math::PI}}}}, - {"hY", "Dplus,Dminus candidates;candidate #it{#y};entries", {HistType::kTH1F, {{100, -5., 5.}}}}, + {"hEta", "Dplus,Dminus candidates;candidate #it{#eta};entries", {HistType::kTH1F, {{yAxisBins, yAxisMin, yAxisMax}}}}, + {"hPhi", "Dplus,Dminus candidates;candidate #it{#varphi};entries", {HistType::kTH1F, {{phiAxisBins, phiAxisMin, phiAxisMax}}}}, + {"hY", "Dplus,Dminus candidates;candidate #it{y};entries", {HistType::kTH1F, {{yAxisBins, yAxisMin, yAxisMax}}}}, {"hMultiplicityPreSelection", "multiplicity prior to selection;multiplicity;entries", {HistType::kTH1F, {{10000, 0., 10000.}}}}, {"hMultiplicity", "multiplicity;multiplicity;entries", {HistType::kTH1F, {{10000, 0., 10000.}}}}, - {"hDDbarVsEtaCut", "Dplus,Dminus pairs vs #eta cut;#eta_{max};entries", {HistType::kTH2F, {{(int)(maxEtaCut / incrementEtaCut), 0., maxEtaCut}, {(int)(ptThresholdForMaxEtaCut / incrementPtThreshold), 0., ptThresholdForMaxEtaCut}}}}}}; + {"hDDbarVsEtaCut", "Dplus,Dminus pairs vs #eta cut;#eta_{max};candidates #it{p}_{T} threshold (GeV/#it{c});entries", {HistType::kTH2F, {{(int)(maxEtaCut / incrementEtaCut), 0., maxEtaCut}, {(int)(ptThresholdForMaxEtaCut / incrementPtThreshold), 0., ptThresholdForMaxEtaCut}}}}, + {"hPtCandMCRec", "Dplus,Dminus candidates - MC reco;candidate #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{ptDAxisBins, ptDAxisMin, ptDAxisMax}}}}, + {"hPtProng0MCRec", "Dplus,Dminus candidates - MC reco;prong 0 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{ptDAxisBins, ptDAxisMin, ptDAxisMax}}}}, + {"hPtProng1MCRec", "Dplus,Dminus candidates - MC reco;prong 1 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{ptDAxisBins, ptDAxisMin, ptDAxisMax}}}}, + {"hPtProng2MCRec", "Dplus,Dminus candidates - MC reco;prong 2 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{ptDAxisBins, ptDAxisMin, ptDAxisMax}}}}, + {"hSelectionStatusMCRec", "Dplus,Dminus candidates - MC reco;selection status;entries", {HistType::kTH1F, {{2, -0.5, 1.5}}}}, + {"hEtaMCRec", "Dplus,Dminus candidates - MC reco;candidate #it{#eta};entries", {HistType::kTH1F, {{yAxisBins, yAxisMin, yAxisMax}}}}, + {"hPhiMCRec", "Dplus,Dminus candidates - MC reco;candidate #it{#varphi};entries", {HistType::kTH1F, {{phiAxisBins, phiAxisMin, phiAxisMax}}}}, + {"hYMCRec", "Dplus,Dminus candidates - MC reco;candidate #it{y};entries", {HistType::kTH1F, {{yAxisBins, yAxisMin, yAxisMax}}}}, + {"hMCEvtCount", "Event counter - MC gen;;entries", {HistType::kTH1F, {{1, -0.5, 0.5}}}}, + {"hPtCandMCGen", "Dplus,Dminus particles - MC gen;particle #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{ptDAxisBins, ptDAxisMin, ptDAxisMax}}}}, + {"hEtaMCGen", "Dplus,Dminus particles - MC gen;particle #it{#eta};entries", {HistType::kTH1F, {{yAxisBins, yAxisMin, yAxisMax}}}}, + {"hPhiMCGen", "Dplus,Dminus particles - MC gen;particle #it{#varphi};entries", {HistType::kTH1F, {{phiAxisBins, phiAxisMin, phiAxisMax}}}}, + {"hYMCGen", "Dplus,Dminus candidates - MC gen;candidate #it{y};entries", {HistType::kTH1F, {{yAxisBins, yAxisMin, yAxisMax}}}}, + {"hCountDplusDminusPerEvent", "Dplus,Dminus particles - MC gen;Number per event;entries", {HistType::kTH1F, {{20, 0., 20.}}}}, + {"hDDbarVsDaughterEtaCut", "Dplus,Dminus pairs vs #eta cut on D daughters;#eta_{max};candidates #it{p}_{T} threshold (GeV/#it{c});entries", {HistType::kTH2F, {{(int)(maxEtaCut / incrementEtaCut), 0., maxEtaCut}, {(int)(ptThresholdForMaxEtaCut / incrementPtThreshold), 0., ptThresholdForMaxEtaCut}}}}, + {"hCountCCbarPerEvent", "c,cbar particles - MC gen;Number per event;entries", {HistType::kTH1F, {{20, 0., 20.}}}}, + {"hCountCCbarPerEventBeforeEtaCut", "c,cbar particles - MC gen;Number per event pre #eta cut;entries", {HistType::kTH1F, {{20, 0., 20.}}}}}}; Configurable selectionFlagDplus{"selectionFlagDplus", 1, "Selection Flag for Dplus,Dminus"}; Configurable cutYCandMax{"cutYCandMax", -1., "max. cand. rapidity"}; @@ -93,14 +114,22 @@ struct HfCorrelatorDplusDminus { void init(o2::framework::InitContext&) { - registry.add("hMass", "Dplus,Dminus candidates;inv. mass (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{120, 1.5848, 2.1848}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); - registry.add("hMassDplus", "Dplus,Dminus candidates;inv. mass Dplus only (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{120, 1.5848, 2.1848}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); - registry.add("hMassDminus", "Dplus,Dminus candidates;inv. mass Dminus only (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{120, 1.5848, 2.1848}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + auto vbins = (std::vector)bins; + registry.add("hMass", "Dplus,Dminus candidates;inv. mass (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{massAxisBins, massAxisMin, massAxisMax}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hMassDplus", "Dplus,Dminus candidates;inv. mass Dplus only (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{massAxisBins, massAxisMin, massAxisMax}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hMassDminus", "Dplus,Dminus candidates;inv. mass Dminus only (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{massAxisBins, massAxisMin, massAxisMax}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hMassDplusMCRecSig", "Dplus signal candidates - MC reco;inv. mass D+ only (#pi K #pi) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{massAxisBins, massAxisMin, massAxisMax}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hMassDminusMCRecSig", "Dminus signal candidates - MC reco;inv. mass D- only (#pi K #pi) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{massAxisBins, massAxisMin, massAxisMax}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hMassDplusMCRecBkg", "Dplus background candidates - MC reco;inv. mass D+ only (#pi K #pi) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{massAxisBins, massAxisMin, massAxisMax}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hMassDminusMCRecBkg", "Dminus background candidates - MC reco;inv. mass D- only (#pi K #pi) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{massAxisBins, massAxisMin, massAxisMax}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hCountDplustriggersMCGen", "Dplus trigger particles - MC gen;;N of trigger D0", {HistType::kTH2F, {{1, -0.5, 0.5}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hCountCtriggersMCGen", "c trigger particles - MC gen;;N of trigger c quark", {HistType::kTH2F, {{1, -0.5, 0.5}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); } Filter filterSelectCandidates = (aod::hf_selcandidate_dplus::isSelDplusToPiKPi >= selectionFlagDplus); - void process(aod::Collision const& collision, soa::Join& tracks, soa::Filtered> const& candidates, aod::BigTracks const& bigtracks) + /// Dplus-Dminus 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& tracks, soa::Filtered> const& candidates, aod::BigTracks const& bigtracks) { int nTracks = 0; if (collision.numContrib() > 1) { @@ -127,8 +156,8 @@ struct HfCorrelatorDplusDminus { if (cutPtCandMin >= 0. && candidate1.pt() < cutPtCandMin) { continue; } - //check decay channel flag for candidate1 - if (!(candidate1.hfflag() & 1 << DecayType::DPlusToPiKPi)) { //probably dummy since already selected? not sure... + // check decay channel flag for candidate1 + if (!(candidate1.hfflag() & 1 << DecayType::DPlusToPiKPi)) { // probably dummy since already selected? not sure... continue; } @@ -137,13 +166,13 @@ struct HfCorrelatorDplusDminus { efficiencyWeight = 1. / efficiencyDmeson->at(o2::analysis::findBin(bins, candidate1.pt())); } - int outerParticleSign = 1; //Dplus + int outerParticleSign = 1; // Dplus auto outerSecondTrack = candidate1.index1_as(); if (outerSecondTrack.sign() == 1) { - outerParticleSign = -1; //Dminus (second daughter track is positive) + outerParticleSign = -1; // Dminus (second daughter track is positive) } - //fill invariant mass plots and generic info from all Dplus/Dminus candidates + // fill invariant mass plots and generic info from all Dplus/Dminus candidates if (outerParticleSign == 1) { registry.fill(HIST("hMass"), InvMassDPlus(candidate1), candidate1.pt(), efficiencyWeight); registry.fill(HIST("hMassDplus"), InvMassDPlus(candidate1), candidate1.pt(), efficiencyWeight); @@ -160,18 +189,18 @@ struct HfCorrelatorDplusDminus { registry.fill(HIST("hY"), YDPlus(candidate1)); registry.fill(HIST("hSelectionStatus"), candidate1.isSelDplusToPiKPi()); - //D-Dbar correlation dedicated section - //if the candidate is a Dplus, search for Dminus and evaluate correlations + // D-Dbar correlation dedicated section + // if the candidate is a Dplus, search for Dminus and evaluate correlations if (outerParticleSign != 1) { continue; } for (auto& candidate2 : candidates) { - //check decay channel flag for candidate2 - if (!(candidate2.hfflag() & 1 << DecayType::DPlusToPiKPi)) { //probably dummy since already selected? not sure... + // check decay channel flag for candidate2 + if (!(candidate2.hfflag() & 1 << DecayType::DPlusToPiKPi)) { // probably dummy since already selected? not sure... continue; } auto innerSecondTrack = candidate2.index1_as(); - if (innerSecondTrack.sign() != 1) { //keep only Dminus (with second daughter track positive) + if (innerSecondTrack.sign() != 1) { // keep only Dminus (with second daughter track positive) continue; } if (cutYCandMax >= 0. && std::abs(YDPlus(candidate2)) > cutYCandMax) { @@ -189,10 +218,10 @@ struct HfCorrelatorDplusDminus { 0); double etaCut = 0.; double ptCut = 0.; - do { //fill pairs vs etaCut plot + do { // fill pairs vs etaCut plot ptCut = 0.; etaCut += incrementEtaCut; - do { //fill pairs vs etaCut plot + 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); } @@ -200,50 +229,13 @@ struct HfCorrelatorDplusDminus { } while (ptCut < ptThresholdForMaxEtaCut - epsilon); } while (etaCut < maxEtaCut - epsilon); } // end inner loop (Dminus) - } //end outer loop (Dplus) - } -}; - -/// Dplus-Dminus correlation pair builder - for MC reco-level analysis (candidates matched to true signal only, but also the various bkg sources are studied) -struct HfCorrelatorDplusDminusMcRec { - Produces entryDplusDminusPair; - Produces entryDplusDminusRecoInfo; - - HistogramRegistry registry{ - "registry", - //NOTE: use hMassDplus for trigger normalisation (S*0.955), and hMass2DCorrelationPairs (in final task) for 2D-sideband-subtraction purposes - {{"hPtCandMCRec", "Dplus,Dminus candidates - MC reco;candidate #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, - {"hPtProng0MCRec", "Dplus,Dminus candidates - MC reco;prong 0 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, - {"hPtProng1MCRec", "Dplus,Dminus candidates - MC reco;prong 1 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, - {"hPtProng2MCRec", "Dplus,Dminus candidates - MC reco;prong 2 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, - {"hSelectionStatusMCRec", "Dplus,Dminus candidates - MC reco;selection status;entries", {HistType::kTH1F, {{2, -0.5, 1.5}}}}, - {"hEtaMCRec", "Dplus,Dminus candidates - MC reco;candidate #it{#eta};entries", {HistType::kTH1F, {{100, -5., 5.}}}}, - {"hPhiMCRec", "Dplus,Dminus candidates - MC reco;candidate #it{#varphi};entries", {HistType::kTH1F, {{32, 0., 2. * o2::constants::math::PI}}}}, - {"hYMCRec", "Dplus,Dminus candidates - MC reco;candidate #it{#y};entries", {HistType::kTH1F, {{100, -5., 5.}}}}, - {"hMultiplicityPreSelection", "multiplicity prior to selection;multiplicity;entries", {HistType::kTH1F, {{10000, 0., 10000.}}}}, - {"hMultiplicity", "multiplicity;multiplicity;entries", {HistType::kTH1F, {{10000, 0., 10000.}}}}, - {"hDDbarVsEtaCut", "Dplus,Dminus pairs vs #eta cut;#eta_{max};entries", {HistType::kTH2F, {{(int)(maxEtaCut / incrementEtaCut), 0., maxEtaCut}, {(int)(ptThresholdForMaxEtaCut / incrementPtThreshold), 0., ptThresholdForMaxEtaCut}}}}}}; - - Configurable selectionFlagDplus{"selectionFlagDplus", 1, "Selection Flag for Dplus,Dminus"}; - Configurable cutYCandMax{"cutYCandMax", -1., "max. cand. rapidity"}; - Configurable cutPtCandMin{"cutPtCandMin", -1., "min. cand. pT"}; - Configurable> bins{"ptBinsForMassAndEfficiency", std::vector{o2::analysis::hf_cuts_dplus_topikpi::pTBins_v}, "pT bin limits for candidate mass plots and efficiency"}; - Configurable> efficiencyDmeson{"efficiencyDmeson", std::vector{efficiencyDmeson_v}, "Efficiency values for D0 meson"}; - Configurable flagApplyEfficiency{"efficiencyFlagD", 1, "Flag for applying D-meson efficiency weights"}; - Configurable multMin{"multMin", 0., "minimum multiplicity accepted"}; - Configurable multMax{"multMax", 10000., "maximum multiplicity accepted"}; - - void init(o2::framework::InitContext&) - { - registry.add("hMassDplusMCRecSig", "Dplus signal candidates - MC reco;inv. mass D+ only (#pi K #pi) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{120, 1.5848, 2.1848}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); - registry.add("hMassDminusMCRecSig", "Dminus signal candidates - MC reco;inv. mass D- only (#pi K #pi) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{120, 1.5848, 2.1848}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); - registry.add("hMassDplusMCRecBkg", "Dplus background candidates - MC reco;inv. mass D+ only (#pi K #pi) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{120, 1.5848, 2.1848}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); - registry.add("hMassDminusMCRecBkg", "Dminus background candidates - MC reco;inv. mass D- only (#pi K #pi) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{120, 1.5848, 2.1848}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + } // end outer loop (Dplus) } - Filter filterSelectCandidates = (aod::hf_selcandidate_dplus::isSelDplusToPiKPi >= selectionFlagDplus); + PROCESS_SWITCH(HfCorrelatorDplusDminus, processData, "Process data", false); - void process(aod::Collision const& collision, soa::Join& tracks, soa::Filtered> const& candidates, aod::BigTracks const& bigtracks) + /// Dplus-Dminus 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& tracks, soa::Filtered> const& candidates, aod::BigTracks const& bigtracks) { int nTracks = 0; if (collision.numContrib() > 1) { @@ -263,11 +255,11 @@ struct HfCorrelatorDplusDminusMcRec { } registry.fill(HIST("hMultiplicity"), nTracks); - //MC reco level + // MC reco level bool flagDplusSignal = false; bool flagDminusSignal = false; for (auto& candidate1 : candidates) { - //check decay channel flag for candidate1 + // check decay channel flag for candidate1 if (!(candidate1.hfflag() & 1 << DecayType::DPlusToPiKPi)) { continue; } @@ -283,16 +275,16 @@ struct HfCorrelatorDplusDminusMcRec { efficiencyWeight = 1. / efficiencyDmeson->at(o2::analysis::findBin(bins, candidate1.pt())); } - int outerParticleSign = 1; //Dplus + int outerParticleSign = 1; // Dplus auto outerSecondTrack = candidate1.index1_as(); if (outerSecondTrack.sign() == 1) { - outerParticleSign = -1; //Dminus (second daughter track is positive) + outerParticleSign = -1; // Dminus (second daughter track is positive) } if (std::abs(candidate1.flagMCMatchRec()) == 1 << DecayType::DPlusToPiKPi) { - //fill invariant mass plots and per-candidate distributions from Dplus/Dminus signal candidates - if (outerParticleSign == 1) { //reco and matched as Dplus + // fill invariant mass plots and per-candidate distributions from Dplus/Dminus signal candidates + if (outerParticleSign == 1) { // reco and matched as Dplus registry.fill(HIST("hMassDplusMCRecSig"), InvMassDPlus(candidate1), candidate1.pt(), efficiencyWeight); - } else { //reco and matched as Dminus + } else { // reco and matched as Dminus registry.fill(HIST("hMassDminusMCRecSig"), InvMassDPlus(candidate1), candidate1.pt(), efficiencyWeight); } registry.fill(HIST("hPtCandMCRec"), candidate1.pt()); @@ -304,36 +296,36 @@ struct HfCorrelatorDplusDminusMcRec { registry.fill(HIST("hYMCRec"), YDPlus(candidate1)); registry.fill(HIST("hSelectionStatusMCRec"), candidate1.isSelDplusToPiKPi()); } else { - //fill invariant mass plots from Dplus/Dminus background candidates - if (outerParticleSign == 1) { //reco as Dplus + // fill invariant mass plots from Dplus/Dminus background candidates + if (outerParticleSign == 1) { // reco as Dplus registry.fill(HIST("hMassDplusMCRecBkg"), InvMassDPlus(candidate1), candidate1.pt(), efficiencyWeight); - } else { //matched as Dminus + } else { // matched as Dminus registry.fill(HIST("hMassDminusMCRecBkg"), InvMassDPlus(candidate1), candidate1.pt(), efficiencyWeight); } } - //D-Dbar correlation dedicated section + // D-Dbar correlation dedicated section if (outerParticleSign == -1) { - continue; //reject Dminus in outer loop + continue; // reject Dminus in outer loop } - flagDplusSignal = std::abs(candidate1.flagMCMatchRec()) == 1 << DecayType::DPlusToPiKPi; //flagDplusSignal 'true' if candidate1 matched to Dplus + flagDplusSignal = std::abs(candidate1.flagMCMatchRec()) == 1 << DecayType::DPlusToPiKPi; // flagDplusSignal 'true' if candidate1 matched to Dplus for (auto& candidate2 : candidates) { - if (!(candidate2.hfflag() & 1 << DecayType::DPlusToPiKPi)) { //check decay channel flag for candidate2 + if (!(candidate2.hfflag() & 1 << DecayType::DPlusToPiKPi)) { // check decay channel flag for candidate2 continue; } auto innerSecondTrack = candidate2.index1_as(); - if (innerSecondTrack.sign() != 1) { //keep only Dminus (with second daughter track positive) + if (innerSecondTrack.sign() != 1) { // keep only Dminus (with second daughter track positive) continue; } - flagDminusSignal = std::abs(candidate2.flagMCMatchRec()) == 1 << DecayType::DPlusToPiKPi; //flagDminusSignal 'true' if candidate2 matched to Dminus + flagDminusSignal = std::abs(candidate2.flagMCMatchRec()) == 1 << DecayType::DPlusToPiKPi; // flagDminusSignal 'true' if candidate2 matched to Dminus if (cutYCandMax >= 0. && std::abs(YDPlus(candidate2)) > cutYCandMax) { continue; } if (cutPtCandMin >= 0. && candidate2.pt() < cutPtCandMin) { continue; } - //choice of options (Dplus/Dminus signal/bkg) - int pairSignalStatus = 0; //0 = bkg/bkg, 1 = bkg/ref, 2 = bkg/sig, 3 = ref/bkg, 4 = ref/ref, 5 = ref/sig, 6 = sig/bkg, 7 = sig/ref, 8 = sig/sig. Of course only 0,2,6,8 are relevant for D+D- + // choice of options (Dplus/Dminus 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. Of course only 0,2,6,8 are relevant for D+D- if (flagDplusSignal) { pairSignalStatus += 6; } @@ -349,10 +341,10 @@ struct HfCorrelatorDplusDminusMcRec { pairSignalStatus); double etaCut = 0.; double ptCut = 0.; - do { //fill pairs vs etaCut plot + do { // fill pairs vs etaCut plot ptCut = 0.; etaCut += incrementEtaCut; - do { //fill pairs vs etaCut plot + 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); } @@ -361,44 +353,19 @@ struct HfCorrelatorDplusDminusMcRec { } while (etaCut < maxEtaCut - epsilon); } // end inner loop (Dbars) - } //end outer loop + } // end outer loop } -}; - -using MCParticlesPlus = soa::Join; -/// Dplus-Dminus correlation pair builder - for MC gen-level analysis (no filter/selection, only true signal) -struct HfCorrelatorDplusDminusMcGen { - - Produces entryDplusDminusPair; + PROCESS_SWITCH(HfCorrelatorDplusDminus, processMcRec, "Process MC Reco mode", true); - HistogramRegistry registry{ - "registry", - {{"hMCEvtCount", "Event counter - MC gen;;entries", {HistType::kTH1F, {{1, -0.5, 0.5}}}}, - {"hPtCandMCGen", "Dplus,Dminus particles - MC gen;particle #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, - {"hEtaMCGen", "Dplus,Dminus particles - MC gen;particle #it{#eta};entries", {HistType::kTH1F, {{100, -5., 5.}}}}, - {"hPhiMCGen", "Dplus,Dminus particles - MC gen;particle #it{#varphi};entries", {HistType::kTH1F, {{32, 0., 2. * o2::constants::math::PI}}}}, - {"hYMCGen", "Dplus,Dminus candidates - MC gen;candidate #it{#y};entries", {HistType::kTH1F, {{100, -5., 5.}}}}, - {"hcountDplusDminusPerEvent", "Dplus,Dminus particles - MC gen;Number per event;entries", {HistType::kTH1F, {{20, 0., 20.}}}}, - {"hDDbarVsEtaCut", "Dplus,Dminus pairs vs #eta cut of D mesons;#eta_{max};entries", {HistType::kTH2F, {{(int)(maxEtaCut / incrementEtaCut), 0., maxEtaCut}, {(int)(ptThresholdForMaxEtaCut / incrementPtThreshold), 0., ptThresholdForMaxEtaCut}}}}, - {"hDDbarVsDaughterEtaCut", "Dplus,Dminus pairs vs #eta cut on D daughters;#eta_{max};entries", {HistType::kTH2F, {{(int)(maxEtaCut / incrementEtaCut), 0., maxEtaCut}, {(int)(ptThresholdForMaxEtaCut / incrementPtThreshold), 0., ptThresholdForMaxEtaCut}}}}}}; - - Configurable cutYCandMax{"cutYCandMax", -1., "max. cand. rapidity"}; - Configurable cutPtCandMin{"cutPtCandMin", -1., "min. cand. pT"}; - Configurable> bins{"ptBinsForMass", std::vector{o2::analysis::hf_cuts_dplus_topikpi::pTBins_v}, "pT bin limits for candidate mass plots"}; - - void init(o2::framework::InitContext&) - { - registry.add("hcountDplustriggersMCGen", "Dplus trigger particles - MC gen;;N of trigger D0", {HistType::kTH2F, {{1, -0.5, 0.5}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); - } - - void process(aod::McCollision const& mccollision, MCParticlesPlus const& particlesMC) + /// Dplus-Dminus correlation pair builder - for MC gen-level analysis (no filter/selection, only true signal) + void processMcGen(aod::McCollision const& mccollision, MCParticlesPlus3Prong const& particlesMC) { int counterDplusDminus = 0; registry.fill(HIST("hMCEvtCount"), 0); - //MC gen level + // MC gen level for (auto& particle1 : particlesMC) { - //check if the particle is Dplus or Dminus (for general plot filling and selection, so both cases are fine) - NOTE: decay channel is not probed! + // check if the particle is Dplus or Dminus (for general plot filling and selection, so both cases are fine) - NOTE: decay channel is not probed! if (std::abs(particle1.pdgCode()) != pdg::Code::kDPlus) { continue; } @@ -415,14 +382,14 @@ struct HfCorrelatorDplusDminusMcGen { registry.fill(HIST("hYMCGen"), yD); counterDplusDminus++; - //D-Dbar correlation dedicated section - //if it's a Dplus particle, search for Dminus and evaluate correlations - if (particle1.pdgCode() != pdg::Code::kDPlus) { //just checking the particle PDG, not the decay channel (differently from Reco: you have a BR factor btw such levels!) + // D-Dbar correlation dedicated section + // if it's a Dplus particle, search for Dminus and evaluate correlations + if (particle1.pdgCode() != pdg::Code::kDPlus) { // just checking the particle PDG, not the decay channel (differently from Reco: you have a BR factor btw such levels!) continue; } - registry.fill(HIST("hcountDplustriggersMCGen"), 0, particle1.pt()); //to count trigger Dplus (for normalisation) + registry.fill(HIST("hCountDplustriggersMCGen"), 0, particle1.pt()); // to count trigger Dplus (for normalisation) for (auto& particle2 : particlesMC) { - if (particle2.pdgCode() != -pdg::Code::kDPlus) { //check that inner particle is a Dminus + if (particle2.pdgCode() != -pdg::Code::kDPlus) { // check that inner particle is a Dminus continue; } if (cutYCandMax >= 0. && std::abs(RecoDecay::Y(array{particle2.px(), particle2.py(), particle2.pz()}, RecoDecay::getMassPDG(particle2.pdgCode()))) > cutYCandMax) { @@ -435,10 +402,13 @@ struct HfCorrelatorDplusDminusMcGen { particle2.eta() - particle1.eta(), particle1.pt(), particle2.pt()); + entryDplusDminusRecoInfo(1.869, + 1.869, + 8); // Dummy double etaCut = 0.; double ptCut = 0.; - //fill pairs vs etaCut plot + // fill pairs vs etaCut plot bool rightDecayChannels = false; if ((std::abs(particle1.flagMCMatchGen()) == 1 << DecayType::DPlusToPiKPi) && (std::abs(particle2.flagMCMatchGen()) == 1 << DecayType::DPlusToPiKPi)) { rightDecayChannels = true; @@ -446,17 +416,17 @@ struct HfCorrelatorDplusDminusMcGen { do { ptCut = 0.; etaCut += incrementEtaCut; - do { //fill pairs vs etaCut plot + do { // fill pairs vs etaCut plot if (std::abs(particle1.eta()) < etaCut && std::abs(particle2.eta()) < etaCut && particle1.pt() > ptCut && particle2.pt() > ptCut) { registry.fill(HIST("hDDbarVsEtaCut"), etaCut - epsilon, ptCut + epsilon); } - if (rightDecayChannels) { //fill with D and Dbar daughter particls acceptance checks - double etaCandidate1Daughter1 = particle1.daughter0_as().eta(); - double etaCandidate1Daughter2 = (particle1.daughter0_as() + 1).eta(); - double etaCandidate1Daughter3 = (particle1.daughter0_as() + 2).eta(); - double etaCandidate2Daughter1 = particle2.daughter0_as().eta(); - double etaCandidate2Daughter2 = (particle2.daughter0_as() + 1).eta(); - double etaCandidate2Daughter3 = (particle2.daughter0_as() + 2).eta(); + if (rightDecayChannels) { // fill with D and Dbar daughter particls acceptance checks + double etaCandidate1Daughter1 = particle1.daughter0_as().eta(); + double etaCandidate1Daughter2 = (particle1.daughter0_as() + 1).eta(); + double etaCandidate1Daughter3 = (particle1.daughter0_as() + 2).eta(); + double etaCandidate2Daughter1 = particle2.daughter0_as().eta(); + double etaCandidate2Daughter2 = (particle2.daughter0_as() + 1).eta(); + double etaCandidate2Daughter3 = (particle2.daughter0_as() + 2).eta(); if (std::abs(etaCandidate1Daughter1) < etaCut && std::abs(etaCandidate1Daughter2) < etaCut && std::abs(etaCandidate1Daughter3) < etaCut && std::abs(etaCandidate2Daughter1) < etaCut && std::abs(etaCandidate2Daughter2) < etaCut && std::abs(etaCandidate2Daughter3) < etaCut && particle1.pt() > ptCut && particle2.pt() > ptCut) { @@ -466,345 +436,30 @@ struct HfCorrelatorDplusDminusMcGen { ptCut += incrementPtThreshold; } while (ptCut < ptThresholdForMaxEtaCut - epsilon); } while (etaCut < maxEtaCut - epsilon); - } //end inner loop - } //end outer loop - registry.fill(HIST("hcountDplusDminusPerEvent"), counterDplusDminus); - } -}; - -/// Dplus-Dminus correlation pair builder - LIKE SIGN - for real data and data-like analysis (i.e. reco-level w/o matching request via MC truth) -/// NOTE: At the moment, both dPhi-symmetrical correlation pairs (part1-part2 and part2-part1) are filled, since we bin in pT and selecting as trigger the largest pT particle would bias the distributions w.r.t. the ULS case. -struct HfCorrelatorDplusDminusLs { - - Produces entryDplusDminusPair; - Produces entryDplusDminusRecoInfo; - - HistogramRegistry registry{ - "registry", - //NOTE: use hMassDplus for trigger normalisation (S*0.955), and hMass2DCorrelationPairs (in final task) for 2D-sideband-subtraction purposes - {{"hPtCand", "Dplus,Dminus candidates;candidate #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, - {"hPtProng0", "Dplus,Dminus candidates;prong 0 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, - {"hPtProng1", "Dplus,Dminus candidates;prong 1 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, - {"hPtProng2", "Dplus,Dminus candidates;prong 2 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, - {"hSelectionStatus", "Dplus,Dminus candidates;selection status;entries", {HistType::kTH1F, {{2, -0.5, 1.5}}}}, - {"hEta", "Dplus,Dminus candidates;candidate #it{#eta};entries", {HistType::kTH1F, {{100, -5., 5.}}}}, - {"hPhi", "Dplus,Dminus candidates;candidate #it{#varphi};entries", {HistType::kTH1F, {{32, 0., 2. * o2::constants::math::PI}}}}, - {"hY", "Dplus,Dminus candidates;candidate #it{#y};entries", {HistType::kTH1F, {{100, -5., 5.}}}}}}; - - Configurable selectionFlagDplus{"dSelectivonFlagDplus", 1, "Selection Flag for Dplus,Dminus"}; - Configurable cutYCandMax{"cutYCandMax", -1., "max. cand. rapidity"}; - Configurable cutPtCandMin{"cutPtCandMin", -1., "min. cand. pT"}; - Configurable> bins{"ptBinsForMass", std::vector{o2::analysis::hf_cuts_dplus_topikpi::pTBins_v}, "pT bin limits for candidate mass plots"}; - - void init(o2::framework::InitContext&) - { - registry.add("hMass", "Dplus,Dminus candidates;inv. mass (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{120, 1.5848, 2.1848}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); - registry.add("hMassDplus", "Dplus,Dminus candidates;inv. mass Dplus only (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{120, 1.5848, 2.1848}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); - registry.add("hMassDminus", "Dplus,Dminus candidates;inv. mass Dminus only (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{120, 1.5848, 2.1848}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); - } - - Filter filterSelectCandidates = (aod::hf_selcandidate_dplus::isSelDplusToPiKPi >= selectionFlagDplus); - - void process(aod::Collision const& collision, soa::Filtered> const& candidates, aod::BigTracks const& bigtracks) - { - for (auto& candidate1 : candidates) { - if (cutYCandMax >= 0. && std::abs(YDPlus(candidate1)) > cutYCandMax) { - continue; - } - if (cutPtCandMin >= 0. && candidate1.pt() < cutPtCandMin) { - continue; - } - //check decay channel flag for candidate1 - if (!(candidate1.hfflag() & 1 << DecayType::DPlusToPiKPi)) { //probably dummy since already selected? not sure... - continue; - } - int outerParticleSign = 1; //Dplus - auto outerSecondTrack = candidate1.index1_as(); - if (outerSecondTrack.sign() == 1) { - outerParticleSign = -1; //Dminus (second daughter track is positive) - } - //fill invariant mass plots and generic info from all Dplus/Dminus candidates - if (outerParticleSign == 1) { - registry.fill(HIST("hMass"), InvMassDPlus(candidate1), candidate1.pt()); - registry.fill(HIST("hMassDplus"), InvMassDPlus(candidate1), candidate1.pt()); - } else { - registry.fill(HIST("hMass"), InvMassDPlus(candidate1), candidate1.pt()); - registry.fill(HIST("hMassDminus"), InvMassDPlus(candidate1), candidate1.pt()); - } - registry.fill(HIST("hPtCand"), candidate1.pt()); - registry.fill(HIST("hPtProng0"), candidate1.ptProng0()); - registry.fill(HIST("hPtProng1"), candidate1.ptProng1()); - registry.fill(HIST("hPtProng2"), candidate1.ptProng2()); - registry.fill(HIST("hEta"), candidate1.eta()); - registry.fill(HIST("hPhi"), candidate1.phi()); - registry.fill(HIST("hY"), YDPlus(candidate1)); - registry.fill(HIST("hSelectionStatus"), candidate1.isSelDplusToPiKPi()); - - //D-Dbar correlation dedicated section - //For like-sign, first loop on both Dplus and Dminus. First candidate is for sure a Dplus/Dminus (checked before, so don't re-check anything on it) - for (auto& candidate2 : candidates) { - //check decay channel flag for candidate2 - if (!(candidate2.hfflag() & 1 << DecayType::DPlusToPiKPi)) { - continue; - } - //check if inner particle is Dplus or Dminus - int innerParticleSign = 1; //Dplus - auto innerSecondTrack = candidate2.index1_as(); - if (innerSecondTrack.sign() == 1) { - innerParticleSign = -1; //Dminus (second daughter track is positive) - } - //for the associated, has to have smaller pT, and has to be Dplus if outer is Dplu, or Dminus if outer is Dminus - if (innerParticleSign == outerParticleSign) { - if (cutYCandMax >= 0. && std::abs(YDPlus(candidate2)) > cutYCandMax) { - continue; - } - if (cutPtCandMin >= 0. && candidate2.pt() < cutPtCandMin) { - continue; - } - //Excluding self-correlations - if (candidate1.mRowIndex == candidate2.mRowIndex) { - continue; - } - entryDplusDminusPair(getDeltaPhi(candidate2.phi(), candidate1.phi()), - candidate2.eta() - candidate1.eta(), - candidate1.pt(), - candidate2.pt()); - entryDplusDminusRecoInfo(InvMassDPlus(candidate1), - InvMassDPlus(candidate2), - 0); - } - } // end inner loop - } //end outer loop - } -}; - -/// Dplus-Dminus correlation pair builder - LIKE SIGN - for MC reco analysis (data-like but matching to true Dplus and Dminus) -/// NOTE: At the moment, both dPhi-symmetrical correlation pairs (part1-part2 and part2-part1) are filled, since we bin in pT and selecting as trigger the largest pT particle would bias the distributions w.r.t. the ULS case. -struct HfCorrelatorDplusDminusMcRecLs { - Produces entryDplusDminusPair; - Produces entryDplusDminusRecoInfo; - - HistogramRegistry registry{ - "registry", - //NOTE: use hMassD0 for trigger normalisation (S*0.955), and hMass2DCorrelationPairs (in final task) for 2D-sideband-subtraction purposes - {{"hPtCandMCRec", "Dplus,Dminus candidates - MC reco;candidate #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, - {"hPtProng0MCRec", "Dplus,Dminus candidates - MC reco;prong 0 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, - {"hPtProng1MCRec", "Dplus,Dminus candidates - MC reco;prong 1 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, - {"hPtProng2MCRec", "Dplus,Dminus candidates - MC reco;prong 2 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, - {"hSelectionStatusMCRec", "Dplus,Dminus candidates - MC reco;selection status;entries", {HistType::kTH1F, {{2, -0.5, 1.5}}}}, - {"hEtaMCRec", "Dplus,Dminus candidates - MC reco;candidate #it{#eta};entries", {HistType::kTH1F, {{100, -5., 5.}}}}, - {"hPhiMCRec", "Dplus,Dminus candidates - MC reco;candidate #it{#varphi};entries", {HistType::kTH1F, {{32, 0., 2. * o2::constants::math::PI}}}}, - {"hYMCRec", "Dplus,Dminus candidates - MC reco;candidate #it{#y};entries", {HistType::kTH1F, {{100, -5., 5.}}}}}}; - - Configurable selectionFlagDplus{"selectionFlagDplus", 1, "Selection Flag for Dplus,Dminus"}; - Configurable cutYCandMax{"cutYCandMax", -1., "max. cand. rapidity"}; - Configurable cutPtCandMin{"cutPtCandMin", -1., "min. cand. pT"}; - Configurable> bins{"ptBinsForMass", std::vector{o2::analysis::hf_cuts_dplus_topikpi::pTBins_v}, "pT bin limits for candidate mass plots"}; - - void init(o2::framework::InitContext&) - { - registry.add("hMassDplusMCRec", "Dplus,Dminus candidates - MC reco;inv. mass (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{120, 1.5848, 2.1848}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); - registry.add("hMassDminusMCRec", "Dplus,Dminus candidates - MC reco;inv. mass D0 only (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{120, 1.5848, 2.1848}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); - } - - Filter filterSelectCandidates = (aod::hf_selcandidate_dplus::isSelDplusToPiKPi >= selectionFlagDplus); - - void process(aod::Collision const& collision, soa::Filtered> const& candidates, aod::BigTracks const& bigtracks) - { - //MC reco level - for (auto& candidate1 : candidates) { - //check decay channel flag for candidate1 - if (!(candidate1.hfflag() & 1 << DecayType::DPlusToPiKPi)) { - continue; - } - if (cutYCandMax >= 0. && std::abs(YDPlus(candidate1)) > cutYCandMax) { - continue; - } - if (cutPtCandMin >= 0. && candidate1.pt() < cutPtCandMin) { - continue; - } - int outerParticleSign = 1; //Dplus - auto outerSecondTrack = candidate1.index1_as(); - if (outerSecondTrack.sign() == 1) { - outerParticleSign = -1; //Dminus (second daughter track is positive) - } - if (std::abs(candidate1.flagMCMatchRec()) != 1 << DecayType::DPlusToPiKPi) { //keep only Dplus/Dminus matched candidates - continue; - } - //fill invariant mass plots and generic info from all Dplus/Dminus candidates - if (outerParticleSign == 1) { //matched as Dplus - registry.fill(HIST("hMassDplusMCRec"), InvMassDPlus(candidate1), candidate1.pt()); - } else { //matched as Dminus - registry.fill(HIST("hMassDminusMCRec"), InvMassDPlus(candidate1), candidate1.pt()); - } - registry.fill(HIST("hPtCandMCRec"), candidate1.pt()); - registry.fill(HIST("hPtProng0MCRec"), candidate1.ptProng0()); - registry.fill(HIST("hPtProng1MCRec"), candidate1.ptProng1()); - registry.fill(HIST("hPtProng2MCRec"), candidate1.ptProng2()); - registry.fill(HIST("hEtaMCRec"), candidate1.eta()); - registry.fill(HIST("hPhiMCRec"), candidate1.phi()); - registry.fill(HIST("hYMCRec"), YDPlus(candidate1)); - registry.fill(HIST("hSelectionStatusMCRec"), candidate1.isSelDplusToPiKPi()); - - //D-Dbar correlation dedicated section - //For like-sign, first loop on both Dplus and Dminus. First candidate is for sure a Dplus/Dminus (looping on filtered) and was already matched, so don't re-check anything on it) - for (auto& candidate2 : candidates) { - //check decay channel flag for candidate2 - if (!(candidate2.hfflag() & 1 << DecayType::DPlusToPiKPi)) { - continue; - } - int innerParticleSign = 1; //Dplus - auto innerSecondTrack = candidate2.index1_as(); - if (innerSecondTrack.sign() == 1) { - innerParticleSign = -1; //Dminus (second daughter track is positive) - } - if (std::abs(candidate2.flagMCMatchRec()) != 1 << DecayType::DPlusToPiKPi) { //reject fake candidates - continue; - } - if (outerParticleSign == innerParticleSign) { //LS pair (of Dplus or Dminus) + pt2= 0. && std::abs(YDPlus(candidate2)) > cutYCandMax) { - continue; - } - if (cutPtCandMin >= 0. && candidate2.pt() < cutPtCandMin) { - continue; - } - //Excluding self-correlations (in principle not possible due to the '<' condition, but could rounding break it?) - if (candidate1.mRowIndex == candidate2.mRowIndex) { - continue; - } - entryDplusDminusPair(getDeltaPhi(candidate2.phi(), candidate1.phi()), - candidate2.eta() - candidate1.eta(), - candidate1.pt(), - candidate2.pt()); - entryDplusDminusRecoInfo(InvMassDPlus(candidate1), - InvMassDPlus(candidate2), - 0); //for LS studies we set a dummy 0 for pairSignalStatus (there are no more the usual 4 possible combinations) - - } //end inner if (MC match) - } // end inner loop - } //end outer loop - } -}; - -/// Dplus-Dminus correlation pair builder - for MC gen-level analysis, like sign particles -/// NOTE: At the moment, both dPhi-symmetrical correlation pairs (part1-part2 and part2-part1) are filled, since we bin in pT and selecting as trigger the largest pT particle would bias the distributions w.r.t. the ULS case. -struct HfCorrelatorDplusDminusMcGenLs { - - Produces entryDplusDminusPair; - - HistogramRegistry registry{ - "registry", - {{"hMCEvtCount", "Event counter - MC gen;;entries", {HistType::kTH1F, {{1, -0.5, 0.5}}}}, - {"hPtCandMCGen", "Dplus,Dminus particles - MC gen;particle #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, - {"hEtaMCGen", "Dplus,Dminus particles - MC gen;particle #it{#eta};entries", {HistType::kTH1F, {{100, -5., 5.}}}}, - {"hPhiMCGen", "Dplus,Dminus particles - MC gen;particle #it{#varphi};entries", {HistType::kTH1F, {{32, 0., 2. * o2::constants::math::PI}}}}, - {"hYMCGen", "Dplus,Dminus candidates - MC gen;candidate #it{#y};entries", {HistType::kTH1F, {{100, -5., 5.}}}}, - {"hcountDplusDminusPerEvent", "Dplus,Dminus particles - MC gen;Number per event;entries", {HistType::kTH1F, {{20, 0., 20.}}}}}}; - - Configurable cutYCandMax{"cutYCandMax", -1., "max. cand. rapidity"}; - Configurable cutPtCandMin{"cutPtCandMin", -1., "min. cand. pT"}; - Configurable> bins{"ptBinsForMass", std::vector{o2::analysis::hf_cuts_dplus_topikpi::pTBins_v}, "pT bin limits for trigger counters"}; - - void init(o2::framework::InitContext&) - { - registry.add("hcountDplustriggersMCGen", "Dplus trigger particles - MC gen;;N of trigger D0", {HistType::kTH2F, {{1, -0.5, 0.5}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); - } - - void process(aod::McCollision const& mccollision, soa::Join const& particlesMC) - { - int counterDplusDminus = 0; - registry.fill(HIST("hMCEvtCount"), 0); - //MC gen level - for (auto& particle1 : particlesMC) { - //check if the particle is Dplus or Dminus (both can be trigger) - NOTE: decay channel is not probed! - if (std::abs(particle1.pdgCode()) != pdg::Code::kDPlus) { - continue; - } - double yD = RecoDecay::Y(array{particle1.px(), particle1.py(), particle1.pz()}, RecoDecay::getMassPDG(particle1.pdgCode())); - if (cutYCandMax >= 0. && std::abs(yD) > cutYCandMax) { - continue; - } - if (cutPtCandMin >= 0. && particle1.pt() < cutPtCandMin) { - continue; - } - - registry.fill(HIST("hPtCandMCGen"), particle1.pt()); - registry.fill(HIST("hEtaMCGen"), particle1.eta()); - registry.fill(HIST("hPhiMCGen"), particle1.phi()); - registry.fill(HIST("hYMCGen"), yD); - counterDplusDminus++; - - //D-Dbar correlation dedicated section - //if it's Dplus/Dminus, search for LS pair companions and evaluate correlations. - registry.fill(HIST("hcountDplustriggersMCGen"), 0, particle1.pt()); //to count trigger Dplus (normalisation) - for (auto& particle2 : particlesMC) { - if (std::abs(particle2.pdgCode()) != pdg::Code::kDPlus) { //check that associated is a Dplus/Dminus (both are fine) - continue; - } - if (cutYCandMax >= 0. && std::abs(RecoDecay::Y(array{particle2.px(), particle2.py(), particle2.pz()}, RecoDecay::getMassPDG(particle2.pdgCode()))) > cutYCandMax) { - continue; - } - if (cutPtCandMin >= 0. && particle2.pt() < cutPtCandMin) { - continue; - } - if (particle2.pdgCode() == particle1.pdgCode()) { //like-sign condition (both 411 or both -411) and pT_Trig>pT_assoc - //Excluding self-correlations (in principle not possible due to the '<' condition, but could rounding break it?) - if (particle1.mRowIndex == particle2.mRowIndex) { - continue; - } - entryDplusDminusPair(getDeltaPhi(particle2.phi(), particle1.phi()), - particle2.eta() - particle1.eta(), - particle1.pt(), - particle2.pt()); - } } // end inner loop - } //end outer loop - registry.fill(HIST("hcountDplusDminusPerEvent"), counterDplusDminus); + } // end outer loop + registry.fill(HIST("hCountDplusDminusPerEvent"), counterDplusDminus); } -}; - -using MCParticlesPlus2 = soa::Join; - -/// c-cbar correlator table builder - for MC gen-level analysis -struct HfCorrelatorCCbarMcGen { - - Produces entryccbarPair; - HistogramRegistry registry{ - "registry", - {{"hMCEvtCount", "Event counter - MC gen;;entries", {HistType::kTH1F, {{1, -0.5, 0.5}}}}, - {"hPtCandMCGen", "c,cbar particles - MC gen;particle #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, - {"hEtaMCGen", "c,cbar particles - MC gen;particle #it{#eta};entries", {HistType::kTH1F, {{100, -5., 5.}}}}, - {"hYMCGen", "c,cbar candidates - MC gen;candidate #it{#y};entries", {HistType::kTH1F, {{100, -5., 5.}}}}, - {"hPhiMCGen", "c,cbar particles - MC gen;particle #it{#varphi};entries", {HistType::kTH1F, {{32, 0., 2. * o2::constants::math::PI}}}}, - {"hcountCCbarPerEvent", "c,cbar particles - MC gen;Number per event;entries", {HistType::kTH1F, {{20, 0., 20.}}}}, - {"hcountCCbarPerEventPreEtaCut", "c,cbar particles - MC gen;Number per event pre #eta cut;entries", {HistType::kTH1F, {{20, 0., 20.}}}}}}; - - Configurable cutYCandMax{"cutYCandMax", -1., "max. cand. rapidity"}; - Configurable cutPtCandMin{"cutPtCandMin", -1., "min. cand. pT"}; - Configurable> bins{"ptBinsForMass", std::vector{o2::analysis::hf_cuts_d0_topik::pTBins_v}, "pT bin limits for trigger counters"}; - - void init(o2::framework::InitContext&) - { - registry.add("hcountCtriggersMCGen", "c trigger particles - MC gen;;N of trigger c quark", {HistType::kTH2F, {{1, -0.5, 0.5}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); - } + PROCESS_SWITCH(HfCorrelatorDplusDminus, processMcGen, "Process MC Gen mode", false); - void process(aod::McCollision const& mccollision, MCParticlesPlus2 const& particlesMC) + /// c-cbar correlator table builder - for MC gen-level analysis + void processccbar(aod::McCollision const& mccollision, MCParticlesPlus2Prong const& particlesMC) { registry.fill(HIST("hMCEvtCount"), 0); - int counterccbar = 0, counterccbarPreEtasel = 0; + int counterCCbar = 0, counterCCbarBeforeEtasel = 0; - //loop over particles at MC gen level + // loop over particles at MC gen level for (auto& particle1 : particlesMC) { - if (std::abs(particle1.pdgCode()) != PDG_t::kCharm) { //search c or cbar particles + if (std::abs(particle1.pdgCode()) != PDG_t::kCharm) { // search c or cbar particles continue; } - int partMothPDG = particle1.mother0_as().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 + int partMothPDG = particle1.mother0_as().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) + counterCCbarBeforeEtasel++; // count c or cbar (before kinematic selection) double yC = RecoDecay::Y(array{particle1.px(), particle1.py(), particle1.pz()}, RecoDecay::getMassPDG(particle1.pdgCode())); if (cutYCandMax >= 0. && std::abs(yC) > cutYCandMax) { continue; @@ -816,14 +471,14 @@ struct HfCorrelatorCCbarMcGen { registry.fill(HIST("hEtaMCGen"), particle1.eta()); registry.fill(HIST("hPhiMCGen"), particle1.phi()); registry.fill(HIST("hYMCGen"), yC); - counterccbar++; //count if c or cbar don't come from themselves during fragmentation (after kinematic selection) + 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. + // c-cbar correlation dedicated section + // if it's c, search for cbar and evaluate correlations. if (particle1.pdgCode() != PDG_t::kCharm) { continue; } - registry.fill(HIST("hcountCtriggersMCGen"), 0, particle1.pt()); //to count trigger c quark (for normalisation) + registry.fill(HIST("hCountCtriggersMCGen"), 0, particle1.pt()); // to count trigger c quark (for normalisation) for (auto& particle2 : particlesMC) { if (particle2.pdgCode() != PDG_t::kCharmBar) { @@ -835,136 +490,27 @@ struct HfCorrelatorCCbarMcGen { if (cutPtCandMin >= 0. && particle2.pt() < cutPtCandMin) { continue; } - //check whether mothers of quark cbar (from associated loop) are still '-4' particles - in that case the cbar quark comes from its own fragmentation, skip it - if (particle2.mother0_as().pdgCode() == PDG_t::kCharmBar) { + // 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 (particle2.mother0_as().pdgCode() == PDG_t::kCharmBar) { continue; } - entryccbarPair(getDeltaPhi(particle2.phi(), particle1.phi()), - particle2.eta() - particle1.eta(), - particle1.pt(), - particle2.pt()); - } // end inner loop - } //end outer loop - registry.fill(HIST("hcountCCbarPerEvent"), counterccbar); - registry.fill(HIST("hcountCCbarPerEventPreEtaCut"), counterccbarPreEtasel); - } -}; - -/// c-cbar correlator table builder - for MC gen-level analysis - Like Sign -struct HfCorrelatorCCbarMcGenLs { - - Produces entryccbarPair; - - HistogramRegistry registry{ - "registry", - {{"hMCEvtCount", "Event counter - MC gen;;entries", {HistType::kTH1F, {{1, -0.5, 0.5}}}}, - {"hPtCandMCGen", "c,cbar particles - MC gen;particle #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{180, 0., 36.}}}}, - {"hEtaMCGen", "c,cbar particles - MC gen;particle #it{#eta};entries", {HistType::kTH1F, {{100, -5., 5.}}}}, - {"hYMCGen", "c,cbar candidates - MC gen;candidate #it{#y};entries", {HistType::kTH1F, {{100, -5., 5.}}}}, - {"hPhiMCGen", "c,cbar particles - MC gen;particle #it{#varphi};entries", {HistType::kTH1F, {{32, 0., 2. * o2::constants::math::PI}}}}, - {"hcountCCbarPerEvent", "c,cbar particles - MC gen;Number per event;entries", {HistType::kTH1F, {{20, 0., 20.}}}}, - {"hcountCCbarPerEventPreEtaCut", "c,cbar particles - MC gen;Number per event pre #eta cut;entries", {HistType::kTH1F, {{20, 0., 20.}}}}}}; - - Configurable cutYCandMax{"cutYCandMax", -1., "max. cand. rapidity"}; - Configurable cutPtCandMin{"cutPtCandMin", -1., "min. cand. pT"}; - Configurable> bins{"ptBinsForMass", std::vector{o2::analysis::hf_cuts_d0_topik::pTBins_v}, "pT bin limits for trigger counters"}; - - void init(o2::framework::InitContext&) - { - registry.add("hcountCtriggersMCGen", "c trigger particles - MC gen;;N of trigger c quark", {HistType::kTH2F, {{1, -0.5, 0.5}, {(std::vector)bins, "#it{p}_{T} (GeV/#it{c})"}}}); + entryDplusDminusPair(getDeltaPhi(particle2.phi(), particle1.phi()), + particle2.eta() - particle1.eta(), + particle1.pt(), + particle2.pt()); + entryDplusDminusRecoInfo(1.869, + 1.869, + 8); // Dummy + } // end inner loop + } // end outer loop + registry.fill(HIST("hCountCCbarPerEvent"), counterCCbar); + registry.fill(HIST("hCountCCbarPerEventBeforeEtaCut"), counterCCbarBeforeEtasel); } - void process(aod::McCollision const& mccollision, MCParticlesPlus2 const& particlesMC) - { - registry.fill(HIST("hMCEvtCount"), 0); - int counterccbar = 0, counterccbarPreEtasel = 0; - - //loop over particles at MC gen level - for (auto& particle1 : particlesMC) { - if (std::abs(particle1.pdgCode()) != PDG_t::kCharm) { //search c or cbar particles - continue; - } - int partMothPDG = particle1.mother0_as().pdgCode(); - //check whether mothers of quark c/cbar are still '4'/'-4' particles - in that case the c/cbar quark comes from its own fragmentation, skip it - if (partMothPDG == particle1.pdgCode()) { - continue; - } - counterccbarPreEtasel++; //count c or cbar (before kinematic selection) - double yC = RecoDecay::Y(array{particle1.px(), particle1.py(), particle1.pz()}, RecoDecay::getMassPDG(particle1.pdgCode())); - if (cutYCandMax >= 0. && std::abs(yC) > cutYCandMax) { - continue; - } - if (cutPtCandMin >= 0. && particle1.pt() < cutPtCandMin) { - continue; - } - registry.fill(HIST("hPtCandMCGen"), particle1.pt()); - registry.fill(HIST("hEtaMCGen"), particle1.eta()); - registry.fill(HIST("hPhiMCGen"), particle1.phi()); - registry.fill(HIST("hYMCGen"), yC); - counterccbar++; //count if c or cbar don't come from themselves during fragmentation (after kinematic selection) - - //c-cbar correlation dedicated section - registry.fill(HIST("hcountCtriggersMCGen"), 0, particle1.pt()); //to count trigger c quark (for normalisation) - - for (auto& particle2 : particlesMC) { - if (std::abs(particle2.pdgCode()) != PDG_t::kCharm) { //search c or cbar for associated particles - continue; - } - if (cutYCandMax >= 0. && std::abs(RecoDecay::Y(array{particle2.px(), particle2.py(), particle2.pz()}, RecoDecay::getMassPDG(particle2.pdgCode()))) > cutYCandMax) { - continue; - } - if (cutPtCandMin >= 0. && particle2.pt() < cutPtCandMin) { - continue; - } - if (particle2.pdgCode() == particle1.pdgCode()) { - //check whether mothers of quark cbar (from associated loop) are still '-4' particles - in that case the cbar quark comes from its own fragmentation, skip it - if (particle2.mother0_as().pdgCode() == particle2.pdgCode()) { - continue; - } - //Excluding self-correlations (in principle not possible due to the '<' condition, but could rounding break it?) - if (particle1.mRowIndex == particle2.mRowIndex) { - continue; - } - entryccbarPair(getDeltaPhi(particle2.phi(), particle1.phi()), - particle2.eta() - particle1.eta(), - particle1.pt(), - particle2.pt()); - } // end outer if (check PDG associate) - } // end inner loop - } //end outer loop - registry.fill(HIST("hcountCCbarPerEvent"), counterccbar); - registry.fill(HIST("hcountCCbarPerEventPreEtaCut"), counterccbarPreEtasel); - } + PROCESS_SWITCH(HfCorrelatorDplusDminus, processccbar, "Process ccbar pairs", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { - WorkflowSpec workflow{}; - const bool doMCccbar = cfgc.options().get("doMCccbar"); - const bool doMCGen = cfgc.options().get("doMCGen"); - const bool doMCRec = cfgc.options().get("doMCRec"); - const bool doLikeSign = cfgc.options().get("doLikeSign"); - if (!doLikeSign) { //unlike-sign analyses - if (doMCGen) { //MC-Gen analysis - workflow.push_back(adaptAnalysisTask(cfgc)); - } else if (doMCRec) { //MC-Reco analysis - workflow.push_back(adaptAnalysisTask(cfgc)); - } else if (doMCccbar) { //MC-Reco analysis - workflow.push_back(adaptAnalysisTask(cfgc)); - } else { //data analysis - workflow.push_back(adaptAnalysisTask(cfgc)); - } - } else { //like-sign analyses - if (doMCGen) { //MC-Gen analysis - workflow.push_back(adaptAnalysisTask(cfgc)); - } else if (doMCRec) { //MC-Reco analysis - workflow.push_back(adaptAnalysisTask(cfgc)); - } else if (doMCccbar) { //MC-Reco analysis - workflow.push_back(adaptAnalysisTask(cfgc)); - } else { //data analysis - workflow.push_back(adaptAnalysisTask(cfgc)); - } - } - - return workflow; + return WorkflowSpec{adaptAnalysisTask(cfgc)}; } diff --git a/PWGHF/Tasks/taskCorrelationDDbar.cxx b/PWGHF/Tasks/taskCorrelationDDbar.cxx index 80827078a37..6208dd18b1e 100644 --- a/PWGHF/Tasks/taskCorrelationDDbar.cxx +++ b/PWGHF/Tasks/taskCorrelationDDbar.cxx @@ -32,14 +32,6 @@ namespace o2::aod using DDbarPairFull = soa::Join; } // namespace o2::aod -void customize(std::vector& workflowOptions) -{ - ConfigParamSpec optionDoMCGen{"doMCGen", VariantType::Bool, false, {"Run MC-Gen dedicated tasks."}}; - ConfigParamSpec optionDoMCRec{"doMCRec", VariantType::Bool, true, {"Run MC-Rec dedicated tasks."}}; - workflowOptions.push_back(optionDoMCGen); - workflowOptions.push_back(optionDoMCRec); -} - #include "Framework/runDataProcessing.h" /// @@ -91,8 +83,6 @@ const int npTBinsEfficiency = o2::analysis::hf_cuts_d0_topik::npTBins; const double efficiencyDmesonDefault[npTBinsEfficiency] = {}; auto efficiencyDmeson_v = std::vector{efficiencyDmesonDefault, efficiencyDmesonDefault + npTBinsEfficiency}; -/// D-Dbar correlation pair filling task, from pair tables - for real data and data-like analysis (i.e. reco-level w/o matching request via MC truth) -/// Works on both USL and LS analyses pair tables struct HfTaskCorrelationDDbar { HistogramRegistry registry{ @@ -110,100 +100,8 @@ struct HfTaskCorrelationDDbar { {"hCorrel2DPtIntSidebands", stringDDbar + stringSideband + stringDeltaPhi + stringDeltaEta + "entries", {HistType::kTH2F, {{64, -o2::constants::math::PI / 2., 3. * o2::constants::math::PI / 2.}, {200, -10., 10.}}}}, {"hCorrel2DVsPtSidebands", stringDDbar + stringSideband + stringDeltaPhi + stringDeltaEta + stringPtD + stringPtDbar + "entries", {HistType::kTHnSparseD, {{64, -o2::constants::math::PI / 2., 3. * o2::constants::math::PI / 2.}, {120, -6., 6.}, {10, 0., 10.}, {10, 0., 10.}}}}, //note: axes 3 and 4 (the pT) are updated in the init() {"hDeltaPtDDbarSidebands", stringDDbar + stringSideband + stringDeltaPt + "entries", {HistType::kTH1F, {{144, -36., 36.}}}}, - {"hDeltaPtMaxMinSidebands", stringDDbar + stringSideband + stringDeltaPtMaxMin + "entries", {HistType::kTH1F, {{72, 0., 36.}}}}}}; - - //pT ranges for correlation plots: the default values are those embedded in hf_cuts_d0_topik (i.e. the mass pT bins), but can be redefined via json files - Configurable> binsCorrelations{"ptBinsForCorrelations", std::vector{pTBinsCorrelations_v}, "pT bin limits for correlation plots"}; - //pT bins for effiencies: same as above - Configurable> binsEfficiency{"ptBinsForEfficiency", std::vector{o2::analysis::hf_cuts_d0_topik::pTBins_v}, "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> efficiencyDmeson{"efficiencyDmeson", std::vector{efficiencyDmeson_v}, "Efficiency values for D meson specie under study"}; - Configurable flagApplyEfficiency{"efficiencyFlagD0", 1, "Flag for applying efficiency weights"}; - - 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); - registry.get(HIST("hMass2DCorrelationPairs"))->Sumw2(); - registry.get(HIST("hCorrel2DVsPtSignalRegion"))->Sumw2(); - registry.get(HIST("hCorrel2DVsPtSidebands"))->Sumw2(); - } - } - - void process(aod::DDbarPairFull const& pairEntries) - { - for (auto& pairEntry : pairEntries) { - //define variables for widely used quantities - double deltaPhi = pairEntry.deltaPhi(); - double deltaEta = pairEntry.deltaEta(); - double ptD = pairEntry.ptD(); - double ptDbar = pairEntry.ptDbar(); - double massD = pairEntry.mD(); - double massDbar = pairEntry.mDbar(); - - int pTBinD = o2::analysis::findBin(binsCorrelations, ptD); - int pTBinDbar = o2::analysis::findBin(binsCorrelations, ptDbar); - - double efficiencyWeight = 1.; - if (flagApplyEfficiency) { - efficiencyWeight = 1. / (efficiencyDmeson->at(o2::analysis::findBin(binsEfficiency, ptD)) * efficiencyDmeson->at(o2::analysis::findBin(binsEfficiency, 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 - } -}; - -/// D-Dbar correlation pair filling task, from pair tables - for MC reco-level analysis (candidates matched to true signal only, but also bkg sources are studied) -/// Works on both USL and LS analyses pair tables -struct HfTaskCorrelationDDbarMcRec { - - HistogramRegistry registry{ - "registry", - //NOTE: use hMassD0 (from correlator task) for normalisation, and hMass2DCorrelationPairs for 2D-sideband-subtraction purposes - {{"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() + {"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() @@ -231,15 +129,21 @@ struct HfTaskCorrelationDDbarMcRec { {"hCorrel2DPtIntSidebandsMCRec", stringMCReco + stringSideband + stringDeltaPhi + stringDeltaEta + "entries", {HistType::kTH2F, {{64, -o2::constants::math::PI / 2., 3. * o2::constants::math::PI / 2.}, {200, -10., 10.}}}}, {"hDeltaPtDDbarSidebandsMCRec", stringMCReco + stringSideband + stringDeltaPt + "entries", {HistType::kTH1F, {{144, -36., 36.}}}}, {"hDeltaPtMaxMinSidebandsMCRec", stringMCReco + stringSideband + stringDeltaPtMaxMin + "entries", {HistType::kTH1F, {{72, 0., 36.}}}}, - {"hCorrel2DVsPtSidebandsMCRecBkgBkg", stringMCReco + "BkgBkg" + stringSideband + stringDeltaPhi + stringDeltaEta + stringPtD + stringPtDbar + "entries", {HistType::kTHnSparseD, {{64, -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() - {"hCorrel2DVsPtSidebandsMCRecBkgRef", stringMCReco + "BkgRef" + stringSideband + stringDeltaPhi + stringDeltaEta + stringPtD + stringPtDbar + "entries", {HistType::kTHnSparseD, {{64, -o2::constants::math::PI / 2., 3. * o2::constants::math::PI / 2.}, {120, -6., 6.}, {10, 0., 10.}, {10, 0., 10.}}}}, //note: axes 3 and 4 (the pT) are updated in the init() - {"hCorrel2DVsPtSidebandsMCRecBkgSig", stringMCReco + "BkgSig" + stringSideband + stringDeltaPhi + stringDeltaEta + stringPtD + stringPtDbar + "entries", {HistType::kTHnSparseD, {{64, -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() - {"hCorrel2DVsPtSidebandsMCRecRefBkg", stringMCReco + "RefBkg" + stringSideband + stringDeltaPhi + stringDeltaEta + stringPtD + stringPtDbar + "entries", {HistType::kTHnSparseD, {{64, -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() - {"hCorrel2DVsPtSidebandsMCRecRefRef", stringMCReco + "RefRef" + stringSideband + stringDeltaPhi + stringDeltaEta + stringPtD + stringPtDbar + "entries", {HistType::kTHnSparseD, {{64, -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() - {"hCorrel2DVsPtSidebandsMCRecRefSig", stringMCReco + "RefSig" + stringSideband + stringDeltaPhi + stringDeltaEta + stringPtD + stringPtDbar + "entries", {HistType::kTHnSparseD, {{64, -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() - {"hCorrel2DVsPtSidebandsMCRecSigBkg", stringMCReco + "SigBkg" + stringSideband + stringDeltaPhi + stringDeltaEta + stringPtD + stringPtDbar + "entries", {HistType::kTHnSparseD, {{64, -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() - {"hCorrel2DVsPtSidebandsMCRecSigRef", stringMCReco + "SigRef" + stringSideband + stringDeltaPhi + stringDeltaEta + stringPtD + stringPtDbar + "entries", {HistType::kTHnSparseD, {{64, -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() - {"hCorrel2DVsPtSidebandsMCRecSigSig", stringMCReco + "SigSig" + stringSideband + stringDeltaPhi + stringDeltaEta + stringPtD + stringPtDbar + "entries", {HistType::kTHnSparseD, {{64, -o2::constants::math::PI / 2., 3. * o2::constants::math::PI / 2.}, {120, -6., 6.}, {10, 0., 10.}, {10, 0., 10.}}}}}}; //note: axes 3 and 4 (the pT) are updated in the init() + {"hCorrel2DVsPtSidebandsMCRecBkgBkg", stringMCReco + "BkgBkg" + stringSideband + stringDeltaPhi + stringDeltaEta + stringPtD + stringPtDbar + "entries", {HistType::kTHnSparseD, {{64, -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() + {"hCorrel2DVsPtSidebandsMCRecBkgRef", stringMCReco + "BkgRef" + stringSideband + stringDeltaPhi + stringDeltaEta + stringPtD + stringPtDbar + "entries", {HistType::kTHnSparseD, {{64, -o2::constants::math::PI / 2., 3. * o2::constants::math::PI / 2.}, {120, -6., 6.}, {10, 0., 10.}, {10, 0., 10.}}}}, //note: axes 3 and 4 (the pT) are updated in the init() + {"hCorrel2DVsPtSidebandsMCRecBkgSig", stringMCReco + "BkgSig" + stringSideband + stringDeltaPhi + stringDeltaEta + stringPtD + stringPtDbar + "entries", {HistType::kTHnSparseD, {{64, -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() + {"hCorrel2DVsPtSidebandsMCRecRefBkg", stringMCReco + "RefBkg" + stringSideband + stringDeltaPhi + stringDeltaEta + stringPtD + stringPtDbar + "entries", {HistType::kTHnSparseD, {{64, -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() + {"hCorrel2DVsPtSidebandsMCRecRefRef", stringMCReco + "RefRef" + stringSideband + stringDeltaPhi + stringDeltaEta + stringPtD + stringPtDbar + "entries", {HistType::kTHnSparseD, {{64, -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() + {"hCorrel2DVsPtSidebandsMCRecRefSig", stringMCReco + "RefSig" + stringSideband + stringDeltaPhi + stringDeltaEta + stringPtD + stringPtDbar + "entries", {HistType::kTHnSparseD, {{64, -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() + {"hCorrel2DVsPtSidebandsMCRecSigBkg", stringMCReco + "SigBkg" + stringSideband + stringDeltaPhi + stringDeltaEta + stringPtD + stringPtDbar + "entries", {HistType::kTHnSparseD, {{64, -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() + {"hCorrel2DVsPtSidebandsMCRecSigRef", stringMCReco + "SigRef" + stringSideband + stringDeltaPhi + stringDeltaEta + stringPtD + stringPtDbar + "entries", {HistType::kTHnSparseD, {{64, -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() + {"hCorrel2DVsPtSidebandsMCRecSigSig", stringMCReco + "SigSig" + stringSideband + stringDeltaPhi + stringDeltaEta + stringPtD + stringPtDbar + "entries", {HistType::kTHnSparseD, {{64, -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() + {"hDeltaEtaPtIntMCGen", stringMCParticles + stringDeltaEta + "entries", {HistType::kTH1F, {{200, -10., 10.}}}}, + {"hDeltaPhiPtIntMCGen", stringMCParticles + stringDeltaPhi + "entries", {HistType::kTH1F, {{64, -o2::constants::math::PI / 2., 3. * o2::constants::math::PI / 2.}}}}, + {"hCorrel2DPtIntMCGen", stringMCParticles + stringDeltaPhi + stringDeltaEta + "entries", {HistType::kTH2F, {{64, -o2::constants::math::PI / 2., 3. * o2::constants::math::PI / 2.}, {200, -10., 10.}}}}, + {"hCorrel2DVsPtMCGen", stringMCParticles + stringDeltaPhi + stringDeltaEta + stringPtD + stringPtDbar + "entries", {HistType::kTHnSparseD, {{64, -o2::constants::math::PI / 2., 3. * o2::constants::math::PI / 2.}, {120, -6., 6.}, {10, 0., 10.}, {10, 0., 10.}}}}, //note: axes 3 and 4 (the pT) are updated in the init() + {"hDeltaPtDDbarMCGen", stringMCParticles + stringDeltaPt + "entries", {HistType::kTH1F, {{144, -36., 36.}}}}, + {"hDeltaPtMaxMinMCGen", stringMCParticles + stringDeltaPtMaxMin + "entries", {HistType::kTH1F, {{72, 0., 36.}}}}}}; //pT ranges for correlation plots: the default values are those embedded in hf_cuts_d0_topik (i.e. the mass pT bins), but can be redefined via json files Configurable> binsCorrelations{"ptBinsForCorrelations", std::vector{pTBinsCorrelations_v}, "pT bin limits for correlation plots"}; @@ -253,7 +157,7 @@ struct HfTaskCorrelationDDbarMcRec { 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> efficiencyDmeson{"efficiencyDmeson", std::vector{efficiencyDmeson_v}, "Efficiency values for D meson specie under study"}; - Configurable flagApplyEfficiency{"efficiencyFlagD", 1, "Flag for applying D-meson efficiency weights"}; + Configurable flagApplyEfficiency{"efficiencyFlagD", 1, "Flag for applying efficiency weights"}; void init(o2::framework::InitContext&) { @@ -262,6 +166,12 @@ struct HfTaskCorrelationDDbarMcRec { 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); + 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); @@ -316,10 +226,71 @@ struct HfTaskCorrelationDDbarMcRec { registry.get(HIST("hCorrel2DVsPtSidebandsMCRecSigBkg"))->Sumw2(); registry.get(HIST("hCorrel2DVsPtSidebandsMCRecSigRef"))->Sumw2(); registry.get(HIST("hCorrel2DVsPtSidebandsMCRecSigSig"))->Sumw2(); + registry.get(HIST("hCorrel2DVsPtMCGen"))->GetAxis(i)->Set(nBinspTaxis, valuespTaxis); + registry.get(HIST("hCorrel2DVsPtMCGen"))->Sumw2(); } } - void process(aod::DDbarPairFull const& pairEntries) + /// 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 (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(binsCorrelations, ptD); + int pTBinDbar = o2::analysis::findBin(binsCorrelations, ptDbar); + + double efficiencyWeight = 1.; + if (flagApplyEfficiency) { + efficiencyWeight = 1. / (efficiencyDmeson->at(o2::analysis::findBin(binsEfficiency, ptD)) * efficiencyDmeson->at(o2::analysis::findBin(binsEfficiency, 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(HfTaskCorrelationDDbar, processData, "Process data", false); + + /// 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 (auto& pairEntry : pairEntries) { //define variables for widely used quantities @@ -461,37 +432,12 @@ struct HfTaskCorrelationDDbarMcRec { } } //end loop } -}; -/// D-Dbar correlation pair filling task, from pair tables - for MC gen-level analysis (no filter/selection, only true signal) - Ok for both USL and LS analyses -/// Works on both USL and LS analyses pair tables (and if tables are filled with quark pairs as well) -struct HfTaskCorrelationDDbarMcGen { + PROCESS_SWITCH(HfTaskCorrelationDDbar, processMcRec, "Process MC Reco mode", true); - HistogramRegistry registry{ - "registry", - {{"hDeltaEtaPtIntMCGen", stringMCParticles + stringDeltaEta + "entries", {HistType::kTH1F, {{200, -10., 10.}}}}, - {"hDeltaPhiPtIntMCGen", stringMCParticles + stringDeltaPhi + "entries", {HistType::kTH1F, {{64, -o2::constants::math::PI / 2., 3. * o2::constants::math::PI / 2.}}}}, - {"hCorrel2DPtIntMCGen", stringMCParticles + stringDeltaPhi + stringDeltaEta + "entries", {HistType::kTH2F, {{64, -o2::constants::math::PI / 2., 3. * o2::constants::math::PI / 2.}, {200, -10., 10.}}}}, - {"hCorrel2DVsPtMCGen", stringMCParticles + stringDeltaPhi + stringDeltaEta + stringPtD + stringPtDbar + "entries", {HistType::kTHnSparseD, {{64, -o2::constants::math::PI / 2., 3. * o2::constants::math::PI / 2.}, {120, -6., 6.}, {10, 0., 10.}, {10, 0., 10.}}}}, //note: axes 3 and 4 (the pT) are updated in the init() - {"hDeltaPtDDbarMCGen", stringMCParticles + stringDeltaPt + "entries", {HistType::kTH1F, {{144, -36., 36.}}}}, - {"hDeltaPtMaxMinMCGen", stringMCParticles + stringDeltaPtMaxMin + "entries", {HistType::kTH1F, {{72, 0., 36.}}}}}}; - - //pT ranges for correlation plots: the default values are those embedded in hf_cuts_d0_topik (i.e. the mass pT bins), but can be redefined via json files - Configurable> binsCorrelations{"ptBinsForCorrelations", std::vector{pTBinsCorrelations_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); - registry.get(HIST("hCorrel2DVsPtMCGen"))->Sumw2(); - } - } - - void process(aod::DDbarPair const& pairEntries) + /// D-Dbar correlation pair filling task, from pair tables - for MC gen-level analysis (no filter/selection, only true signal) - Ok for both USL and LS analyses + /// Works on both USL and LS analyses pair tables (and if tables are filled with quark pairs as well) + void processMcGen(aod::DDbarPair const& pairEntries) { for (auto& pairEntry : pairEntries) { //define variables for widely used quantities @@ -513,20 +459,11 @@ struct HfTaskCorrelationDDbarMcGen { registry.fill(HIST("hDeltaPtMaxMinMCGen"), std::abs(ptDbar - ptD)); } //end loop } + + PROCESS_SWITCH(HfTaskCorrelationDDbar, processMcGen, "Process MC Gen mode", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { - WorkflowSpec workflow{}; - //MC-based tasks - const bool doMCGen = cfgc.options().get("doMCGen"); - const bool doMCRec = cfgc.options().get("doMCRec"); - if (doMCGen) { //MC-Gen analysis - workflow.push_back(adaptAnalysisTask(cfgc)); - } else if (doMCRec) { //MC-Rec analysis - workflow.push_back(adaptAnalysisTask(cfgc)); - } else { //data analysis - workflow.push_back(adaptAnalysisTask(cfgc)); - } - return workflow; + return WorkflowSpec{adaptAnalysisTask(cfgc, TaskName{"hf-task-correlation-d-dbar"})}; }