Skip to content

Commit b573546

Browse files
author
fcolamar
committed
ALICE 3: adding support for MC reco mode in DDbar correlation via OTF sims
1 parent 28a13a8 commit b573546

File tree

4 files changed

+98
-109
lines changed

4 files changed

+98
-109
lines changed

‎ALICE3/DataModel/A3DecayFinderTables.h‎

Lines changed: 19 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -85,10 +85,7 @@ DECLARE_SOA_DYNAMIC_COLUMN(E, e, //!
8585
[](float px, float py, float pz, double m) -> float { return RecoDecay::e(px, py, pz, m); });
8686
DECLARE_SOA_COLUMN(Eta, eta, float); //!
8787
DECLARE_SOA_COLUMN(Phi, phi, float); //!
88-
DECLARE_SOA_DYNAMIC_COLUMN(Y, y, //!
89-
[](float px, float py, float pz, double m) -> float { return RecoDecay::y(std::array{px, py, pz}, m); });
90-
DECLARE_SOA_COLUMN(IsSelD0, isSelD0, int); //!
91-
DECLARE_SOA_COLUMN(IsSelD0bar, isSelD0bar, int); //!
88+
DECLARE_SOA_COLUMN(Y, y, float);
9289
} // namespace a3D0meson
9390
DECLARE_SOA_TABLE(Alice3D0Meson, "AOD", "ALICE3D0MESON", //!
9491
o2::soa::Index<>,
@@ -103,9 +100,24 @@ DECLARE_SOA_TABLE(Alice3D0Meson, "AOD", "ALICE3D0MESON", //!
103100
a3D0meson::E<a3D0meson::Px, a3D0meson::Py, a3D0meson::Pz, a3D0meson::M>,
104101
a3D0meson::Eta,
105102
a3D0meson::Phi,
106-
a3D0meson::Y<a3D0meson::Px, a3D0meson::Py, a3D0meson::Pz, a3D0meson::M>,
107-
a3D0meson::IsSelD0,
108-
a3D0meson::IsSelD0bar);
103+
a3D0meson::Y);
104+
105+
namespace a3D0Selection
106+
{
107+
DECLARE_SOA_COLUMN(IsSelD0, isSelD0, int); //!
108+
DECLARE_SOA_COLUMN(IsSelD0bar, isSelD0bar, int); //!
109+
} // namespace a3D0Selection
110+
DECLARE_SOA_TABLE(Alice3D0Sel, "AOD", "ALICE3D0SEL", //!
111+
a3D0Selection::IsSelD0,
112+
a3D0Selection::IsSelD0bar);
113+
114+
namespace a3D0MCTruth
115+
{
116+
DECLARE_SOA_COLUMN(McTruthInfo, mcTruthInfo, int); //! 0 for bkg, 1 for true D0, 2 for true D0bar
117+
} // namespace a3D0MCTruth
118+
DECLARE_SOA_TABLE(Alice3D0MCTruth, "AOD", "ALICE3D0MCTRUTH", //!
119+
a3D0MCTruth::McTruthInfo); //!
120+
109121
} // namespace o2::aod
110122

111123
#endif // ALICE3_DATAMODEL_A3DECAYFINDERTABLES_H_

‎ALICE3/TableProducer/alice3-correlatorDDbar.cxx‎

