Skip to content

Commit b0a9d21

Browse files
authored
[PWGHF] KFParticle based 3 prong decay reconstruction development (AliceO2Group#10610)
1 parent b2e8973 commit b0a9d21

File tree

4 files changed

+482
-162
lines changed

4 files changed

+482
-162
lines changed

PWGHF/Core/SelectorCuts.h

Lines changed: 15 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -465,6 +465,7 @@ namespace hf_cuts_lc_to_p_k_pi
465465
{
466466
static constexpr int NBinsPt = 10;
467467
static constexpr int NCutVars = 11;
468+
static constexpr int NCutKfVars = 12;
468469
// default values for the pT bin edges (can be used to configure histogram axis)
469470
// offset by 1 from the bin numbers in cuts array
470471
constexpr double BinsPt[NBinsPt + 1] = {
@@ -493,6 +494,19 @@ constexpr double Cuts[NBinsPt][NCutVars] = {{0.4, 0.4, 0.4, 0.4, 0., 0.005, 0.,
493494
{0.4, 0.4, 0.4, 0.4, 0., 0.005, 0., 0., 0., 1e+10, -1.}, /* 12 < pT < 24 */
494495
{0.4, 0.4, 0.4, 0.4, 0., 0.005, 0., 0., 0., 1e+10, -1.}}; /* 24 < pT < 36 */
495496

497+
// default value for the cuts Chi2Prim Chi2Geo DCA, cm Chi2Geo Chi2Topo
498+
// P K Pi KPi PPi PK KPi PPi PK ↓ LdL ↓
499+
constexpr double CutsKf[NBinsPt][NCutKfVars] = {{3., 3., 3., 3., 3., 3., 0.01, 0.01, 0.01, 3., 5., 5.}, /* 0 < pT < 1 */
500+
{3., 3., 3., 3., 3., 3., 0.01, 0.01, 0.01, 3., 5., 5.}, /* 1 < pT < 2 */
501+
{3., 3., 3., 3., 3., 3., 0.01, 0.01, 0.01, 3., 5., 5.}, /* 2 < pT < 3 */
502+
{3., 3., 3., 3., 3., 3., 0.01, 0.01, 0.01, 3., 5., 5.}, /* 3 < pT < 4 */
503+
{3., 3., 3., 3., 3., 3., 0.01, 0.01, 0.01, 3., 5., 5.}, /* 4 < pT < 5 */
504+
{3., 3., 3., 3., 3., 3., 0.01, 0.01, 0.01, 3., 5., 5.}, /* 5 < pT < 6 */
505+
{3., 3., 3., 3., 3., 3., 0.01, 0.01, 0.01, 3., 5., 5.}, /* 6 < pT < 8 */
506+
{3., 3., 3., 3., 3., 3., 0.01, 0.01, 0.01, 3., 5., 5.}, /* 8 < pT < 12 */
507+
{3., 3., 3., 3., 3., 3., 0.01, 0.01, 0.01, 3., 5., 5.}, /* 12 < pT < 24 */
508+
{3., 3., 3., 3., 3., 3., 0.01, 0.01, 0.01, 3., 5., 5.}}; /* 24 < pT < 36 */
509+
496510
// row labels
497511
static const std::vector<std::string> labelsPt = {
498512
"pT bin 0",
@@ -508,6 +522,7 @@ static const std::vector<std::string> labelsPt = {
508522

509523
// column labels
510524
static const std::vector<std::string> labelsCutVar = {"m", "pT p", "pT K", "pT Pi", "Chi2PCA", "decay length", "cos pointing angle", "decLengthXY", "normDecLXY", "impParXY", "mass (Kpi)"};
525+
static const std::vector<std::string> labelsCutKfVar = {"kfChi2PrimPr", "kfChi2PrimKa", "kfChi2PrimPi", "kfChi2GeoKaPi", "kfChi2GeoPrPi", "kfChi2GeoPrKa", "kfDcaKaPi", "kfDcaPrPi", "kfDcaPrKa", "kfChi2Geo", "kfLdL", "kfChi2Topo"};
511526
} // namespace hf_cuts_lc_to_p_k_pi
512527

513528
namespace hf_cuts_lc_to_k0s_p

PWGHF/TableProducer/candidateCreator3Prong.cxx

Lines changed: 28 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -85,6 +85,9 @@ struct HfCandidateCreator3Prong {
8585
Configurable<bool> createDs{"createDs", false, "enable Ds+/- candidate creation"};
8686
Configurable<bool> createLc{"createLc", false, "enable Lc+/- candidate creation"};
8787
Configurable<bool> createXic{"createXic", false, "enable Xic+/- candidate creation"};
88+
// KF
89+
Configurable<bool> applyTopoConstraint{"applyTopoConstraint", false, "apply origin from PV hypothesis for created candidate, works only in KF mode"};
90+
Configurable<bool> applyInvMassConstraint{"applyInvMassConstraint", false, "apply particle type hypothesis to recalculate created candidate's momentum, works only in KF mode"};
8891

8992
HfEventSelection hfEvSel; // event selection and monitoring
9093
o2::vertexing::DCAFitterN<3> df; // 3-prong vertex fitter
@@ -149,6 +152,10 @@ struct HfCandidateCreator3Prong {
149152
LOGP(fatal, "At least one particle specie should be enabled for the creation.");
150153
}
151154

155+
if (createLc && createXic && applyInvMassConstraint) {
156+
LOGP(fatal, "Unable to apply invariant mass constraint due to ambiguity of mass hypothesis: only one of Lc and Xic can be reconstructed.");
157+
}
158+
152159
// histograms
153160
registry.add("hMass3PKPi", "3-prong candidates;inv. mass (pK#pi) (GeV/#it{c}^{2});entries", {HistType::kTH1D, {{1200, 1.8, 3.0}}});
154161
registry.add("hMass3PiKP", "3-prong candidates;inv. mass (#pi Kp) (GeV/#it{c}^{2});entries", {HistType::kTH1D, {{1200, 1.8, 3.0}}});
@@ -487,6 +494,15 @@ struct HfCandidateCreator3Prong {
487494
kfCandPiKK.SetConstructMethod(2);
488495
kfCandPiKK.Construct(kfDaughtersPiKK, 3);
489496

497+
const float chi2topo = kfCalculateChi2ToPrimaryVertex(kfCandPKPi, kfpV);
498+
499+
if (applyTopoConstraint) { // constraints applied after chi2topo getter - to preserve unbiased value of chi2topo
500+
for (const auto& kfCand : std::array<KFParticle*, 5>{&kfCandPKPi, &kfCandPiKP, &kfCandPiKPi, &kfCandKKPi, &kfCandPiKK}) {
501+
kfCand->SetProductionVertex(kfpV);
502+
kfCand->TransportToDecayVertex();
503+
}
504+
}
505+
490506
KFParticle kfPairKPi;
491507
const KFParticle* kfDaughtersKPi[3] = {&kfSecondKaon, &kfThirdPion};
492508
kfPairKPi.SetConstructMethod(2);
@@ -505,8 +521,15 @@ struct HfCandidateCreator3Prong {
505521
const float massKPi = kfPairKPi.GetMass();
506522
const float massPiK = kfPairPiK.GetMass();
507523

524+
if (applyInvMassConstraint) { // constraints applied after minv getters - to preserve unbiased values of minv
525+
kfCandPKPi.SetNonlinearMassConstraint(createLc ? MassLambdaCPlus : MassXiCPlus);
526+
kfCandPiKP.SetNonlinearMassConstraint(createLc ? MassLambdaCPlus : MassXiCPlus);
527+
kfCandPiKPi.SetNonlinearMassConstraint(MassDPlus);
528+
kfCandKKPi.SetNonlinearMassConstraint(MassDS);
529+
kfCandPiKK.SetNonlinearMassConstraint(MassDS);
530+
}
531+
508532
const float chi2geo = kfCandPKPi.Chi2() / kfCandPKPi.NDF();
509-
const float chi2topo = kfCalculateChi2ToPrimaryVertex(kfCandPKPi, kfpV);
510533
const std::pair<float, float> ldl = kfCalculateLdL(kfCandPKPi, kfpV);
511534

512535
std::array<float, 3> pProng0 = kfCalculateProngMomentumInSecondaryVertex(kfFirstProton, kfCandPiKP);
@@ -787,6 +810,8 @@ struct HfCandidateCreator3ProngExpressions {
787810
Configurable<bool> matchKinkedDecayTopology{"matchKinkedDecayTopology", false, "Match also candidates with tracks that decay with kinked topology"};
788811
Configurable<bool> matchInteractionsWithMaterial{"matchInteractionsWithMaterial", false, "Match also candidates with tracks that interact with material"};
789812

813+
constexpr static std::size_t NDaughtersResonant{2u};
814+
790815
HfEventSelectionMc hfEvSelMc; // mc event selection and monitoring
791816

792817
using BCsInfo = soa::Join<aod::BCs, aod::Timestamps, aod::BcSels>;
@@ -924,7 +949,7 @@ struct HfCandidateCreator3ProngExpressions {
924949
swapping = int8_t(std::abs(arrayDaughters[0].mcParticle().pdgCode()) == kPiPlus);
925950
}
926951
RecoDecay::getDaughters(mcParticles.rawIteratorAt(indexRec), &arrDaughIndex, std::array{0}, 1);
927-
if (arrDaughIndex.size() == 2) {
952+
if (arrDaughIndex.size() == NDaughtersResonant) {
928953
for (auto iProng = 0u; iProng < arrDaughIndex.size(); ++iProng) {
929954
auto daughI = mcParticles.rawIteratorAt(arrDaughIndex[iProng]);
930955
arrPDGDaugh[iProng] = std::abs(daughI.pdgCode());
@@ -970,7 +995,7 @@ struct HfCandidateCreator3ProngExpressions {
970995
swapping = int8_t(std::abs(arrayDaughters[0].mcParticle().pdgCode()) == kPiPlus);
971996
}
972997
RecoDecay::getDaughters(mcParticles.rawIteratorAt(indexRec), &arrDaughIndex, std::array{0}, 1);
973-
if (arrDaughIndex.size() == 2) {
998+
if (arrDaughIndex.size() == NDaughtersResonant) {
974999
for (auto iProng = 0u; iProng < arrDaughIndex.size(); ++iProng) {
9751000
auto daughI = mcParticles.rawIteratorAt(arrDaughIndex[iProng]);
9761001
arrPDGDaugh[iProng] = std::abs(daughI.pdgCode());

PWGHF/TableProducer/candidateSelectorLc.cxx

Lines changed: 152 additions & 51 deletions
Original file line numberDiff line numberDiff line change
@@ -73,6 +73,9 @@ struct HfCandidateSelectorLc {
7373
// topological cuts
7474
Configurable<std::vector<double>> binsPt{"binsPt", std::vector<double>{hf_cuts_lc_to_p_k_pi::vecBinsPt}, "pT bin limits"};
7575
Configurable<LabeledArray<double>> cuts{"cuts", {hf_cuts_lc_to_p_k_pi::Cuts[0], hf_cuts_lc_to_p_k_pi::NBinsPt, hf_cuts_lc_to_p_k_pi::NCutVars, hf_cuts_lc_to_p_k_pi::labelsPt, hf_cuts_lc_to_p_k_pi::labelsCutVar}, "Lc candidate selection per pT bin"};
76+
Configurable<LabeledArray<double>> kfCuts{"kfCuts", {hf_cuts_lc_to_p_k_pi::CutsKf[0], hf_cuts_lc_to_p_k_pi::NBinsPt, hf_cuts_lc_to_p_k_pi::NCutKfVars, hf_cuts_lc_to_p_k_pi::labelsPt, hf_cuts_lc_to_p_k_pi::labelsCutKfVar}, "Lc candidate selection per pT bin with KF-associated variables"};
77+
Configurable<bool> applyNonKfCuts{"applyNonKfCuts", true, "Whether to apply non-KF cuts when running in KF mode. In DCAFitter mode this field is not effective"};
78+
Configurable<bool> applyKfCuts{"applyKfCuts", true, "Whether to apply KF cuts when running in KF mode. In DCAFitter mode this field is not effective"};
7679
// QA switch
7780
Configurable<bool> activateQA{"activateQA", false, "Flag to enable QA histogram"};
7881
// ML inference
@@ -175,7 +178,7 @@ struct HfCandidateSelectorLc {
175178
/// Conjugate-independent topological cuts
176179
/// \param candidate is candidate
177180
/// \return true if candidate passes all cuts
178-
template <typename T>
181+
template <aod::hf_cand::VertexerType reconstructionType, typename T>
179182
bool selectionTopol(const T& candidate)
180183
{
181184
auto candpT = candidate.pt();
@@ -190,33 +193,54 @@ struct HfCandidateSelectorLc {
190193
return false;
191194
}
192195

193-
// cosine of pointing angle
194-
if (candidate.cpa() <= cuts->get(pTBin, "cos pointing angle")) {
195-
return false;
196-
}
196+
if (reconstructionType == aod::hf_cand::VertexerType::DCAFitter || (reconstructionType == aod::hf_cand::VertexerType::KfParticle && applyNonKfCuts)) {
197+
// cosine of pointing angle
198+
if (candidate.cpa() <= cuts->get(pTBin, "cos pointing angle")) {
199+
return false;
200+
}
197201

198-
// candidate chi2PCA
199-
if (candidate.chi2PCA() > cuts->get(pTBin, "Chi2PCA")) {
200-
return false;
201-
}
202+
// candidate chi2PCA
203+
if (candidate.chi2PCA() > cuts->get(pTBin, "Chi2PCA")) {
204+
return false;
205+
}
202206

203-
if (candidate.decayLength() <= cuts->get(pTBin, "decay length")) {
204-
return false;
205-
}
207+
if (candidate.decayLength() <= cuts->get(pTBin, "decay length")) {
208+
return false;
209+
}
206210

207-
// candidate decay length XY
208-
if (candidate.decayLengthXY() <= cuts->get(pTBin, "decLengthXY")) {
209-
return false;
210-
}
211+
// candidate decay length XY
212+
if (candidate.decayLengthXY() <= cuts->get(pTBin, "decLengthXY")) {
213+
return false;
214+
}
211215

212-
// candidate normalized decay length XY
213-
if (candidate.decayLengthXYNormalised() < cuts->get(pTBin, "normDecLXY")) {
214-
return false;
216+
// candidate normalized decay length XY
217+
if (candidate.decayLengthXYNormalised() < cuts->get(pTBin, "normDecLXY")) {
218+
return false;
219+
}
220+
221+
// candidate impact parameter XY
222+
if (std::abs(candidate.impactParameterXY()) > cuts->get(pTBin, "impParXY")) {
223+
return false;
224+
}
215225
}
216226

217-
// candidate impact parameter XY
218-
if (std::abs(candidate.impactParameterXY()) > cuts->get(pTBin, "impParXY")) {
219-
return false;
227+
if constexpr (reconstructionType == aod::hf_cand::VertexerType::KfParticle) {
228+
if (applyKfCuts) {
229+
// candidate chi2geo of the triplet of prongs
230+
if (candidate.kfChi2Geo() > kfCuts->get(pTBin, "kfChi2Geo")) {
231+
return false;
232+
}
233+
234+
// candidate's decay point separation from the PV in terms of its error
235+
if (candidate.kfDecayLength() / candidate.kfDecayLengthError() < kfCuts->get(pTBin, "kfLdL")) {
236+
return false;
237+
}
238+
239+
// candidate's chi2 to the PV
240+
if (candidate.kfChi2Topo() > kfCuts->get(pTBin, "kfChi2Topo")) {
241+
return false;
242+
}
243+
}
220244
}
221245

222246
if (!isSelectedCandidateProngDca(candidate)) {
@@ -242,40 +266,117 @@ struct HfCandidateSelectorLc {
242266
return false;
243267
}
244268

245-
// cut on daughter pT
246-
if (trackProton.pt() < cuts->get(pTBin, "pT p") || trackKaon.pt() < cuts->get(pTBin, "pT K") || trackPion.pt() < cuts->get(pTBin, "pT Pi")) {
247-
return false;
248-
}
269+
if (reconstructionType == aod::hf_cand::VertexerType::DCAFitter || (reconstructionType == aod::hf_cand::VertexerType::KfParticle && applyNonKfCuts)) {
249270

250-
// cut on Lc->pKpi, piKp mass values
251-
/// cut on the Kpi pair invariant mass, to study Lc->pK*(->Kpi)
252-
float massLc, massKPi;
253-
if constexpr (reconstructionType == aod::hf_cand::VertexerType::DCAFitter) {
254-
if (trackProton.globalIndex() == candidate.prong0Id()) {
255-
massLc = hfHelper.invMassLcToPKPi(candidate);
256-
massKPi = hfHelper.invMassKPiPairLcToPKPi(candidate);
257-
} else {
258-
massLc = hfHelper.invMassLcToPiKP(candidate);
259-
massKPi = hfHelper.invMassKPiPairLcToPiKP(candidate);
271+
// cut on daughter pT
272+
if (trackProton.pt() < cuts->get(pTBin, "pT p") || trackKaon.pt() < cuts->get(pTBin, "pT K") || trackPion.pt() < cuts->get(pTBin, "pT Pi")) {
273+
return false;
260274
}
261-
} else if constexpr (reconstructionType == aod::hf_cand::VertexerType::KfParticle) {
262-
if (trackProton.globalIndex() == candidate.prong0Id()) {
263-
massLc = candidate.kfMassPKPi();
264-
massKPi = candidate.kfMassKPi();
265-
} else {
266-
massLc = candidate.kfMassPiKP();
267-
massKPi = candidate.kfMassPiK();
275+
276+
float massLc, massKPi;
277+
if constexpr (reconstructionType == aod::hf_cand::VertexerType::DCAFitter) {
278+
if (trackProton.globalIndex() == candidate.prong0Id()) {
279+
massLc = hfHelper.invMassLcToPKPi(candidate);
280+
massKPi = hfHelper.invMassKPiPairLcToPKPi(candidate);
281+
} else {
282+
massLc = hfHelper.invMassLcToPiKP(candidate);
283+
massKPi = hfHelper.invMassKPiPairLcToPiKP(candidate);
284+
}
285+
} else if constexpr (reconstructionType == aod::hf_cand::VertexerType::KfParticle) {
286+
if (trackProton.globalIndex() == candidate.prong0Id()) {
287+
massLc = candidate.kfMassPKPi();
288+
massKPi = candidate.kfMassKPi();
289+
} else {
290+
massLc = candidate.kfMassPiKP();
291+
massKPi = candidate.kfMassPiK();
292+
}
268293
}
269-
}
270294

271-
if (std::abs(massLc - o2::constants::physics::MassLambdaCPlus) > cuts->get(pTBin, "m")) {
272-
return false;
295+
// cut on Lc->pKpi, piKp mass values
296+
if (std::abs(massLc - o2::constants::physics::MassLambdaCPlus) > cuts->get(pTBin, "m")) {
297+
return false;
298+
}
299+
300+
/// cut on the Kpi pair invariant mass, to study Lc->pK*(->Kpi)
301+
const double cutMassKPi = cuts->get(pTBin, "mass (Kpi)");
302+
if (cutMassKPi > 0 && std::abs(massKPi - massK0Star892) > cutMassKPi) {
303+
return false;
304+
}
273305
}
274306

275-
/// cut on the Kpi pair invariant mass, to study Lc->pK*(->Kpi)
276-
const double cutMassKPi = cuts->get(pTBin, "mass (Kpi)");
277-
if (cutMassKPi > 0 && std::abs(massKPi - massK0Star892) > cutMassKPi) {
278-
return false;
307+
if constexpr (reconstructionType == aod::hf_cand::VertexerType::KfParticle) {
308+
if (applyKfCuts) {
309+
const float chi2PrimProng0 = candidate.kfChi2PrimProng0();
310+
const float chi2PrimProng1 = candidate.kfChi2PrimProng1();
311+
const float chi2PrimProng2 = candidate.kfChi2PrimProng2();
312+
313+
const float chi2GeoProng1Prong2 = candidate.kfChi2GeoProng1Prong2();
314+
const float chi2GeoProng0Prong2 = candidate.kfChi2GeoProng0Prong2();
315+
const float chi2GeoProng0Prong1 = candidate.kfChi2GeoProng0Prong1();
316+
317+
const float dcaProng1Prong2 = candidate.kfDcaProng1Prong2();
318+
const float dcaProng0Prong2 = candidate.kfDcaProng0Prong2();
319+
const float dcaProng0Prong1 = candidate.kfDcaProng0Prong1();
320+
321+
const double cutChi2PrimPr = kfCuts->get(pTBin, "kfChi2PrimPr");
322+
const double cutChi2PrimKa = kfCuts->get(pTBin, "kfChi2PrimKa");
323+
const double cutChi2PrimPi = kfCuts->get(pTBin, "kfChi2PrimPi");
324+
325+
const double cutChi2GeoKaPi = kfCuts->get(pTBin, "kfChi2GeoKaPi");
326+
const double cutChi2GeoPrPi = kfCuts->get(pTBin, "kfChi2GeoPrPi");
327+
const double cutChi2GeoPrKa = kfCuts->get(pTBin, "kfChi2GeoPrKa");
328+
329+
const double cutDcaKaPi = kfCuts->get(pTBin, "kfDcaKaPi");
330+
const double cutDcaPrPi = kfCuts->get(pTBin, "kfDcaPrPi");
331+
const double cutDcaPrKa = kfCuts->get(pTBin, "kfDcaPrKa");
332+
333+
// kaon's chi2 to the PV
334+
if (chi2PrimProng1 < cutChi2PrimKa) {
335+
return false;
336+
}
337+
338+
// chi2geo between proton and pion
339+
if (chi2GeoProng0Prong2 > cutChi2GeoPrPi) {
340+
return false;
341+
}
342+
343+
// dca between proton and pion
344+
if (dcaProng0Prong2 > cutDcaPrPi) {
345+
return false;
346+
}
347+
348+
const bool isPKPi = trackProton.globalIndex() == candidate.prong0Id();
349+
350+
// 0-th prong chi2 to the PV
351+
if (chi2PrimProng0 < (isPKPi ? cutChi2PrimPr : cutChi2PrimPi)) {
352+
return false;
353+
}
354+
355+
// 2-nd prong chi2 to the PV
356+
if (chi2PrimProng2 < (isPKPi ? cutChi2PrimPi : cutChi2PrimPr)) {
357+
return false;
358+
}
359+
360+
// chi2geo between 1-st and 2-nd prongs
361+
if (chi2GeoProng1Prong2 > (isPKPi ? cutChi2GeoKaPi : cutChi2GeoPrKa)) {
362+
return false;
363+
}
364+
365+
// chi2geo between 0-th and 1-st prongs
366+
if (chi2GeoProng0Prong1 > (isPKPi ? cutChi2GeoPrKa : cutChi2GeoKaPi)) {
367+
return false;
368+
}
369+
370+
// dca between 1-st and 2-nd prongs
371+
if (dcaProng1Prong2 > (isPKPi ? cutDcaKaPi : cutDcaPrKa)) {
372+
return false;
373+
}
374+
375+
// dca between 0-th and 1-st prongs
376+
if (dcaProng0Prong1 > (isPKPi ? cutDcaPrKa : cutDcaKaPi)) {
377+
return false;
378+
}
379+
}
279380
}
280381

281382
return true;
@@ -358,7 +459,7 @@ struct HfCandidateSelectorLc {
358459
}
359460

360461
// conjugate-independent topological selection
361-
if (!selectionTopol(candidate)) {
462+
if (!selectionTopol<reconstructionType>(candidate)) {
362463
hfSelLcCandidate(statusLcToPKPi, statusLcToPiKP);
363464
if (applyMl) {
364465
hfMlLcToPKPiCandidate(outputMlLcToPKPi, outputMlLcToPiKP);

0 commit comments

Comments
 (0)