Skip to content

Commit c4f948b

Browse files
MarcellocostiMarcellocosti
committed
Fix Dstar matching and add other channels
1 parent f72db07 commit c4f948b

File tree

4 files changed

+191
-79
lines changed

4 files changed

+191
-79
lines changed

‎Common/Core/RecoDecay.h‎

Lines changed: 40 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -24,6 +24,7 @@
2424
#include <cstdio>
2525
#include <utility> // std::move
2626
#include <vector> // std::vector
27+
#include <iostream>
2728

2829
// ROOT includes
2930
#include <TMCProcess.h> // for VMC Particle Production Process
@@ -904,7 +905,11 @@ struct RecoDecay {
904905
}
905906
// Check the PDG code of the particle.
906907
auto pdgCandidate = candidate.pdgCode();
908+
bool print = false; // Print debug messages only if print is true.
907909
// Printf("MC Gen: Candidate PDG: %d", pdgCandidate);
910+
if (std::abs(pdgCandidate) == 413) {
911+
print = true; // Print debug messages for D* candidates.
912+
}
908913
if (pdgCandidate == pdgParticle) { // exact PDG match
909914
sgn = 1;
910915
} else if (acceptAntiParticles && pdgCandidate == -pdgParticle) { // antiparticle PDG match
@@ -915,29 +920,56 @@ struct RecoDecay {
915920
}
916921
// Check the PDG codes of the decay products.
917922
if (N > 0) {
923+
if (print) {
924+
std::cout << "MC Gen: Checking decay products of " << N << " daughters" << std::endl;
925+
// Printf("MC Gen: Checking decay products of %ld daughters", N);
926+
}
918927
// Printf("MC Gen: Checking %d daughters", N);
919928
std::vector<int> arrAllDaughtersIndex; // vector of indices of all daughters
920929
// Check the daughter indices.
921930
if (!candidate.has_daughters()) {
931+
if (print) {
932+
std::cout << "MC Gen: Rejected: bad daughter index range: " << candidate.daughtersIds().front() << "-" << candidate.daughtersIds().back() << std::endl;
933+
// Printf("MC Gen: Rejected: bad daughter index range: %d-%d", candidate.daughtersIds().front(), candidate.daughtersIds().back());
934+
}
922935
// Printf("MC Gen: Rejected: bad daughter index range: %d-%d", candidate.daughtersIds().front(), candidate.daughtersIds().back());
923936
return false;
924937
}
925938
// Check that the number of direct daughters is not larger than the number of expected final daughters.
926939
if constexpr (!checkProcess) {
927940
if (candidate.daughtersIds().back() - candidate.daughtersIds().front() + 1 > static_cast<int>(N)) {
941+
if (print) {
942+
std::cout << "MC Gen: Rejected: too many direct daughters: " << candidate.daughtersIds().back() - candidate.daughtersIds().front() + 1
943+
<< " (expected " << N << " final)" << std::endl;
944+
// Printf("MC Gen: Rejected: too many direct daughters: %d (expected %ld final)", candidate.daughtersIds().back() - candidate.daughtersIds().front() + 1, N);
945+
}
928946
// Printf("MC Gen: Rejected: too many direct daughters: %d (expected %ld final)", candidate.daughtersIds().back() - candidate.daughtersIds().front() + 1, N);
929947
return false;
930948
}
931949
}
932950
// Get the list of actual final daughters.
933951
getDaughters<checkProcess>(candidate, &arrAllDaughtersIndex, arrPdgDaughters, depthMax);
952+
if (print) {
953+
std::cout << "MC Gen: Mother " << candidate.globalIndex() << " has " << arrAllDaughtersIndex.size() << " final daughters:";
954+
// Printf("MC Gen: Mother %ld has %ld final daughters:", candidate.globalIndex(), arrAllDaughtersIndex.size());
955+
for (auto i : arrAllDaughtersIndex) {
956+
std::cout << " " << i;
957+
// Printf(" %d", i);
958+
}
959+
std::cout << std::endl;
960+
// Printf("\n");
961+
}
934962
// printf("MC Gen: Mother %ld has %ld final daughters:", candidate.globalIndex(), arrAllDaughtersIndex.size());
935963
// for (auto i : arrAllDaughtersIndex) {
936964
// printf(" %d", i);
937965
// }
938966
// printf("\n");
939967
// Check whether the number of final daughters is equal to the required number.
940968
if (arrAllDaughtersIndex.size() != N) {
969+
if (print) {
970+
std::cout << "MC Gen: Rejected: incorrect number of final daughters: " << arrAllDaughtersIndex.size() << " (expected " << N << ")" << std::endl;
971+
// Printf("MC Gen: Rejected: incorrect number of final daughters: %ld (expected %ld)", arrAllDaughtersIndex.size(), N);
972+
}
941973
// Printf("MC Gen: Rejected: incorrect number of final daughters: %ld (expected %ld)", arrAllDaughtersIndex.size(), N);
942974
return false;
943975
}
@@ -955,6 +987,10 @@ struct RecoDecay {
955987
for (auto indexDaughterI : arrAllDaughtersIndex) { // o2-linter: disable=const-ref-in-for-loop (int elements)
956988
auto candidateDaughterI = particlesMC.rawIteratorAt(indexDaughterI - particlesMC.offset()); // ith daughter particle
957989
auto pdgCandidateDaughterI = candidateDaughterI.pdgCode(); // PDG code of the ith daughter
990+
if (print) {
991+
std::cout << "MC Gen: Daughter " << indexDaughterI << " PDG: " << pdgCandidateDaughterI << std::endl;
992+
// Printf("MC Gen: Daughter %d PDG: %d", indexDaughterI, pdgCandidateDaughterI);
993+
}
958994
// Printf("MC Gen: Daughter %d PDG: %d", indexDaughterI, pdgCandidateDaughterI);
959995
bool isPdgFound = false; // Is the PDG code of this daughter among the remaining expected PDG codes?
960996
for (std::size_t iProngCp = 0; iProngCp < N; ++iProngCp) {
@@ -973,6 +1009,10 @@ struct RecoDecay {
9731009
*listIndexDaughters = arrAllDaughtersIndex;
9741010
}
9751011
}
1012+
if (print) {
1013+
std::cout << "MC Gen: Accepted: m: " << candidate.globalIndex() << std::endl;
1014+
// Printf("MC Gen: Accepted: m: %ld", candidate.globalIndex());
1015+
}
9761016
// Printf("MC Gen: Accepted: m: %d", candidate.globalIndex());
9771017
if (sign) {
9781018
*sign = sgn;

‎PWGHF/Core/CorrelatedBkgs.h‎

Lines changed: 90 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -31,43 +31,53 @@ namespace o2::hf_corrbkg
3131
enum FinalStates2Prongs {
3232
KPi = 1,
3333
KK,
34-
KPiPi0,
34+
KMinusPiPi0,
3535
PiPi,
3636
PiPiPi0,
3737
NFinalStates2P
3838
};
3939

4040
std::unordered_map<FinalStates2Prongs, std::vector<int> > finalStates2Prongs =
4141
{
42-
{FinalStates2Prongs::KPi, std::vector<int>{+kKMinus, +kPiPlus}},
43-
{FinalStates2Prongs::KK, std::vector<int>{+kKMinus, +kKPlus}},
44-
{FinalStates2Prongs::KPiPi0, std::vector<int>{+kKMinus, +kPiPlus, +kPi0}},
45-
{FinalStates2Prongs::PiPi, std::vector<int>{+kPiMinus, +kPiPlus}},
46-
{FinalStates2Prongs::PiPiPi0, std::vector<int>{+kPiMinus, +kPiPlus, +kPi0}}
42+
{FinalStates2Prongs::KPi, std::vector<int>{+kKMinus, +kPiPlus}},
43+
{FinalStates2Prongs::KK, std::vector<int>{+kKMinus, +kKPlus}},
44+
{FinalStates2Prongs::KMinusPiPi0, std::vector<int>{+kKMinus, +kPiPlus, +kPi0}},
45+
{FinalStates2Prongs::PiPi, std::vector<int>{+kPiMinus, +kPiPlus}},
46+
{FinalStates2Prongs::PiPiPi0, std::vector<int>{+kPiMinus, +kPiPlus, +kPi0}}
4747
};
4848

4949
enum FinalStates3Prongs {
5050
KKPi = 1,
5151
KKPiPi0,
52-
KPiPi,
53-
KPiPiPi0,
52+
KMinusPiPi,
53+
KMinusPiPiPi0,
54+
KPlusPiPi,
55+
KPlusPiPiPi0,
5456
PiPiPi,
5557
PiPiPiPi0,
56-
ProtonPiPi,
5758
ProtonKPi,
59+
ProtonKPiPi0,
60+
ProtonPiPi,
61+
ProtonKK,
62+
SigmaPiPi,
5863
NFinalStates3P
5964
};
6065

6166
std::unordered_map<FinalStates3Prongs, std::vector<int> > finalStates3Prongs =
6267
{
63-
{FinalStates3Prongs::KKPi, std::vector<int>{+kKMinus, +kKPlus, +kPiPlus}},
64-
{FinalStates3Prongs::KKPiPi0, std::vector<int>{+kKMinus, +kKPlus, +kPiPlus, +kPi0}},
65-
{FinalStates3Prongs::KPiPi, std::vector<int>{+kKMinus, +kPiPlus, +kPiPlus}},
66-
{FinalStates3Prongs::KPiPiPi0, std::vector<int>{+kKMinus, +kPiPlus, +kPiPlus, +kPi0}},
67-
{FinalStates3Prongs::PiPiPi, std::vector<int>{+kPiMinus, +kPiPlus, +kPiPlus}},
68-
{FinalStates3Prongs::PiPiPiPi0, std::vector<int>{+kPiMinus, +kPiPlus, +kPiPlus, +kPi0}},
69-
{FinalStates3Prongs::ProtonPiPi, std::vector<int>{+kProton, +kPiMinus, +kPiPlus}},
70-
{FinalStates3Prongs::ProtonKPi, std::vector<int>{+kProton, +kKMinus, +kPiPlus}}
68+
{FinalStates3Prongs::KKPi, std::vector<int>{+kKMinus, +kKPlus, +kPiPlus}},
69+
{FinalStates3Prongs::KKPiPi0, std::vector<int>{+kKMinus, +kKPlus, +kPiPlus, +kPi0}},
70+
{FinalStates3Prongs::KMinusPiPi, std::vector<int>{+kKMinus, +kPiPlus, +kPiPlus}},
71+
{FinalStates3Prongs::KMinusPiPiPi0, std::vector<int>{+kKMinus, +kPiPlus, +kPiPlus, +kPi0}},
72+
{FinalStates3Prongs::KPlusPiPi, std::vector<int>{+kPlus, +kPiPlus, +kPiMinus}},
73+
{FinalStates3Prongs::KPlusPiPiPi0, std::vector<int>{+kPlus, +kPiPlus, +kPiMinus, +kPi0}},
74+
{FinalStates3Prongs::PiPiPi, std::vector<int>{+kPiMinus, +kPiPlus, +kPiPlus}},
75+
{FinalStates3Prongs::PiPiPiPi0, std::vector<int>{+kPiMinus, +kPiPlus, +kPiPlus, +kPi0}},
76+
{FinalStates3Prongs::ProtonKPi, std::vector<int>{+kProton, +kKMinus, +kPiPlus}},
77+
{FinalStates3Prongs::ProtonKPiPi0, std::vector<int>{+kProton, +kKMinus, +kPiPlus, +kPi0}},
78+
{FinalStates3Prongs::ProtonPiPi, std::vector<int>{+kProton, +kPiMinus, +kPiPlus}},
79+
{FinalStates3Prongs::ProtonKK, std::vector<int>{+kProton, +kKMinus, +kKPlus}},
80+
{FinalStates3Prongs::SigmaPiPi, std::vector<int>{+kSigmaPlus, +kPiMinus, +kPiPlus}},
7181
};
7282

7383
// Dstar → K± K∓ π±
@@ -106,6 +116,15 @@ namespace o2::hf_corrbkg
106116
{DecayChannelResoDplus::RhoPi, std::array<int, 2>{+113, +kPiPlus}},
107117
{DecayChannelResoDplus::f2_1270Pi, std::array<int, 2>{+225, +kPiPlus}},
108118
};
119+
120+
std::unordered_map<FinalStates3Prongs, std::vector<int> > finalStatesDPlus =
121+
{
122+
{FinalStates3Prongs::KKPi, std::vector<int>{+kKMinus, +kKPlus, +kPiPlus}},
123+
{FinalStates3Prongs::KMinusPiPi, std::vector<int>{+kKMinus, +kPiPlus, +kPiPlus}},
124+
{FinalStates3Prongs::KMinusPiPiPi0, std::vector<int>{+kKMinus, +kPiPlus, +kPiPlus, +kPi0}},
125+
{FinalStates3Prongs::PiPiPi, std::vector<int>{+kPiMinus, +kPiPlus, +kPiPlus}},
126+
{FinalStates3Prongs::PiPiPiPi0, std::vector<int>{+kPiMinus, +kPiPlus, +kPiPlus, +kPi0}},
127+
};
109128
}
110129

111130
// Ds± → K± K∓ π±
@@ -134,6 +153,15 @@ namespace o2::hf_corrbkg
134153
{DecayChannelResoDs::f2_1270Pi, std::array<int, 2>{225, +kPiPlus}},
135154
{DecayChannelResoDs::f2_1370K, std::array<int, 2>{10221, +kKPlus}},
136155
};
156+
157+
std::unordered_map<FinalStates3Prongs, std::vector<int> > finalStatesDs =
158+
{
159+
{FinalStates3Prongs::KKPi, std::vector<int>{+kKMinus, +kKPlus, +kPiPlus}},
160+
{FinalStates3Prongs::KKPiPi0, std::vector<int>{+kKMinus, +kKPlus, +kPiPlus, +kPi0}},
161+
{FinalStates3Prongs::KPlusPiPi, std::vector<int>{+kKPlus, +kPiPlus, +kPiMinus}},
162+
{FinalStates3Prongs::PiPiPi, std::vector<int>{+kPiMinus, +kPiPlus, +kPiPlus}},
163+
{FinalStates3Prongs::PiPiPiPi0, std::vector<int>{+kPiMinus, +kPiPlus, +kPiPlus, +kPi0}},
164+
};
137165
}
138166

139167
// Dstar → K± K∓ π±
@@ -152,6 +180,15 @@ namespace o2::hf_corrbkg
152180
{DecayChannelResoDStarD0::K0starPi0, std::array<int, 2>{-kK0Star892, +kPi0}},
153181
{DecayChannelResoDStarD0::K0starPiPlus, std::array<int, 2>{-kKPlusStar892, +kPiPlus}},
154182
};
183+
184+
std::unordered_map<FinalStates3Prongs, std::vector<int> > finalStatesDStar =
185+
{
186+
{FinalStates3Prongs::KKPi, std::vector<int>{+kKMinus, +kKPlus, +kPiPlus}},
187+
{FinalStates3Prongs::KMinusPiPi, std::vector<int>{+kKMinus, +kPiPlus, +kPiPlus}},
188+
{FinalStates3Prongs::KMinusPiPiPi0, std::vector<int>{+kKMinus, +kPiPlus, +kPiPlus, +kPi0}},
189+
{FinalStates3Prongs::PiPiPi, std::vector<int>{+kPiMinus, +kPiPlus, +kPiPlus}},
190+
{FinalStates3Prongs::PiPiPiPi0, std::vector<int>{+kPiMinus, +kPiPlus, +kPiPlus, +kPi0}},
191+
};
155192
}
156193