Lines changed: 42 additions & 81 deletions
Original file line numberDiff line numberDiff line change
@@ -79,15 +79,15 @@ struct alice3correlatorddbar {
7979
Configurable<int> selectionFlagD0{"selectionFlagD0", 1, "Selection Flag for D0"};
8080
Configurable<int> selectionFlagD0bar{"selectionFlagD0bar", 1, "Selection Flag for D0bar"};
8181
Configurable<int> applyEfficiency{"applyEfficiency", 1, "Flag for applying D-meson efficiency weights"};
82-
Configurable<double> yCandMax{"yCandMax", -1., "max. cand. rapidity"};
83-
Configurable<double> ptCandMin{"ptCandMin", -1., "min. cand. pT"};
82+
Configurable<float> yCandMax{"yCandMax", 999., "max. cand. rapidity"};
83+
Configurable<float> ptCandMin{"ptCandMin", -1., "min. cand. pT"};
8484
Configurable<std::vector<double>> binsPt{"binsPt", std::vector<double>{o2::analysis::hf_cuts_d0_to_pi_k::vecBinsPt}, "pT bin limits for candidate mass plots and efficiency"};
8585
Configurable<std::vector<double>> efficiencyD{"efficiencyD", std::vector<double>{efficiencyDmeson_v}, "Efficiency values for D0 meson"};
8686

8787
//HfHelper hfHelper; //not needed for now
8888

89-
Partition<aod::Alice3D0Meson> selectedCandidates = aod::a3D0meson::isSelD0 >= selectionFlagD0 || aod::a3D0meson::isSelD0bar >= selectionFlagD0bar;
90-
//Partition<soa::Join<aod::Alice3D0Meson, aod::HfCand2ProngMcRec>> selectedCandidatesMC = aod::a3D0meson::isSelD0 >= selectionFlagD0 || aod::a3D0meson::isSelD0bar >= selectionFlagD0bar; //MCRec case (and related columns) to be implemented
89+
Partition<soa::Join<aod::Alice3D0Meson, aod::Alice3D0Sel>> selectedCandidates = aod::a3D0meson::y > -yCandMax && aod::a3D0meson::y < yCandMax && aod::a3D0meson::pt > ptCandMin && (aod::a3D0Selection::isSelD0 >= selectionFlagD0 || aod::a3D0Selection::isSelD0bar >= selectionFlagD0bar);
90+
Partition<soa::Join<aod::Alice3D0Meson, aod::Alice3D0Sel, aod::Alice3D0MCTruth>> selectedCandidatesMC = aod::a3D0meson::y > -yCandMax && aod::a3D0meson::y < yCandMax && aod::a3D0meson::pt > ptCandMin && (aod::a3D0Selection::isSelD0 >= selectionFlagD0 || aod::a3D0Selection::isSelD0bar >= selectionFlagD0bar);
9191

9292
HistogramRegistry registry{
9393
"registry",
@@ -106,14 +106,7 @@ struct alice3correlatorddbar {
106106
{"hSelectionStatusMCRec", "D0,D0bar candidates - MC reco;selection status;entries", {HistType::kTH1F, {{4, -0.5, 3.5}}}},
107107
{"hEtaMCRec", "D0,D0bar candidates - MC reco;candidate #it{#eta};entries", {HistType::kTH1F, {{yAxisBins, yAxisMin, yAxisMax}}}},
108108
{"hPhiMCRec", "D0,D0bar candidates - MC reco;candidate #it{#varphi};entries", {HistType::kTH1F, {{phiAxisBins, phiAxisMin, phiAxisMax}}}},
109-
{"hYMCRec", "D0,D0bar candidates - MC reco;candidate #it{y};entries", {HistType::kTH1F, {{yAxisBins, yAxisMin, yAxisMax}}}},
110-
{"hMCEvtCount", "Event counter - MC gen;;entries", {HistType::kTH1F, {{1, -0.5, 0.5}}}},
111-
{"hPtCandMCGen", "D0,D0bar particles - MC gen;particle #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{ptDAxisBins, ptDAxisMin, ptDAxisMax}}}},
112-
{"hEtaMCGen", "D0,D0bar particles - MC gen;particle #it{#eta};entries", {HistType::kTH1F, {{yAxisBins, yAxisMin, yAxisMax}}}},
113-
{"hPhiMCGen", "D0,D0bar particles - MC gen;particle #it{#varphi};entries", {HistType::kTH1F, {{phiAxisBins, phiAxisMin, phiAxisMax}}}},
114-
{"hYMCGen", "D0,D0bar candidates - MC gen;candidate #it{y};entries", {HistType::kTH1F, {{yAxisBins, yAxisMin, yAxisMax}}}},
115-
{"hCountD0D0barPerEvent", "D0,D0bar particles - MC gen;Number per event;entries", {HistType::kTH1F, {{20, 0., 20.}}}},
116-
{"hDDbarVsDaughterEtaCut", "D0,D0bar pairs vs #eta cut on D daughters;#eta_{max};candidates #it{p}_{T} threshold (GeV/#it{c});entries", {HistType::kTH2F, {{static_cast<int>(maxEtaCut / incrementEtaCut), 0., maxEtaCut}, {static_cast<int>(ptThresholdForMaxEtaCut / incrementPtThreshold), 0., ptThresholdForMaxEtaCut}}}}}};
109+
{"hYMCRec", "D0,D0bar candidates - MC reco;candidate #it{y};entries", {HistType::kTH1F, {{yAxisBins, yAxisMin, yAxisMax}}}}}};
117110

