Skip to content

Commit 23e5a7d

Browse files
njacazionjacazio
authored
[ALICE3] Fix 3 prong decays (#14548)
1 parent 333e4ff commit 23e5a7d

File tree

3 files changed

+84
-38
lines changed

3 files changed

+84
-38
lines changed

‎ALICE3/TableProducer/alice3-decayfinder.cxx‎

Lines changed: 17 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -661,7 +661,7 @@ struct alice3decayFinder {
661661
if (doprocessFindLc) {
662662
for (auto const& mcParticle : mcParticles) {
663663
if (std::abs(mcParticle.pdgCode()) != motherPdgCode) {
664-
mcGenFlags(-1, -1, -1);
664+
mcGenFlags(-1, 0, 0);
665665
continue;
666666
}
667667
std::vector<int> idxBhadMothers{};
@@ -671,7 +671,7 @@ struct alice3decayFinder {
671671
auto bHadMother = mcParticles.rawIteratorAt(idxBhadMothers[0]);
672672
ptBMotherGen = bHadMother.pt();
673673
}
674-
mcGenFlags(origin, ptBMotherGen, mcParticle.pdgCode() ? charmHadFlag : -charmHadFlag);
674+
mcGenFlags(origin, ptBMotherGen, mcParticle.pdgCode() > 0 ? charmHadFlag : -charmHadFlag);
675675
if (mcParticle.pdgCode() > 0) {
676676
histos.fill(HIST("h2dGen3Prong"), mcParticle.pt(), mcParticle.eta());
677677
} else {
@@ -1040,9 +1040,9 @@ struct alice3decayFinder {
10401040
mCandidate3Prong.eta,
10411041
mCandidate3Prong.phi,
10421042
mCandidate3Prong.pt,
1043-
mCandidate3Prong.Pdaug2[0], mCandidate3Prong.Pdaug2[1], mCandidate3Prong.Pdaug2[2],
1044-
mCandidate3Prong.Pdaug1[0], mCandidate3Prong.Pdaug1[1], mCandidate3Prong.Pdaug1[2],
10451043
mCandidate3Prong.Pdaug0[0], mCandidate3Prong.Pdaug0[1], mCandidate3Prong.Pdaug0[2],
1044+
mCandidate3Prong.Pdaug1[0], mCandidate3Prong.Pdaug1[1], mCandidate3Prong.Pdaug1[2],
1045+
mCandidate3Prong.Pdaug2[0], mCandidate3Prong.Pdaug2[1], mCandidate3Prong.Pdaug2[2],
10461046
mCandidate3Prong.impactParameterY0, mCandidate3Prong.impactParameterY1, mCandidate3Prong.impactParameterY2,
10471047
std::sqrt(mCandidate3Prong.errorImpactParameterY0),
10481048
std::sqrt(mCandidate3Prong.errorImpactParameterY1),
@@ -1066,15 +1066,27 @@ struct alice3decayFinder {
10661066

10671067
void processFindLc(aod::Collision const& collision,
10681068
aod::McParticles const& mcParticles,
1069-
Alice3TracksWPid const&)
1069+
Alice3TracksWPid const& tracks)
10701070
{
1071+
LOG(debug) << "Processing Lc candidates for collision " << collision.globalIndex() << " with " << tracks.size() << " tracks";
1072+
for (auto const& track : tracks) {
1073+
if (track.has_mcParticle()) {
1074+
LOG(debug) << " - track index: " << track.globalIndex() << ", pT: " << track.pt() << " (MC pt " << track.mcParticle().pt() << "), PID: " << track.mcParticle().pdgCode();
1075+
}
1076+
}
10711077

10721078
auto tracksPiPlus = tracksPiPlusFromLc->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache);
1079+
LOG(debug) << " - found " << tracksPiPlus.size() << " pi+ from Lc";
10731080
auto tracksKaPlus = tracksKaPlusFromLc->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache);
1081+
LOG(debug) << " - found " << tracksKaPlus.size() << " K+ from Lc";
10741082
auto tracksPrPlus = tracksPrPlusFromLc->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache);
1083+
LOG(debug) << " - found " << tracksPrPlus.size() << " p+ from Lc";
10751084
auto tracksPiMinus = tracksPiMinusFromLc->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache);
1085+
LOG(debug) << " - found " << tracksPiMinus.size() << " pi- from Lc";
10761086
auto tracksKaMinus = tracksKaMinusFromLc->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache);
1087+
LOG(debug) << " - found " << tracksKaMinus.size() << " K- from Lc";
10771088
auto tracksPrMinus = tracksPrMinusFromLc->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache);
1089+
LOG(debug) << " - found " << tracksPrMinus.size() << " p- from Lc";
10781090

10791091
if (doDCAplots3Prong) {
10801092
for (auto const& track : tracksPiPlus)

‎ALICE3/TableProducer/alice3HfSelector3Prong.cxx‎

Lines changed: 9 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -313,8 +313,6 @@ struct Alice3HfSelector3Prong {
313313
template <CharmHadAlice3 CharmHad, typename CandType>
314314
void runSelect3Prong(CandType const& cands)
315315
{
316-
bool isSelMassHypo0{false};
317-
bool isSelMassHypo1{false};
318316
std::vector<float> outputMl{-1.f, -1.f, -1.f};
319317
uint32_t pidMask = 0;
320318

@@ -324,21 +322,21 @@ struct Alice3HfSelector3Prong {
324322
outputMl = {-1.f, -1.f, -1.f};
325323
pidMask = 0;
326324

327-
auto ptCand = cand.pt();
328-
int const ptBin = findBin(binsPt, ptCand);
325+
const float ptCand = cand.pt();
326+
const int ptBin = findBin(binsPt, ptCand);
329327
if (ptBin == -1) {
330-
candSelFlags(isSelMassHypo0, isSelMassHypo1, pidMask);
328+
candSelFlags(false, false, pidMask);
331329
if (applyMl) {
332330
candMlScores(outputMl[0], outputMl[1], outputMl[2]);
333331
}
334332
continue;
335333
}
336334

337335
// Here all cands pass the cut on the mass selection
338-
bool const selMassHypo0 = selectionCandidateMass<CharmHad, false>(ptBin, cand);
339-
bool const selMassHypo1 = selectionCandidateMass<CharmHad, true>(ptBin, cand);
336+
const bool selMassHypo0 = selectionCandidateMass<CharmHad, false>(ptBin, cand);
337+
const bool selMassHypo1 = selectionCandidateMass<CharmHad, true>(ptBin, cand);
340338
if (!selMassHypo0 && !selMassHypo1) {
341-
candSelFlags(isSelMassHypo0, isSelMassHypo1, pidMask);
339+
candSelFlags(false, false, pidMask);
342340
if (applyMl) {
343341
candMlScores(outputMl[0], outputMl[1], outputMl[2]);
344342
}
@@ -348,7 +346,7 @@ struct Alice3HfSelector3Prong {
348346

349347
// Topological selection (TODO: track quality selection)
350348
if (!selectionTopol(cand, ptCand)) {
351-
candSelFlags(isSelMassHypo0, isSelMassHypo1, pidMask);
349+
candSelFlags(false, false, pidMask);
352350
if (applyMl) {
353351
candMlScores(outputMl[0], outputMl[1], outputMl[2]);
354352
}
@@ -359,7 +357,7 @@ struct Alice3HfSelector3Prong {
359357
// PID selection
360358
configurePidMask<CharmHad>(cand, pidMask);
361359
if (pidMask == 0) {
362-
candSelFlags(isSelMassHypo0, isSelMassHypo1, pidMask);
360+
candSelFlags(false, false, pidMask);
363361
if (applyMl) {
364362
candMlScores(outputMl[0], outputMl[1], outputMl[2]);
365363
}
@@ -375,7 +373,7 @@ struct Alice3HfSelector3Prong {
375373
isSelectedMl = mlResponse.isSelectedMl(inputFeaturesMassHypo0, ptCand, outputMl);
376374
candMlScores(outputMl[0], outputMl[1], outputMl[2]);
377375
if (!isSelectedMl) {
378-
candSelFlags(isSelMassHypo0, isSelMassHypo1, pidMask);
376+
candSelFlags(false, false, pidMask);
379377
continue;
380378
}
381379

‎ALICE3/Tasks/alice3HfTask3Prong.cxx‎

Lines changed: 58 additions & 22 deletions
Original file line numberDiff line numberDiff line change
@@ -87,14 +87,15 @@ struct Alice3HfTask3Prong {
8787
HistogramRegistry registry{"registry", {}};
8888

8989
// Names of folders and suffixes for MC signal histograms
90-
constexpr static std::string_view SignalFolders[] = {"signal", "prompt", "nonprompt"};
91-
constexpr static std::string_view SignalSuffixes[] = {"", "Prompt", "NonPrompt"};
92-
9390
enum SignalClasses : int {
9491
Signal = 0,
9592
Prompt,
96-
NonPrompt
93+
NonPrompt,
94+
Bkg,
95+
NSignalClasses
9796
};
97+
constexpr static std::string_view SignalFolders[SignalClasses::NSignalClasses] = {"signal", "prompt", "nonprompt", "background"};
98+
constexpr static std::string_view SignalSuffixes[SignalClasses::NSignalClasses] = {"", "Prompt", "NonPrompt", "Bkg"};
9899

99100
void init(InitContext&)
100101
{
@@ -109,9 +110,11 @@ struct Alice3HfTask3Prong {
109110
}
110111

111112
auto addHistogramsRec = [&](const std::string& histoName, const std::string& xAxisTitle, const std::string& yAxisTitle, const HistogramConfigSpec& configSpec) {
112-
registry.add(("MC/rec/signal/" + histoName + "RecSig").c_str(), ("3-prong cands (matched);" + xAxisTitle + ";" + yAxisTitle).c_str(), configSpec);
113-
registry.add(("MC/rec/prompt/" + histoName + "RecSigPrompt").c_str(), ("3-prong cands (matched, prompt);" + xAxisTitle + ";" + yAxisTitle).c_str(), configSpec);
114-
registry.add(("MC/rec/nonprompt/" + histoName + "RecSigNonPrompt").c_str(), ("3-prong cands (matched, non-prompt);" + xAxisTitle + ";" + yAxisTitle).c_str(), configSpec);
113+
const char* basePath = "MC/rec";
114+
registry.add(Form("%s/signal/%sRecSig%s", basePath, histoName.c_str(), SignalSuffixes[SignalClasses::Signal].data()), ("3-prong cands (matched);" + xAxisTitle + ";" + yAxisTitle).c_str(), configSpec);
115+
registry.add(Form("%s/prompt/%sRecSig%s", basePath, histoName.c_str(), SignalSuffixes[SignalClasses::Prompt].data()), ("3-prong cands (matched, prompt);" + xAxisTitle + ";" + yAxisTitle).c_str(), configSpec);
116+
registry.add(Form("%s/nonprompt/%sRecSig%s", basePath, histoName.c_str(), SignalSuffixes[SignalClasses::NonPrompt].data()), ("3-prong cands (matched, non-prompt);" + xAxisTitle + ";" + yAxisTitle).c_str(), configSpec);
117+
registry.add(Form("%s/background/%sRecSig%s", basePath, histoName.c_str(), SignalSuffixes[SignalClasses::Bkg].data()), ("3-prong cands (unmatched);" + xAxisTitle + ";" + yAxisTitle).c_str(), configSpec);
115118
};
116119

117120
auto addHistogramsGen = [&](const std::string& histoName, const std::string& xAxisTitle, const std::string& yAxisTitle, const HistogramConfigSpec& configSpec) {
@@ -174,9 +177,16 @@ struct Alice3HfTask3Prong {
174177
auto h2 = registry.add<TH2>("hSelectionStatus", "3-prong cands;selection status;entries", {HistType::kTH2F, {{5, -0.5, 4.5}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}});
175178
h2->GetXaxis()->SetBinLabel(1, "mass hypo 0");
176179
h2->GetXaxis()->SetBinLabel(2, "mass hypo 1");
177-
auto h = registry.add<TH1>("MC/rec/hCandidateCounter", "Candidate counter;entries", {HistType::kTH1D, {{2, -0.5, 1.5}}});
180+
auto h = registry.add<TH1>("MC/rec/hCandidateCounter", "Candidate counter;entries", {HistType::kTH1D, {{4, -0.5, 3.5}}});
178181
h->GetXaxis()->SetBinLabel(1, "Calls");
179182
h->GetXaxis()->SetBinLabel(2, "Candidates");
183+
h->GetXaxis()->SetBinLabel(3, "Passed Y cut");
184+
h->GetXaxis()->SetBinLabel(4, "Has MC match");
185+
186+
registry.add("MC/rec/hPtDeltaProng0", ";prong 0 (#it{p}_{T}-#it{p}_{T, gen}) (GeV/#it{c});entries", {HistType::kTH1F, {{100, -5, 5.}}});
187+
registry.add("MC/rec/hPtDeltaProng1", ";prong 1 (#it{p}_{T}-#it{p}_{T, gen}) (GeV/#it{c});entries", {HistType::kTH1F, {{100, -5, 5.}}});
188+
registry.add("MC/rec/hPtDeltaProng2", ";prong 2 (#it{p}_{T}-#it{p}_{T, gen}) (GeV/#it{c});entries", {HistType::kTH1F, {{100, -5, 5.}}});
189+
180190
// Number of events processed
181191
h = registry.add<TH1>("hNEventsProcessed", "number of events processed;entries;", {HistType::kTH1F, {{2, 0.5, 2.5}}});
182192
h->GetXaxis()->SetBinLabel(1, "Generated");
@@ -207,7 +217,7 @@ struct Alice3HfTask3Prong {
207217
/// Helper function for filling MC reconstructed histograms for prompt, nonpromt and common (signal)
208218
/// \param candidate is a reconstructed candidate
209219
/// \tparam SignalType is an enum defining which histogram in which folder (signal, prompt or nonpromt) to fill
210-
template <CharmHadAlice3 CharmHad, int SignalType, typename CandidateType>
220+
template <CharmHadAlice3 CharmHad, SignalClasses SignalType, typename CandidateType>
211221
void fillHistogramsRecSig(CandidateType const& candidate, float mass, bool isSwapped = false)
212222
{
213223
static constexpr auto histoPrefix = HIST("MC/rec/") + HIST(SignalFolders[SignalType]) + HIST("/");
@@ -268,33 +278,50 @@ struct Alice3HfTask3Prong {
268278
if (yCandRecoMax >= 0. && std::abs(hfHelper.getCandY<CharmHad>(candidate)) > yCandRecoMax) {
269279
continue;
270280
}
271-
auto mcParticle = allParticles.iteratorAt(candidate.particleMcRec());
272-
if (candidate.particleMcRec() > 0) {
281+
registry.fill(HIST("MC/rec/hCandidateCounter"), 2.);
282+
if (candidate.particleMcRec() >= 0) {
283+
registry.fill(HIST("MC/rec/hCandidateCounter"), 3.);
284+
auto mcParticle = allParticles.iteratorAt(candidate.particleMcRec());
273285
if (mcParticle.has_daughters()) {
274286
auto daughters = mcParticle.daughtersIds();
275-
LOG(info) << "Reco candidate matched to MC particle with PDG " << mcParticle.pdgCode() << " daughters: " << daughters.size();
276-
for (auto dauId : daughters) {
287+
LOG(debug) << "Reco candidate matched to MC particle with PDG " << mcParticle.pdgCode() << " daughters: " << daughters.size();
288+
int prongIdx = 0;
289+
for (int dauId = daughters[0]; dauId <= daughters[1]; dauId++) {
277290
auto dau = allParticles.iteratorAt(dauId);
278-
LOG(info) << " dauId: " << dauId << " PDG: " << dau.pdgCode();
291+
LOG(debug) << " dauId: " << dauId << " PDG: " << dau.pdgCode() << " with pT: " << dau.pt();
292+
switch (prongIdx) {
293+
case 0:
294+
registry.fill(HIST("MC/rec/hPtDeltaProng0"), candidate.ptProng0() - dau.pt());
295+
break;
296+
case 1:
297+
registry.fill(HIST("MC/rec/hPtDeltaProng1"), candidate.ptProng1() - dau.pt());
298+
break;
299+
case 2:
300+
registry.fill(HIST("MC/rec/hPtDeltaProng2"), candidate.ptProng2() - dau.pt());
301+
break;
302+
default:
303+
break;
304+
}
305+
prongIdx++;
279306
}
280307
}
281308
}
282309

283-
if (candidate.flagMcRec() != 0) {
284-
// Get the corresponding MC particle.
310+
if (candidate.flagMcRec() != 0) { // Particle is matched to MC truth
285311

312+
// Get the corresponding MC particle.
286313
const auto pt = candidate.pt();
287314
const auto originType = candidate.originMcRec();
288315

289316
if (candidate.isSelMassHypo0()) {
290317
registry.fill(HIST("hSelectionStatus"), 0., pt);
291318
double mass = hfHelper.getCandMass<CharmHad, false>(candidate);
292319
/// Fill histograms
293-
fillHistogramsRecSig<CharmHad, Signal>(candidate, mass, false);
320+
fillHistogramsRecSig<CharmHad, SignalClasses::Signal>(candidate, mass, false);
294321
if (originType == RecoDecay::OriginType::Prompt) {
295-
fillHistogramsRecSig<CharmHad, Prompt>(candidate, mass, false);
322+
fillHistogramsRecSig<CharmHad, SignalClasses::Prompt>(candidate, mass, false);
296323
} else if (originType == RecoDecay::OriginType::NonPrompt) {
297-
fillHistogramsRecSig<CharmHad, NonPrompt>(candidate, mass, false);
324+
fillHistogramsRecSig<CharmHad, SignalClasses::NonPrompt>(candidate, mass, false);
298325
}
299326
if (fillThn) {
300327
std::vector<double> valuesToFill{mass, pt};
@@ -312,11 +339,11 @@ struct Alice3HfTask3Prong {
312339
registry.fill(HIST("hSelectionStatus"), 1., pt);
313340
double mass = hfHelper.getCandMass<CharmHad, true>(candidate);
314341
/// Fill histograms
315-
fillHistogramsRecSig<CharmHad, Signal>(candidate, mass, true);
342+
fillHistogramsRecSig<CharmHad, SignalClasses::Signal>(candidate, mass, true);
316343
if (originType == RecoDecay::OriginType::Prompt) {
317-
fillHistogramsRecSig<CharmHad, Prompt>(candidate, mass, true);
344+
fillHistogramsRecSig<CharmHad, SignalClasses::Prompt>(candidate, mass, true);
318345
} else if (originType == RecoDecay::OriginType::NonPrompt) {
319-
fillHistogramsRecSig<CharmHad, NonPrompt>(candidate, mass, true);
346+
fillHistogramsRecSig<CharmHad, SignalClasses::NonPrompt>(candidate, mass, true);
320347
}
321348
if (fillThn) {
322349
std::vector<double> valuesToFill{mass, pt};
@@ -330,6 +357,15 @@ struct Alice3HfTask3Prong {
330357
registry.get<THnSparse>(HIST("hSparseRec"))->Fill(valuesToFill.data());
331358
}
332359
}
360+
} else { // Background
361+
if (candidate.isSelMassHypo0()) {
362+
double mass = hfHelper.getCandMass<CharmHad, false>(candidate);
363+
fillHistogramsRecSig<CharmHad, SignalClasses::Bkg>(candidate, mass, false);
364+
}
365+
if (candidate.isSelMassHypo1()) {
366+
double mass = hfHelper.getCandMass<CharmHad, true>(candidate);
367+
fillHistogramsRecSig<CharmHad, SignalClasses::Bkg>(candidate, mass, true);
368+
}
333369
}
334370
}
335371
}

0 commit comments

Comments
 (0)