Skip to content

Commit 43c76ca

Browse files
MarcellocostiMarcellocosti
committed
[WIP] D+ matching
1 parent 791e592 commit 43c76ca

File tree

5 files changed

+233
-32
lines changed

5 files changed

+233
-32
lines changed

‎Common/Core/RecoDecay.h‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -789,6 +789,8 @@ struct RecoDecay {
789789
// Check that the number of direct daughters is not larger than the number of expected final daughters.
790790
if constexpr (!acceptIncompleteReco && !checkProcess) {
791791
if (particleMother.daughtersIds().back() - particleMother.daughtersIds().front() + 1 > static_cast<int>(N)) {
792+
// LOG(info) << "MC Rec: Rejected: too many direct daughters: " << particleMother.daughtersIds().back() - particleMother.daughtersIds().front() + 1
793+
// << " (expected " << N << " final)";
792794
// Printf("MC Rec: Rejected: too many direct daughters: %d (expected %ld final)", particleMother.daughtersIds().back() - particleMother.daughtersIds().front() + 1, N);
793795
return -1;
794796
}
@@ -801,7 +803,9 @@ struct RecoDecay {
801803
// }
802804
// printf("\n");
803805
// Check whether the number of actual final daughters is equal to the number of expected final daughters (i.e. the number of provided prongs).
806+
// Printf("MC Rec: Number of final daughters: %ld (expected %ld)", arrAllDaughtersIndex.size(), N);
804807
if (!acceptIncompleteReco && arrAllDaughtersIndex.size() != N) {
808+
// LOG(info) << "MC Rec: Rejected: incorrect number of final daughters: " << arrAllDaughtersIndex.size() << " (expected " << N << ")";
805809
// Printf("MC Rec: Rejected: incorrect number of final daughters: %ld (expected %ld)", arrAllDaughtersIndex.size(), N);
806810
return -1;
807811
}

‎Common/TableProducer/centralityTable.cxx‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -240,6 +240,7 @@ struct CentralityTable {
240240
mRunNumber = 0;
241241
listCalib.setObject(new TList);
242242
if (!produceHistograms.value) {
243+
LOG(info) << "Debug histograms are disabled. No histograms will be produced.";
243244
return;
244245
}
245246

@@ -259,6 +260,7 @@ struct CentralityTable {
259260
histos.addClone("FT0C/", "sel8FT0C/");
260261
histos.addClone("FT0A/", "sel8FT0A/");
261262

263+
LOG(info) << "Debug histograms are enabled. Histograms will be produced.";
262264
histos.print();
263265
}
264266

@@ -569,6 +571,7 @@ struct CentralityTable {
569571
mftInfo.mCalibrationStored = false;
570572
if (callst != nullptr) {
571573
if (produceHistograms) {
574+
LOG(info) << "Adding calibration list to the list of calibration lists";
572575
listCalib->Add(callst->Clone(Form("%i", bc.runNumber())));
573576
}
574577
LOGF(info, "Getting new histograms with %d run number for %d run number", mRunNumber, bc.runNumber());

‎PWGHF/Core/CorrelatedBkgs.h‎

Lines changed: 194 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -13,25 +13,40 @@
1313
/// \brief Definitions of channels for correlated bkgs
1414
/// \author Marcello Di Costanzo <marcello.di.costanzo@cern.ch>, Polytechnic of Turin and INFN
1515

16+
#include <unordered_map>
17+
#include <array>
18+
#include <variant>
19+
1620
#ifndef PWGHF_CORE_CORRELATEDBKGTAGGER_H_
1721
#define PWGHF_CORE_CORRELATEDBKGTAGGER_H_
1822

23+
using namespace o2::constants::physics;
24+
1925
namespace o2::hf_corrbkg
2026
{
2127

28+
enum FinalStatesDMesons {
29+
KKPi = 0,
30+
KPiPi,
31+
KPiPiPi0,
32+
PiPiPi,
33+
PiPiPiPi0,
34+
NFinalStatesDMesons
35+
};
36+
2237
// D± → K± K∓ π±
2338
enum DecayChannelResoDplus {
24-
PhiPi = 0,
25-
K0starK,
39+
PhiPiDplus = 0,
40+
K0starKDplus,
2641
// K01430K,
2742
// RhoPi,
2843
// f2_1270Pi
2944
};
3045

3146
// Ds± → K± K∓ π±
3247
enum DecayChannelResoDs {
33-
PhiPi = DecayChannelResoDplus::PhiPi,
34-
K0starK,
48+
PhiPiDs = 0,
49+
K0starKDs,
3550
// PhiRho,
3651
// f2_1270Pi,
3752
// K0starPi,
@@ -60,15 +75,183 @@ namespace o2::hf_corrbkg
6075
// PhiP,
6176
// };
6277

63-
// Function to associate vectors with enum values
64-
std::vector<int> getFinalDaughters(DecayType finalState) {
65-
switch(finalState) {
66-
case DecayChannelCorrBkg::DstarToPiKPiBkg:
67-
return std::array{+kPiPlus, +kPiPlus, -kKPlus};
68-
default:
69-
return std::array{};
78+
79+
// using FinalStateVariant = std::variant<
80+
// std::array<int, 3>,
81+
// std::array<int, 4>
82+
// >;
83+
84+
// std::unordered_map<int, FinalStateVariant> finalStates =
85+
// {
86+
// {FinalStatesDMesons::KKPi, std::array<int,3>{+kKPlus, +kKPlus, +kPiPlus}},
87+
// {FinalStatesDMesons::KPiPiPi0, std::array<int,4>{+kKPlus, +kKPlus, +kPiPlus, +kPi0}}
88+
// };
89+
90+
std::unordered_map<FinalStatesDMesons, std::pair<std::array<int,3>, bool>> finalStates =
91+
{
92+
{FinalStatesDMesons::KKPi, {std::array<int,3>{-kKPlus, +kKPlus, +kPiPlus}, false}},
93+
{FinalStatesDMesons::KPiPi, {std::array<int,3>{+kPiPlus, -kKPlus, +kPiPlus}, false}},
94+
{FinalStatesDMesons::KPiPiPi0, {std::array<int,3>{+kPiPlus, -kKPlus, +kPiPlus}, true}},
95+
{FinalStatesDMesons::PiPiPi, {std::array<int,3>{+kPiMinus, +kPiPlus, +kPiPlus}, false}},
96+
{FinalStatesDMesons::PiPiPiPi0, {std::array<int,3>{+kPiMinus, +kPiPlus, +kPiPlus}, true}},
97+
};
98+
99+
100+
template <bool matchKinkedDecayTopology = false, bool matchInteractionsWithMaterial = false>
101+
int matchFinalStateDMeson(aod::McParticles const& mcParticles, auto arrayDaughters, int pdg, int8_t* flag, int8_t* sign, int depth, int8_t* nKinkedTracks, int8_t* nInteractionsWithMaterial) {
102+
int indexRec = -1;
103+
LOG(info) << "matchFinalStateDMeson: " << pdg;
104+
for (const std::pair<int, std::pair<std::array<int,3>, bool>>& finalState : finalStates) {
105+
if (finalState.second.second) { // 4 prong decay, partly reconstructed
106+
LOG(info) << "[4prong] Taking channel with final state: " << finalState.second.first[0] << " " << finalState.second.first[1] << " " << finalState.second.first[2] << " " << finalState.second.first[3] << ", has 4 daughters: " << finalState.second.second;
107+
LOG(info) << "[4prong] Trying 4-prong decay ... ";
108+
LOG(info) << "[4prong] Matching channel with final state: " << finalState.second.first[0] << " " << finalState.second.first[1] << " " << finalState.second.first[2];
109+
if constexpr (matchKinkedDecayTopology && matchInteractionsWithMaterial) {
110+
indexRec = RecoDecay::getMatchedMCRec<false, false, true, true, true>(mcParticles, arrayDaughters, pdg, finalState.second.first, true, sign, depth, nKinkedTracks, nInteractionsWithMaterial);
111+
} else if constexpr (matchKinkedDecayTopology && !matchInteractionsWithMaterial) {
112+
indexRec = RecoDecay::getMatchedMCRec<false, false, true, true, false>(mcParticles, arrayDaughters, pdg, finalState.second.first, true, sign, depth, nKinkedTracks);
113+
} else if constexpr (!matchKinkedDecayTopology && matchInteractionsWithMaterial) {
114+
indexRec = RecoDecay::getMatchedMCRec<false, false, true, false, true>(mcParticles, arrayDaughters, pdg, finalState.second.first, true, sign, depth, nullptr, nInteractionsWithMaterial);
115+
} else {
116+
indexRec = RecoDecay::getMatchedMCRec<false, false, true, false, false>(mcParticles, arrayDaughters, pdg, finalState.second.first, true, &sign, depth);
117+
}
118+
} else {
119+
LOG(info) << "[3prong] Trying 3-prong decay ... ";
120+
LOG(info) << "[3prong] Taking channel with final state: " << finalState.second.first[0] << " " << finalState.second.first[1] << " " << finalState.second.first[2] << ", has 4 daughters: " << finalState.second.second;
121+
if constexpr (matchKinkedDecayTopology && matchInteractionsWithMaterial) {
122+
indexRec = RecoDecay::getMatchedMCRec<false, false, false, true, true>(mcParticles, arrayDaughters, pdg, finalState.second.first, true, sign, depth, nKinkedTracks, nInteractionsWithMaterial);
123+
} else if constexpr (matchKinkedDecayTopology && !matchInteractionsWithMaterial) {
124+
indexRec = RecoDecay::getMatchedMCRec<false, false, false, true, false>(mcParticles, arrayDaughters, pdg, finalState.second.first, true, sign, depth, nKinkedTracks);
125+
} else if constexpr (!matchKinkedDecayTopology && matchInteractionsWithMaterial) {
126+
indexRec = RecoDecay::getMatchedMCRec<false, false, false, false, true>(mcParticles, arrayDaughters, pdg, finalState.second.first, true, sign, depth, nullptr, nInteractionsWithMaterial);
127+
} else {
128+
indexRec = RecoDecay::getMatchedMCRec<false, false, false, false, false>(mcParticles, arrayDaughters, pdg, finalState.second.first, true, &sign, depth);
129+
}
130+
}
131+
if (indexRec > -1) {
132+
LOG(info) << "Found channel with final state: " << finalState.second.first[0] << " " << finalState.second.first[1] << " " << finalState.second.first[2] << ", has 4 prongs: " << finalState.second.second;
133+
*flag = (*sign) * (1 << finalState.first);
134+
break;
135+
}
70136
}
137+
return indexRec;
71138
}
139+
140+
141+
142+
143+
144+
// Function to associate vectors with enum values
145+
// std::vector<int> getFinalDaughters(DecayType finalState) {
146+
// switch(finalState) {
147+
// case DecayChannelCorrBkg::DstarToPiKPiBkg:
148+
// return std::array{+kPiPlus, +kPiPlus, -kKPlus};
149+
// default:
150+
// return std::array{};
151+
// }
152+
// }
153+
154+
// std::vector<Pdg> getFinalStateDMeson(int index) {
155+
// switch (index) {
156+
// case KKPi:
157+
// return std::vector<Pdg>{+kKPlus, -kKPlus, +kPiPlus};
158+
// case KPiPi:
159+
// return std::vector<Pdg>{+kKPlus, +kPiPlus, +kPiPlus};
160+
// case KPiPiPi0:
161+
// return std::vector<Pdg>{+kKPlus, +kPiPlus, +kPiPlus, +kPi0};
162+
// case PiPiPi:
163+
// return std::vector<Pdg>{+kPiPlus, +kPiPlus, +kPiPlus};
164+
// case PiPiPiPi0:
165+
// return std::vector<Pdg>{+kPiPlus, +kPiPlus, +kPiPlus, +kPi0};
166+
// default:
167+
// throw std::out_of_range("Invalid index for final state D meson");
168+
// }
169+
// }
170+
171+
172+
// template <int N>
173+
// std::array<int, N> getFinalStateDMeson(int index) {
174+
// if constexpr (N == 3) {
175+
// switch (index) {
176+
// case FinalStatesDMesons::KKPi:
177+
// return std::array<int, N>{+kKPlus, -kKPlus, +kPiPlus};
178+
// case FinalStatesDMesons::KPiPi:
179+
// return std::array<int, N>{+kKPlus, +kPiPlus, +kPiPlus};
180+
// case FinalStatesDMesons::PiPiPi:
181+
// return std::array<int, N>{+kPiPlus, +kPiPlus, +kPiPlus};
182+
// default:
183+
// throw std::out_of_range("Invalid index for final state D meson");
184+
// }
185+
// } else if constexpr (N == 4) {
186+
// switch (index) {
187+
// case FinalStatesDMesons::KPiPiPi0:
188+
// return std::array<int, N>{+kKPlus, +kPiPlus, +kPiPlus, +kPi0};
189+
// case FinalStatesDMesons::PiPiPiPi0:
190+
// return std::array<int, N>{+kPiPlus, +kPiPlus, +kPiPlus, +kPi0};
191+
// default:
192+
// throw std::out_of_range("Invalid index for final state D meson");
193+
// }
194+
// } else {
195+
// throw std::invalid_argument("Invalid size for final state D meson");
196+
// }
197+
// }
198+
199+
// int getFinalStateParts(int index) {
200+
// switch (index) {
201+
// case FinalStatesDMesons::KKPi:
202+
// return 3;
203+
// case FinalStatesDMesons::KPiPi:
204+
// return 3;
205+
// case FinalStatesDMesons::KPiPiPi0:
206+
// return 4;
207+
// case FinalStatesDMesons::PiPiPi:
208+
// return 3;
209+
// case FinalStatesDMesons::PiPiPiPi0:
210+
// return 4;
211+
// default:
212+
// throw std::out_of_range("Invalid index for final state D meson");
213+
// }
214+
// }
215+
216+
// template <bool matchKinkedDecayTopology = false, bool matchInteractionsWithMaterial = false>
217+
// int matchFinalStateDMeson(int iChannel, aod::McParticles const& mcParticles, auto arrayDaughters, int pdg, int8_t* sign, int depth, int8_t* nKinkedTracks, int8_t* nInteractionsWithMaterial) {
218+
// int indexRec = -1;
219+
// int nFinalStateParts = getFinalStateParts(iChannel);
220+
// if (nFinalStateParts == 3) {
221+
// auto finalState = getFinalStateDMeson<3>(iChannel);
222+
// if constexpr (matchKinkedDecayTopology && matchInteractionsWithMaterial) {
223+
// indexRec = RecoDecay::getMatchedMCRec<false, false, false, true, true>(mcParticles, arrayDaughters, pdg, finalState, true, sign, depth, nKinkedTracks, nInteractionsWithMaterial);
224+
// } else if constexpr (matchKinkedDecayTopology && !matchInteractionsWithMaterial) {
225+
// indexRec = RecoDecay::getMatchedMCRec<false, false, false, true, false>(mcParticles, arrayDaughters, pdg, finalState, true, sign, depth, nKinkedTracks);
226+
// } else if constexpr (!matchKinkedDecayTopology && matchInteractionsWithMaterial) {
227+
// indexRec = RecoDecay::getMatchedMCRec<false, false, false, false, true>(mcParticles, arrayDaughters, pdg, finalState, true, sign, depth, nullptr, nInteractionsWithMaterial);
228+
// } else {
229+
// indexRec = RecoDecay::getMatchedMCRec(mcParticles, arrayDaughters, pdg, finalState, true, &sign, depth);
230+
// }
231+
// if (indexRec > -1) {
232+
// LOG(info) << "Matched final state: " << iChannel;
233+
// return indexRec;
234+
// }
235+
// } else if (nFinalStateParts == 4) {
236+
// auto finalState = getFinalStateDMeson<4>(iChannel);
237+
// if constexpr (matchKinkedDecayTopology && matchInteractionsWithMaterial) {
238+
// indexRec = RecoDecay::getMatchedMCRec<false, false, false, true, true>(mcParticles, arrayDaughters, pdg, finalState, true, sign, depth, nKinkedTracks, nInteractionsWithMaterial);
239+
// } else if constexpr (matchKinkedDecayTopology && !matchInteractionsWithMaterial) {
240+
// indexRec = RecoDecay::getMatchedMCRec<false, false, false, true, false>(mcParticles, arrayDaughters, pdg, finalState, true, sign, depth, nKinkedTracks);
241+
// } else if constexpr (!matchKinkedDecayTopology && matchInteractionsWithMaterial) {
242+
// indexRec = RecoDecay::getMatchedMCRec<false, false, false, false, true>(mcParticles, arrayDaughters, pdg, finalState, true, sign, depth, nullptr, nInteractionsWithMaterial);
243+
// } else {
244+
// indexRec = RecoDecay::getMatchedMCRec(mcParticles, arrayDaughters, pdg, finalState, true, &sign, depth);
245+
// }
246+
// if (indexRec > -1) {
247+
// LOG(info) << "Matched final state: " << iChannel;
248+
// return indexRec;
249+
// }
250+
// } else {
251+
// throw std::invalid_argument("Invalid size for final state D meson");
252+
// }
253+
// }
254+
72255
} // namespace o2::hf_corrbkg
73256

74257
#endif // PWGHF_CORE_CORRELATEDBKGTAGGER_H_

‎PWGHF/TableProducer/candidateCreator3Prong.cxx‎

Lines changed: 18 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -47,6 +47,7 @@
4747

4848
#include "PWGHF/Core/CentralityEstimation.h"
4949
#include "PWGHF/Core/DecayChannels.h"
50+
#include "PWGHF/Core/CorrelatedBkgs.h"
5051
#include "PWGHF/DataModel/CandidateReconstructionTables.h"
5152
#include "PWGHF/Utils/utilsBfieldCCDB.h"
5253
#include "PWGHF/Utils/utilsEvSelHf.h"
@@ -61,6 +62,7 @@ using namespace o2::aod::hf_cand_3prong;
6162
using namespace o2::hf_decay::hf_cand_3prong;
6263
using namespace o2::hf_centrality;
6364
using namespace o2::hf_occupancy;
65+
using namespace o2::hf_corrbkg;
6466
using namespace o2::constants::physics;
6567
using namespace o2::framework;
6668
using namespace o2::framework::expressions;
@@ -941,7 +943,10 @@ struct HfCandidateCreator3ProngExpressions {
941943
}
942944
}
943945

946+
944947
// D± → π± K∓ π±
948+
LOG(info) << "--------------------------------------------";
949+
LOG(info) << "D+ -> Pi K Pi";
945950
if (flag == 0) {
946951
if (matchKinkedDecayTopology && matchInteractionsWithMaterial) {
947952
indexRec = RecoDecay::getMatchedMCRec<false, false, false, true, true>(mcParticles, arrayDaughters, Pdg::kDPlus, std::array{+kPiPlus, -kKPlus, +kPiPlus}, true, &sign, 2, &nKinkedTracks, &nInteractionsWithMaterial);
@@ -956,6 +961,17 @@ struct HfCandidateCreator3ProngExpressions {
956961
flag = sign * DecayChannelMain::DplusToPiKPi;
957962
}
958963
}
964+
LOG(info) << "D+ matching ended";
965+
966+
// if (matchKinkedDecayTopology && matchInteractionsWithMaterial) {
967+
// indexRec = matchFinalStateDMeson<true, true>(mcParticles, arrayDaughters, Pdg::kDPlus, &sign, 2, &nKinkedTracks, &nInteractionsWithMaterial);
968+
// } else if (matchKinkedDecayTopology && !matchInteractionsWithMaterial) {
969+
// indexRec = matchFinalStateDMeson<true, false>(mcParticles, arrayDaughters, Pdg::kDPlus, &sign, 2, &nKinkedTracks, &nInteractionsWithMaterial);
970+
// } else if (!matchKinkedDecayTopology && matchInteractionsWithMaterial) {
971+
// indexRec = matchFinalStateDMeson<true, false>(mcParticles, arrayDaughters, Pdg::kDPlus, &sign, 2, &nKinkedTracks, &nInteractionsWithMaterial);
972+
// } else {
973+
// indexRec = matchFinalStateDMeson<false, false>(mcParticles, arrayDaughters, Pdg::kDPlus, &sign, 2, &nKinkedTracks, &nInteractionsWithMaterial);
974+
// }
959975

960976
// Ds± → K± K∓ π± and D± → K± K∓ π±
961977
if (flag == 0) {
@@ -1028,27 +1044,6 @@ struct HfCandidateCreator3ProngExpressions {
10281044
if (indexRec > -1) {
10291045
flag = sign * DecayChannelMain::LcToPKPi;
10301046

1031-
// Flagging the different Λc± → p± K∓ π± decay channels
1032-
if (arrayDaughters[0].has_mcParticle()) {
1033-
swapping = int8_t(std::abs(arrayDaughters[0].mcParticle().pdgCode()) == kPiPlus);
1034-
}
1035-
RecoDecay::getDaughters(mcParticles.rawIteratorAt(indexRec), &arrDaughIndex, std::array{0}, 1);
1036-
if (arrDaughIndex.size() == NDaughtersResonant) {
1037-
for (auto iProng = 0u; iProng < arrDaughIndex.size(); ++iProng) {
1038-
auto daughI = mcParticles.rawIteratorAt(arrDaughIndex[iProng]);
1039-
arrPDGDaugh[iProng] = std::abs(daughI.pdgCode());
1040-
}
1041-
if ((arrPDGDaugh[0] == arrPDGResonant1[0] && arrPDGDaugh[1] == arrPDGResonant1[1]) || (arrPDGDaugh[0] == arrPDGResonant1[1] && arrPDGDaugh[1] == arrPDGResonant1[0])) {
1042-
channel = 1;
1043-
} else if ((arrPDGDaugh[0] == arrPDGResonant2[0] && arrPDGDaugh[1] == arrPDGResonant2[1]) || (arrPDGDaugh[0] == arrPDGResonant2[1] && arrPDGDaugh[1] == arrPDGResonant2[0])) {
1044-
channel = 2;
1045-
} else if ((arrPDGDaugh[0] == arrPDGResonant3[0] && arrPDGDaugh[1] == arrPDGResonant3[1]) || (arrPDGDaugh[0] == arrPDGResonant3[1] && arrPDGDaugh[1] == arrPDGResonant3[0])) {
1046-
channel = 3;
1047-
}
1048-
}
1049-
}
1050-
}
1051-
10521047
// Ξc± → p± K∓ π±
10531048
if (flag == 0) {
10541049
if (matchKinkedDecayTopology && matchInteractionsWithMaterial) {
@@ -1070,6 +1065,7 @@ struct HfCandidateCreator3ProngExpressions {
10701065

10711066
// Check whether the particle is non-prompt (from a b quark).
10721067
if (flag != 0) {
1068+
LOG(info) << "Setting origin";
10731069
auto particle = mcParticles.rawIteratorAt(indexRec);
10741070
origin = RecoDecay::getCharmHadronOrigin(mcParticles, particle, false, &idxBhadMothers);
10751071
}
@@ -1079,6 +1075,7 @@ struct HfCandidateCreator3ProngExpressions {
10791075
} else {
10801076
rowMcMatchRec(flag, origin, swapping, channel, -1.f, 0, nKinkedTracks, nInteractionsWithMaterial);
10811077
}
1078+
LOG(info) << "--------------------------------------------";
10821079
}
10831080

10841081
for (const auto& mcCollision : mcCollisions) {

0 commit comments

Comments
 (0)