118111
void init(InitContext&)
119112
{
@@ -127,24 +120,15 @@ struct alice3correlatorddbar {
127120
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})"}}});
128121
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})"}}});
129122
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})"}}});
130-
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})"}}});
131-
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})"}}});
132123
}
133124

134125
/// D0-D0bar correlation pair builder - for real data and data-like analysis (i.e. reco-level w/o matching request via MC truth)
135126
void processData(aod::Collision const& collision,
136-
aod::Alice3D0Meson const&)
127+
soa::Join<aod::Alice3D0Meson, aod::Alice3D0Sel> const&)
137128
{
138129
auto selectedCandidatesGrouped = selectedCandidates->sliceByCached(aod::a3D0meson::collisionId, collision.globalIndex(), cache);
139130

140131
for (const auto& candidate1 : selectedCandidatesGrouped) { //loop over reconstructed and selected D0 and D0bar (together, to fill mass plots first)
141-
if (yCandMax >= 0. && std::abs(candidate1.y()) > yCandMax) { //no ambiguity on Y since each candidate is either D0 or D0bar, never both
142-
continue;
143-
}
144-
if (ptCandMin >= 0. && candidate1.pt() < ptCandMin) {
145-
continue;
146-
}
147-
148132
double efficiencyWeight = 1.;
149133
if (applyEfficiency) {
150134
efficiencyWeight = 1. / efficiencyD->at(o2::analysis::findBin(binsPt, candidate1.pt()));
@@ -165,7 +149,7 @@ struct alice3correlatorddbar {
165149
registry.fill(HIST("hEta"), candidate1.eta());
166150
registry.fill(HIST("hPhi"), candidate1.phi());
167151
registry.fill(HIST("hY"), candidate1.y());
168-
registry.fill(HIST("hSelectionStatus"), candidate1.isSelD0bar() + (candidate1.isSelD0() * 2));
152+
registry.fill(HIST("hSelectionStatus"), candidate1.isSelD0() + (candidate1.isSelD0bar() * 2));
169153

170154
// D-Dbar correlation dedicated section
171155

@@ -177,19 +161,11 @@ struct alice3correlatorddbar {
177161
if (candidate2.isSelD0bar() < selectionFlagD0bar) { // keep only D0bar candidates passing the selection
178162
continue;
179163
}
180-
// kinematic selection on D0bar candidates
181-
if (yCandMax >= 0. && std::abs(candidate2.y()) > yCandMax) {
182-
continue;
183-
}
184-
if (ptCandMin >= 0. && candidate2.pt() < ptCandMin) {
185-
continue;
186-
}
187164
// excluding trigger self-correlations (possible in case of both mass hypotheses accepted)
188165
if (candidate1.mRowIndex == candidate2.mRowIndex) { //this by definition should never happen, since each candidate is either D0 or D0bar
189166
continue;
190167
}
191-
if (candidate1.mRowIndex == candidate2.mRowIndex) { //revised, temporary condition to avoid self-correlations (the best would be check the daughterIDs, but we don't store them at the moment)
192-
if((candidate1.pt() - candidate2.pt()) < 1e-5 && (candidate1.eta() - candidate2.eta()) < 1e-5 && (candidate1.phi() - candidate2.phi()) < 1e-5)
168+
if((candidate1.pt() - candidate2.pt()) < 1e-5 && (candidate1.eta() - candidate2.eta()) < 1e-5 && (candidate1.phi() - candidate2.phi()) < 1e-5) { //revised, temporary condition to avoid self-correlations (the best would be check the daughterIDs, but we don't store them at the moment)
193169
continue;
194170
}
195171
entryD0D0barPair(getDeltaPhi(candidate2.phi(), candidate1.phi()),
@@ -220,63 +196,51 @@ struct alice3correlatorddbar {
220196
}
221197

222198
PROCESS_SWITCH(alice3correlatorddbar, processData, "Process data", true);
223-
/*
199+
224200
/// D0-D0bar correlation pair builder - for MC reco-level analysis (candidates matched to true signal only, but also the various bkg sources are studied)
225201
void processMcRec(aod::Collision const& collision,
226-
aod::TracksWDca const& tracks,
227-
soa::Join<aod::HfCand2Prong, aod::HfSelD0, aod::HfCand2ProngMcRec> const&)
202+
soa::Join<aod::Alice3D0Meson, aod::Alice3D0Sel, aod::Alice3D0MCTruth> const&)
228203
{
229-
auto selectedD0CandidatesGroupedMC = selectedD0candidatesMC->sliceByCached(aod::hf_cand::collisionId, collision.globalIndex(), cache);
204+
auto selectedCandidatesGroupedMC = selectedCandidatesMC->sliceByCached(aod::a3D0meson::collisionId, collision.globalIndex(), cache);
230205

231206
// MC reco level
232207
bool flagD0Signal = false;
233208
bool flagD0Reflection = false;
234209
bool flagD0barSignal = false;
235210
bool flagD0barReflection = false;
236-
for (const auto& candidate1 : selectedD0CandidatesGroupedMC) {
237-
// check decay channel flag for candidate1
238-
if (!(candidate1.hfflag() & 1 << aod::hf_cand_2prong::DecayType::D0ToPiK)) {
239-
continue;
240-
}
241-
if (yCandMax >= 0. && std::abs(hfHelper.yD0(candidate1)) > yCandMax) {
242-
continue;
243-
}
244-
if (ptCandMin >= 0. && candidate1.pt() < ptCandMin) {
245-
continue;
246-
}
247-
211+
for (const auto& candidate1 : selectedCandidatesGroupedMC) {
248212
double efficiencyWeight = 1.;
249213
if (applyEfficiency) {
250214
efficiencyWeight = 1. / efficiencyD->at(o2::analysis::findBin(binsPt, candidate1.pt()));
251215
}
252216

253-
if (std::abs(candidate1.flagMcMatchRec()) == 1 << aod::hf_cand_2prong::DecayType::D0ToPiK) {
217+
if (candidate1.mcTruthInfo()) { //1 or 2, i.e. true D0 or D0bar
254218
// fill per-candidate distributions from D0/D0bar true candidates
255219
registry.fill(HIST("hPtCandMCRec"), candidate1.pt());
256220
registry.fill(HIST("hPtProng0MCRec"), candidate1.ptProng0());
257221
registry.fill(HIST("hPtProng1MCRec"), candidate1.ptProng1());
258222
registry.fill(HIST("hEtaMCRec"), candidate1.eta());
259223
registry.fill(HIST("hPhiMCRec"), candidate1.phi());
260-
registry.fill(HIST("hYMCRec"), hfHelper.yD0(candidate1));
261-
registry.fill(HIST("hSelectionStatusMCRec"), candidate1.isSelD0bar() + (candidate1.isSelD0() * 2));
224+
registry.fill(HIST("hYMCRec"), candidate1.y());
225+
registry.fill(HIST("hSelectionStatusMCRec"), candidate1.isSelD0() + (candidate1.isSelD0bar() * 2));
262226
}
263227
// fill invariant mass plots from D0/D0bar signal and background candidates
264-
if (candidate1.isSelD0() >= selectionFlagD0) { // only reco as D0
265-
if (candidate1.flagMcMatchRec() == 1 << aod::hf_cand_2prong::DecayType::D0ToPiK) { // also matched as D0
266-
registry.fill(HIST("hMassD0MCRecSig"), hfHelper.invMassD0ToPiK(candidate1), candidate1.pt(), efficiencyWeight);
267-
} else if (candidate1.flagMcMatchRec() == -(1 << aod::hf_cand_2prong::DecayType::D0ToPiK)) {
268-
registry.fill(HIST("hMassD0MCRecRefl"), hfHelper.invMassD0ToPiK(candidate1), candidate1.pt(), efficiencyWeight);
228+
if (candidate1.isSelD0() >= selectionFlagD0) { // only reco as D0
229+
if (candidate1.mcTruthInfo() == 1) { // also matched as D0
230+
registry.fill(HIST("hMassD0MCRecSig"), candidate1.m(), candidate1.pt(), efficiencyWeight); //here m is univoque, since a given candidate passes the selection with only a single mass option
231+
} else if (candidate1.mcTruthInfo() == 2) {
232+
registry.fill(HIST("hMassD0MCRecRefl"), candidate1.m(), candidate1.pt(), efficiencyWeight);
269233
} else {
270-
registry.fill(HIST("hMassD0MCRecBkg"), hfHelper.invMassD0ToPiK(candidate1), candidate1.pt(), efficiencyWeight);
234+
registry.fill(HIST("hMassD0MCRecBkg"), candidate1.m(), candidate1.pt(), efficiencyWeight);
271235
}
272236
}
273-
if (candidate1.isSelD0bar() >= selectionFlagD0bar) { // only reco as D0bar
274-
if (candidate1.flagMcMatchRec() == -(1 << aod::hf_cand_2prong::DecayType::D0ToPiK)) { // also matched as D0bar
275-
registry.fill(HIST("hMassD0barMCRecSig"), hfHelper.invMassD0barToKPi(candidate1), candidate1.pt(), efficiencyWeight);
276-
} else if (candidate1.flagMcMatchRec() == 1 << aod::hf_cand_2prong::DecayType::D0ToPiK) {
277-
registry.fill(HIST("hMassD0barMCRecRefl"), hfHelper.invMassD0barToKPi(candidate1), candidate1.pt(), efficiencyWeight);
237+
if (candidate1.isSelD0bar() >= selectionFlagD0bar) { // only reco as D0bar
238+
if (candidate1.mcTruthInfo() == 2) { // also matched as D0bar
239+
registry.fill(HIST("hMassD0barMCRecSig"), candidate1.m(), candidate1.pt(), efficiencyWeight); //here m is univoque, since a given candidate passes the selection with only a single mass option
240+
} else if (candidate1.mcTruthInfo() == 1) {
241+
registry.fill(HIST("hMassD0barMCRecRefl"), candidate1.m(), candidate1.pt(), efficiencyWeight);
278242
} else {
279-
registry.fill(HIST("hMassD0barMCRecBkg"), hfHelper.invMassD0barToKPi(candidate1), candidate1.pt(), efficiencyWeight);
243+
registry.fill(HIST("hMassD0barMCRecBkg"), candidate1.m(), candidate1.pt(), efficiencyWeight);
280244
}
281245
}
282246

@@ -285,27 +249,24 @@ struct alice3correlatorddbar {
285249
if (candidate1.isSelD0() < selectionFlagD0) { // discard candidates not selected as D0 in outer loop
286250
continue;
287251
}
288-
flagD0Signal = candidate1.flagMcMatchRec() == 1 << aod::hf_cand_2prong::DecayType::D0ToPiK; // flagD0Signal 'true' if candidate1 matched to D0 (particle)
289-
flagD0Reflection = candidate1.flagMcMatchRec() == -(1 << aod::hf_cand_2prong::DecayType::D0ToPiK); // flagD0Reflection 'true' if candidate1, selected as D0 (particle), is matched to D0bar (antiparticle)
290-
for (const auto& candidate2 : selectedD0CandidatesGroupedMC) {
291-
if (!(candidate2.hfflag() & 1 << aod::hf_cand_2prong::DecayType::D0ToPiK)) { // check decay channel flag for candidate2
292-
continue;
293-
}
252+
253+
flagD0Signal = candidate1.mcTruthInfo() == 1; // flagD0Signal 'true' if candidate1 matched to D0 (particle)
254+
flagD0Reflection = candidate1.mcTruthInfo() == 2; // flagD0Reflection 'true' if candidate1, selected as D0 (particle), is matched to D0bar (antiparticle)
255+
256+
for (const auto& candidate2 : selectedCandidatesGroupedMC) {
294257
if (candidate2.isSelD0bar() < selectionFlagD0bar) { // discard candidates not selected as D0bar in inner loop
295258
continue;
296259
}
297-
flagD0barSignal = candidate2.flagMcMatchRec() == -(1 << aod::hf_cand_2prong::DecayType::D0ToPiK); // flagD0barSignal 'true' if candidate2 matched to D0bar (antiparticle)
298-
flagD0barReflection = candidate2.flagMcMatchRec() == 1 << aod::hf_cand_2prong::DecayType::D0ToPiK; // flagD0barReflection 'true' if candidate2, selected as D0bar (antiparticle), is matched to D0 (particle)
299-
if (yCandMax >= 0. && std::abs(hfHelper.yD0(candidate2)) > yCandMax) {
300-
continue;
301-
}
302-
if (ptCandMin >= 0. && candidate2.pt() < ptCandMin) {
260+
flagD0barSignal = candidate2.mcTruthInfo() == 2; // flagD0barSignal 'true' if candidate2 matched to D0bar (antiparticle)
261+
flagD0barReflection = candidate2.mcTruthInfo() == 1; // flagD0barReflection 'true' if candidate2, selected as D0bar (antiparticle), is matched to D0 (particle)
262+
263+
// Excluding trigger self-correlations (possible in case of both mass hypotheses of the same real particle, reconstructed as candidate1 for D0 and candidate2 for D0bar)
264+
if (candidate1.mRowIndex == candidate2.mRowIndex) { //this by definition should never happen, since each candidate is either D0 or D0bar
303265
continue;
304266
}
305-
// Excluding trigger self-correlations (possible in case of both mass hypotheses accepted)
306-
if (candidate1.mRowIndex == candidate2.mRowIndex) {
267+
if((candidate1.pt() - candidate2.pt()) < 1e-5 && (candidate1.eta() - candidate2.eta()) < 1e-5 && (candidate1.phi() - candidate2.phi()) < 1e-5) { //revised, temporary condition to avoid self-correlations (the best would be check the daughterIDs, but we don't store them at the moment)
307268
continue;
308-
}
269+
}
309270
// choice of options (D0/D0bar signal/bkg)
310271
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
311272
if (flagD0Signal) {
@@ -324,8 +285,8 @@ struct alice3correlatorddbar {
324285
candidate2.eta() - candidate1.eta(),
325286
candidate1.pt(),
326287
candidate2.pt());
327-
entryD0D0barRecoInfo(hfHelper.invMassD0ToPiK(candidate1),
328-
hfHelper.invMassD0barToKPi(candidate2),
288+
entryD0D0barRecoInfo(candidate1.m(), //mD0
289+
candidate2.m(), //mD0bar
329290
pairSignalStatus);
330291
double etaCut = 0.;
331292
double ptCut = 0.;
@@ -345,7 +306,7 @@ struct alice3correlatorddbar {
345306
}
346307

347308
PROCESS_SWITCH(alice3correlatorddbar, processMcRec, "Process MC Reco mode", false);
348-
*/
309+
349310
};
350311

351312
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)

0 commit comments

Comments
 (0)