157194
// Lc → p K∓ π±
@@ -170,6 +207,14 @@ namespace o2::hf_corrbkg
170207
{DecayChannelResoLambdaC::DeltaK, std::array<int, 2>{+2224, +kKMinus}},
171208
{DecayChannelResoLambdaC::Lambda1520Pi, std::array<int, 2>{+102134, +kPiPlus}},
172209
};
210+
211+
std::unordered_map<FinalStates3Prongs, std::vector<int> > finalStatesLc =
212+
{
213+
{FinalStates3Prongs::ProtonKPi, std::vector<int>{+kProton, +kKMinus, +kPiPlus}},
214+
{FinalStates3Prongs::ProtonKPiPi0, std::vector<int>{+kProton, +kKMinus, +kPiPlus, +kPi0}},
215+
{FinalStates3Prongs::ProtonPiPi, std::vector<int>{+kProton, +kPiMinus, +kPiPlus}},
216+
{FinalStates3Prongs::ProtonKK, std::vector<int>{+kProton, +kKMinus, +kKPlus}}
217+
};
173218
}
174219

175220
// Xic → p K∓ π±
@@ -186,9 +231,36 @@ namespace o2::hf_corrbkg
186231
{DecayChannelResoXiC::ProtonPhi, std::array<int, 2>{+kProton, +kPhi}},
187232
{DecayChannelResoXiC::SigmaPiPi, std::array<int, 2>{+kSigmaPlus, +kPiPlus}},
188233
};
234+
235+
std::unordered_map<FinalStates3Prongs, std::vector<int> > finalStatesXic =
236+
{
237+
{FinalStates3Prongs::ProtonKPi, std::vector<int>{+kProton, +kKMinus, +kPiPlus}},
238+
{FinalStates3Prongs::ProtonKK, std::vector<int>{+kProton, +kKMinus, +kKPlus}},
239+
{FinalStates3Prongs::SigmaPiPi, std::vector<int>{+kSigmaPlus, +kPiMinus, +kPiPlus}},
240+
};
241+
}
242+
243+
std::unordered_map<FinalStates3Prongs, std::vector<int> > getParticleFinalStates3Prongs(int pdgMother) {
244+
switch (pdgMother) {
245+
case Pdg::kDPlus:
246+
return DPlus::finalStatesDPlus;
247+
case Pdg::kDS:
248+
return DS::finalStatesDs;
249+
case Pdg::kDStar:
250+
return DStar::finalStatesDStar;
251+
case Pdg::kLambdaCPlus:
252+
return LambdaC::finalStatesLc;
253+
case Pdg::kXiCPlus:
254+
return XiC::finalStatesXic;
255+
default:
256+
LOG(error) << "Unknown PDG code for 3-prong final states: " << pdgMother;
257+
return {};
258+
}
189259
}
190260

