Skip to content

Commit 964ccfb

Browse files
cterrevoalibuild
cterrevo
and
alibuild
authored
PWGHF: adding THnSparse for BDT and cutVariation to Xic task (#5580)
* adding THnSparse for BDT and cutVariation * fix to thnsparse axes * Please consider the following formatting changes --------- Co-authored-by: ALICE Action Bot <alibuild@cern.ch>
1 parent 5535e6d commit 964ccfb

File tree

1 file changed

+161
-21
lines changed

1 file changed

+161
-21
lines changed

‎PWGHF/D2H/Tasks/taskXic.cxx‎

Lines changed: 161 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -45,6 +45,18 @@ struct HfTaskXic {
4545
Configurable<float> dcaZTrackMax{"dcaZTrackMax", 0.0025, "max. DCAz for track"};
4646
Configurable<std::vector<double>> binsPt{"binsPt", std::vector<double>{hf_cuts_xic_to_p_k_pi::vecBinsPt}, "pT bin limits"};
4747

48+
// THnSparse for ML outputScores and Vars
49+
Configurable<bool> enableTHn{"enableTHn", false, "enable THn for Xic"};
50+
ConfigurableAxis thnConfigAxisPt{"thnConfigAxisPt", {36, 0, 36}, ""};
51+
ConfigurableAxis thnConfigAxisMass{"thnConfigAxisMass", {300, 1.98, 2.58}, ""};
52+
ConfigurableAxis thnConfigAxisPtProng{"thnConfigAxisPtProng", {100, 0, 20}, ""};
53+
ConfigurableAxis thnConfigAxisChi2PCA{"thnConfigAxisChi2PCA", {100, 0, 20}, ""};
54+
ConfigurableAxis thnConfigAxisDecLength{"thnConfigAxisDecLength", {10, 0, 0.05}, ""};
55+
ConfigurableAxis thnConfigAxisDecLengthXY{"thnConfigAxisDecLengthXY", {10, 0, 0.05}, ""};
56+
ConfigurableAxis thnConfigAxisCPA{"thnConfigAxisCPA", {20, 0.8, 1}, ""};
57+
ConfigurableAxis thnConfigAxisBdtScoreBkg{"thnConfigAxisBdtScoreBkg", {100, 0., 1.}, ""};
58+
ConfigurableAxis thnConfigAxisBdtScoreSignal{"thnConfigAxisBdtScoreSignal", {100, 0., 1.}, ""};
59+
//
4860
Service<o2::framework::O2DatabasePDG> pdg;
4961
HfHelper hfHelper;
5062

@@ -79,6 +91,11 @@ struct HfTaskXic {
7991

8092
void init(InitContext&)
8193
{
94+
std::array<bool, 5> doprocess{doprocessDataStd, doprocessDataWithMl, doprocessMcStd, doprocessMcWithMl};
95+
if ((std::accumulate(doprocess.begin(), doprocess.end(), 0)) != 1) {
96+
LOGP(fatal, "no or more than one process function enabled! Please check your configuration!");
97+
}
98+
8299
AxisSpec axisPPid = {100, 0.f, 10.0f, "#it{p} (GeV/#it{c})"};
83100
AxisSpec axisNSigmaPr = {100, -6.f, 6.f, "n#it{#sigma}_{p}"};
84101
AxisSpec axisNSigmaPi = {100, -6.f, 6.f, "n#it{#sigma}_{#pi}"};
@@ -176,7 +193,27 @@ struct HfTaskXic {
176193
registry.add("MC/reconstructed/background/hPtProng1RecBg", "3-prong candidates;prong 1 #it{p}_{T} (GeV/#it{c});;entries", {HistType::kTH2F, {{100, 0., 10.}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}});
177194
registry.add("MC/reconstructed/background/hPtProng2RecBg", "3-prong candidates;prong 2 #it{p}_{T} (GeV/#it{c});;entries", {HistType::kTH2F, {{100, 0., 10.}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}});
178195
registry.add("MC/reconstructed/background/hChi2PCARecBg", "3-prong candidates;prong Chi2PCA to sec. vertex (cm);; entries", {HistType::kTH2F, {{100, 0, 120}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}});
179-
}
196+
197+
// THn for candidates Xic+ cut variation
198+
if (enableTHn) {
199+
const AxisSpec thnAxisMass{thnConfigAxisMass, "inv. mass (p K #pi) (GeV/#it{c}^{2})"};
200+
const AxisSpec thnAxisPt{thnConfigAxisPt, "#it{p}_{T}(#Xi_{c}^{+}) (GeV/#it{c})"};
201+
const AxisSpec thnAxisChi2PCA{thnConfigAxisChi2PCA, "Chi2PCA to sec. vertex (cm)"};
202+
const AxisSpec thnAxisDecLength{thnConfigAxisDecLength, "decay length (cm)"};
203+
const AxisSpec thnAxisDecLengthXY{thnConfigAxisDecLengthXY, "decay length (cm)"};
204+
const AxisSpec thnAxisCPA{thnConfigAxisCPA, "cosine of pointing angle"};
205+
const AxisSpec thnAxisBdtScoreXicBkg{thnConfigAxisBdtScoreBkg, "BDT bkg score (Xic)"};
206+
const AxisSpec thnAxisBdtScoreXicPrompt{thnConfigAxisBdtScoreSignal, "BDT prompt score (Xic)"};
207+
const AxisSpec thnAxisBdtScoreXicNonPrompt{thnConfigAxisBdtScoreSignal, "BDT non-prompt score (Xic)"};
208+
const AxisSpec thnAxisMcOrigin{3, -0.5, 2.5, "MC origin"};
209+
210+
if (doprocessDataWithMl || doprocessMcWithMl) { // with ML
211+
registry.add("hnXicVarsWithBdt", "THn for Xic candidates with BDT scores", HistType::kTHnSparseF, {thnAxisMass, thnAxisPt, thnAxisBdtScoreXicBkg, thnAxisBdtScoreXicPrompt, thnAxisBdtScoreXicNonPrompt, thnAxisMcOrigin});
212+
} else {
213+
registry.add("hnXicVars", "THn for Xic candidates", HistType::kTHnSparseF, {thnAxisMass, thnAxisPt, thnAxisChi2PCA, thnAxisDecLength, thnAxisDecLengthXY, thnAxisCPA, thnAxisMcOrigin});
214+
}
215+
}
216+
} // end init
180217

181218
/// Selection of Xic daughters in geometrical acceptanc
182219
/// \param etaProng is the pseudorapidity of Xic prong
@@ -187,10 +224,10 @@ struct HfTaskXic {
187224
{
188225
return std::abs(etaProng) <= etaMaxAcceptance && ptProng >= ptMinAcceptance;
189226
}
190-
191-
void process(aod::Collision const& collision,
192-
TracksWPid const& tracks,
193-
soa::Filtered<soa::Join<aod::HfCand3Prong, aod::HfSelXicToPKPi>> const& candidates)
227+
template <bool useMl, typename Cands>
228+
void analysisData(aod::Collision const& collision,
229+
Cands const& candidates,
230+
TracksWPid const& tracks)
194231
{
195232
int nTracks = 0;
196233

@@ -207,15 +244,15 @@ struct HfTaskXic {
207244
}
208245

209246
registry.fill(HIST("Data/hMultiplicity"), nTracks); // filling the histo for multiplicity
247+
210248
for (const auto& candidate : candidates) {
211-
auto ptCandidate = candidate.pt();
212249
if (!(candidate.hfflag() & 1 << aod::hf_cand_3prong::DecayType::XicToPKPi)) {
213250
continue;
214251
}
215-
216252
if (yCandRecoMax >= 0. && std::abs(hfHelper.yXic(candidate)) > yCandRecoMax) {
217253
continue;
218254
}
255+
auto ptCandidate = candidate.pt();
219256

220257
if (candidate.isSelXicToPKPi() >= selectionFlagXic) { // pKpi
221258
registry.fill(HIST("Data/hMassVsPt"), hfHelper.invMassXicToPKPi(candidate), ptCandidate);
@@ -250,9 +287,9 @@ struct HfTaskXic {
250287
registry.fill(HIST("Data/hChi2PCA"), candidate.chi2PCA(), ptCandidate);
251288

252289
// PID histos
253-
const auto& trackProng0 = candidate.prong0_as<TracksWPid>(); // bachelor track
254-
const auto& trackProng1 = candidate.prong1_as<TracksWPid>(); // bachelor track
255-
const auto& trackProng2 = candidate.prong2_as<TracksWPid>(); // bachelor track
290+
auto trackProng0 = candidate.template prong0_as<TracksWPid>();
291+
auto trackProng1 = candidate.template prong1_as<TracksWPid>();
292+
auto trackProng2 = candidate.template prong2_as<TracksWPid>();
256293

257294
auto momentumProng0 = trackProng0.p();
258295
auto momentumProng1 = trackProng1.p();
@@ -283,26 +320,77 @@ struct HfTaskXic {
283320
registry.fill(HIST("Data/hPVsTOFNSigmaPr_Prong2"), momentumProng2, trackProng2.tofNSigmaPr());
284321
registry.fill(HIST("Data/hPVsTOFNSigmaPi_Prong2"), momentumProng2, trackProng2.tofNSigmaPi());
285322
registry.fill(HIST("Data/hPVsTOFNSigmaKa_Prong2"), momentumProng2, trackProng2.tofNSigmaKa());
286-
}
323+
324+
// THnSparse
325+
if (enableTHn) {
326+
double massXic(-1);
327+
double outputBkg(-1), outputPrompt(-1), outputFD(-1);
328+
if (candidate.isSelXicToPKPi() >= selectionFlagXic) {
329+
massXic = hfHelper.invMassXicToPKPi(candidate);
330+
if constexpr (useMl) {
331+
if (candidate.mlProbXicToPKPi().size() == 3) {
332+
outputBkg = candidate.mlProbXicToPKPi()[0]; /// bkg score
333+
outputPrompt = candidate.mlProbXicToPKPi()[1]; /// prompt score
334+
outputFD = candidate.mlProbXicToPKPi()[2]; /// non-prompt score
335+
}
336+
/// Fill the ML outputScores and variables of candidate Xic
337+
registry.get<THnSparse>(HIST("hnXicVarsWithBdt"))->Fill(massXic, ptCandidate, outputBkg, outputPrompt, outputFD, 0);
338+
} else {
339+
registry.get<THnSparse>(HIST("hnXicVars"))->Fill(massXic, ptCandidate, candidate.chi2PCA(), candidate.decayLength(), candidate.decayLengthXY(), candidate.cpa(), 0);
340+
}
341+
}
342+
if (candidate.isSelXicToPiKP() >= selectionFlagXic) {
343+
massXic = hfHelper.invMassXicToPiKP(candidate);
344+
if constexpr (useMl) {
345+
if (candidate.mlProbXicToPiKP().size() == 3) {
346+
outputBkg = candidate.mlProbXicToPiKP()[0]; /// bkg score
347+
outputPrompt = candidate.mlProbXicToPiKP()[1]; /// prompt score
348+
outputFD = candidate.mlProbXicToPiKP()[2]; /// non-prompt score
349+
}
350+
/// Fill the ML outputScores and variables of candidate
351+
registry.get<THnSparse>(HIST("hnXicVarsWithBdt"))->Fill(massXic, ptCandidate, outputBkg, outputPrompt, outputFD, 0);
352+
} else {
353+
registry.get<THnSparse>(HIST("hnXicVars"))->Fill(massXic, ptCandidate, candidate.chi2PCA(), candidate.decayLength(), candidate.decayLengthXY(), candidate.cpa(), 0);
354+
}
355+
}
356+
} // thn for Xic
357+
} // loop candidates
358+
} // end process data
359+
360+
void processDataStd(aod::Collision const& collision,
361+
soa::Filtered<soa::Join<aod::HfCand3Prong, aod::HfSelXicToPKPi>> const& candidates,
362+
TracksWPid const& tracks)
363+
{
364+
analysisData<false>(collision, candidates, tracks);
287365
}
288-
// Fill MC histograms
289-
void processMc(soa::Filtered<soa::Join<aod::HfCand3Prong, aod::HfSelXicToPKPi, aod::HfCand3ProngMcRec>> const& candidates,
290-
soa::Join<aod::McParticles, aod::HfCand3ProngMcGen> const& mcParticles)
366+
PROCESS_SWITCH(HfTaskXic, processDataStd, "Process Data with the standard method", true);
367+
368+
void processDataWithMl(aod::Collision const& collision,
369+
soa::Filtered<soa::Join<aod::HfCand3Prong, aod::HfSelXicToPKPi, aod::HfMlXicToPKPi>> const& candidatesMl, TracksWPid const& tracks)
291370
{
371+
analysisData<true>(collision, candidatesMl, tracks);
372+
}
373+
PROCESS_SWITCH(HfTaskXic, processDataWithMl, "Process Data with the ML method", false);
292374

375+
// Fill MC histograms
376+
template <bool useMl, typename Cands>
377+
void analysisMc(Cands const& candidates,
378+
soa::Join<aod::McParticles, aod::HfCand3ProngMcGen> const& mcParticles,
379+
aod::TracksWMc const&)
380+
{
293381
// MC rec.
294382
for (const auto& candidate : candidates) {
383+
auto ptCandidate = candidate.pt();
295384
// Selected Xic
296385
if (!(candidate.hfflag() & 1 << aod::hf_cand_3prong::DecayType::XicToPKPi)) {
297386
continue;
298-
} // rapidity selection
387+
}
299388
if (yCandRecoMax >= 0. && std::abs(hfHelper.yXic(candidate)) > yCandRecoMax) {
300389
continue;
301390
}
302391

303392
auto massXicToPKPi = 0.;
304393
auto massXicToPiKP = 0.;
305-
auto ptCandidate = candidate.pt();
306394

307395
if (candidate.isSelXicToPKPi() >= selectionFlagXic) {
308396
massXicToPKPi = hfHelper.invMassXicToPKPi(candidate);
@@ -312,8 +400,11 @@ struct HfTaskXic {
312400
}
313401

314402
if (std::abs(candidate.flagMcMatchRec()) == 1 << aod::hf_cand_3prong::DecayType::XicToPKPi) {
403+
// Get the corresponding MC particle.
404+
auto mcParticleProng0 = candidate.template prong0_as<aod::TracksWMc>().template mcParticle_as<soa::Join<aod::McParticles, aod::HfCand3ProngMcGen>>();
405+
auto pdgCodeProng0 = std::abs(mcParticleProng0.pdgCode());
315406
// Signal
316-
registry.fill(HIST("MC/reconstructed/signal/hPtRecSig"), ptCandidate); // rec. level pT
407+
registry.fill(HIST("MC/reconstructed/signal/hPtRecSig"), ptCandidate); // rec. level pT
317408

318409
if (massXicToPKPi != 0.) {
319410
registry.fill(HIST("MC/reconstructed/signal/hMassSig"), massXicToPKPi, ptCandidate);
@@ -339,6 +430,39 @@ struct HfTaskXic {
339430
registry.fill(HIST("MC/reconstructed/signal/hImpParErrSig"), candidate.errorImpactParameter0(), ptCandidate);
340431
registry.fill(HIST("MC/reconstructed/signal/hDecLenErrSig"), candidate.errorDecayLength(), ptCandidate);
341432
registry.fill(HIST("MC/reconstructed/signal/hChi2PCARecSig"), candidate.chi2PCA(), ptCandidate);
433+
434+
if (enableTHn) {
435+
double massXic(-1);
436+
double outputBkg(-1), outputPrompt(-1), outputFD(-1);
437+
if ((candidate.isSelXicToPKPi() >= selectionFlagXic) && pdgCodeProng0 == kProton) {
438+
massXic = hfHelper.invMassXicToPKPi(candidate);
439+
if constexpr (useMl) {
440+
if (candidate.mlProbXicToPKPi().size() == 3) {
441+
outputBkg = candidate.mlProbXicToPKPi()[0]; /// bkg score
442+
outputPrompt = candidate.mlProbXicToPKPi()[1]; /// prompt score
443+
outputFD = candidate.mlProbXicToPKPi()[2]; /// non-prompt score
444+
}
445+
/// Fill the ML outputScores and variables of candidate (todo: add multiplicity)
446+
registry.get<THnSparse>(HIST("hnXicVarsWithBdt"))->Fill(massXic, ptCandidate, outputBkg, outputPrompt, outputFD, candidate.originMcRec());
447+
} else {
448+
registry.get<THnSparse>(HIST("hnXicVars"))->Fill(massXic, ptCandidate, candidate.chi2PCA(), candidate.decayLength(), candidate.decayLengthXY(), candidate.cpa(), candidate.originMcRec());
449+
}
450+
}
451+
if ((candidate.isSelXicToPiKP() >= selectionFlagXic) && pdgCodeProng0 == kPiPlus) {
452+
massXic = hfHelper.invMassXicToPiKP(candidate);
453+
if constexpr (useMl) {
454+
if (candidate.mlProbXicToPiKP().size() == 3) {
455+
outputBkg = candidate.mlProbXicToPiKP()[0]; /// bkg score
456+
outputPrompt = candidate.mlProbXicToPiKP()[1]; /// prompt score
457+
outputFD = candidate.mlProbXicToPiKP()[2]; /// non-prompt score
458+
}
459+
/// Fill the ML outputScores and variables of candidate (todo: add multiplicity)
460+
registry.get<THnSparse>(HIST("hnXicVarsWithBdt"))->Fill(massXic, ptCandidate, outputBkg, outputPrompt, outputFD, candidate.originMcRec());
461+
} else {
462+
registry.get<THnSparse>(HIST("hnXicVars"))->Fill(massXic, ptCandidate, candidate.chi2PCA(), candidate.decayLength(), candidate.decayLengthXY(), candidate.cpa(), candidate.originMcRec());
463+
}
464+
}
465+
} // enable THn
342466
} else {
343467
// Background
344468
registry.fill(HIST("MC/reconstructed/background/hPtRecBg"), ptCandidate);
@@ -367,8 +491,9 @@ struct HfTaskXic {
367491
registry.fill(HIST("MC/reconstructed/background/hImpParErrBg"), candidate.errorImpactParameter0(), ptCandidate);
368492
registry.fill(HIST("MC/reconstructed/background/hDecLenErrBg"), candidate.errorDecayLength(), ptCandidate);
369493
registry.fill(HIST("MC/reconstructed/background/hChi2PCARecBg"), candidate.chi2PCA(), ptCandidate);
370-
}
371-
}
494+
} // Xic background
495+
} // candidate loop
496+
372497
// MC gen.
373498
for (const auto& particle : mcParticles) {
374499
if (std::abs(particle.flagMcMatchGen()) == 1 << aod::hf_cand_3prong::DecayType::XicToPKPi) {
@@ -381,7 +506,7 @@ struct HfTaskXic {
381506
std::array<float, 3> yProngs;
382507
std::array<float, 3> etaProngs;
383508
int counter = 0;
384-
for (const auto& daught : particle.daughters_as<aod::McParticles>()) {
509+
for (const auto& daught : particle.daughters_as<soa::Join<aod::McParticles, aod::HfCand3ProngMcGen>>()) {
385510
ptProngs[counter] = daught.pt();
386511
etaProngs[counter] = daught.eta();
387512
yProngs[counter] = RecoDecay::y(std::array{daught.px(), daught.py(), daught.pz()}, pdg->Mass(daught.pdgCode()));
@@ -401,7 +526,22 @@ struct HfTaskXic {
401526
}
402527
}
403528
}
404-
PROCESS_SWITCH(HfTaskXic, processMc, "Process MC", false);
529+
530+
void processMcStd(soa::Filtered<soa::Join<aod::HfCand3Prong, aod::HfSelXicToPKPi, aod::HfCand3ProngMcRec>> const& selectedCandidatesMc,
531+
soa::Join<aod::McParticles, aod::HfCand3ProngMcGen> const& mcParticles,
532+
aod::TracksWMc const& tracksWithMc)
533+
{
534+
analysisMc<false>(selectedCandidatesMc, mcParticles, tracksWithMc);
535+
}
536+
PROCESS_SWITCH(HfTaskXic, processMcStd, "Process MC with the standard method", false);
537+
538+
void processMcWithMl(soa::Filtered<soa::Join<aod::HfCand3Prong, aod::HfSelXicToPKPi, aod::HfMlXicToPKPi, aod::HfCand3ProngMcRec>> const& selectedCandidatesMlMc,
539+
soa::Join<aod::McParticles, aod::HfCand3ProngMcGen> const& mcParticles,
540+
aod::TracksWMc const& tracksWithMc)
541+
{
542+
analysisMc<true>(selectedCandidatesMlMc, mcParticles, tracksWithMc);
543+
}
544+
PROCESS_SWITCH(HfTaskXic, processMcWithMl, "Process Mc with the ML method", false);
405545
};
406546

407547
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)

0 commit comments

Comments
 (0)