191261
bool checkResonantDecay(std::vector<int> arrDaughIndex, std::array<int, 2> arrPDGResonant) {
262+
// LOG(info) << "Entered checkResonantDecay with daughters: " << arrDaughIndex[0] << ", " << arrDaughIndex[1] << " and resonant PDG codes: " << arrPDGResonant[0] << ", " << arrPDGResonant[1];
263+
// LOG(info) << "arrDaughIndex.size(): " << arrDaughIndex.size() << ", arrPDGResonant.size(): " << arrPDGResonant.size();
192264
std::array<int, 2> arrPDGResonantAbs = {std::abs(arrPDGResonant[0]), std::abs(arrPDGResonant[1])};
193265
LOG(info) << "Testing: " << arrDaughIndex[0] << ", " << arrDaughIndex[1] << " matching PDG codes: " << arrPDGResonant[0] << ", " << arrPDGResonant[1];
194266
for (int i = 0; i < 2; i++) {
@@ -250,7 +322,7 @@ namespace o2::hf_corrbkg
250322
break;
251323
case Pdg::kDStar:
252324
for (const auto& [flag, pdgCodes] : DStar::resoStatesDStarD0) {
253-
// std::cout << "Checking DStar resonant decay with flag: " << flag << ", pdgCodes: " << pdgCodes[0] << ", " << pdgCodes[1] << std::endl;
325+
std::cout << "Checking DStar resonant decay with flag: " << flag << ", pdgCodes: " << pdgCodes[0] << ", " << pdgCodes[1] << std::endl;
254326
if (checkResonantDecay(arrDaughIndex, pdgCodes)) {
255327
*channel = flag;
256328
LOG(info) << "Dstar resonant decay found with channel: " << static_cast<int>(*channel);

0 commit comments

Comments
 (0)