Skip to content

Commit 803a34a

Browse files
MarcellocostiMarcellocosti
committed
Fix Dstar matching and add other channels
1 parent 795f547 commit 803a34a

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
@@ -25,6 +25,7 @@
2525
#include <tuple> // std::apply
2626
#include <utility> // std::move
2727
#include <vector> // std::vector
28+
#include <iostream>
2829

2930
// ROOT includes
3031
#include <TMCProcess.h> // for VMC Particle Production Process
@@ -918,7 +919,11 @@ struct RecoDecay {
918919
}
919920
// Check the PDG code of the particle.
920921
auto pdgCandidate = candidate.pdgCode();
922+
bool print = false; // Print debug messages only if print is true.
921923
// Printf("MC Gen: Candidate PDG: %d", pdgCandidate);
924+
if (std::abs(pdgCandidate) == 413) {
925+
print = true; // Print debug messages for D* candidates.
926+
}
922927
if (pdgCandidate == pdgParticle) { // exact PDG match
923928
sgn = 1;
924929
} else if (acceptAntiParticles && pdgCandidate == -pdgParticle) { // antiparticle PDG match
@@ -929,29 +934,56 @@ struct RecoDecay {
929934
}
930935
// Check the PDG codes of the decay products.
931936
if (N > 0) {
937+
if (print) {
938+
std::cout << "MC Gen: Checking decay products of " << N << " daughters" << std::endl;
939+
// Printf("MC Gen: Checking decay products of %ld daughters", N);
940+
}
932941
// Printf("MC Gen: Checking %d daughters", N);
933942
std::vector<int> arrAllDaughtersIndex; // vector of indices of all daughters
934943
// Check the daughter indices.
935944
if (!candidate.has_daughters()) {
945+
if (print) {
946+
std::cout << "MC Gen: Rejected: bad daughter index range: " << candidate.daughtersIds().front() << "-" << candidate.daughtersIds().back() << std::endl;
947+
// Printf("MC Gen: Rejected: bad daughter index range: %d-%d", candidate.daughtersIds().front(), candidate.daughtersIds().back());
948+
}
936949
// Printf("MC Gen: Rejected: bad daughter index range: %d-%d", candidate.daughtersIds().front(), candidate.daughtersIds().back());
937950
return false;
938951
}
939952
// Check that the number of direct daughters is not larger than the number of expected final daughters.
940953
if constexpr (!checkProcess) {
941954
if (candidate.daughtersIds().back() - candidate.daughtersIds().front() + 1 > static_cast<int>(N)) {
955+
if (print) {
956+
std::cout << "MC Gen: Rejected: too many direct daughters: " << candidate.daughtersIds().back() - candidate.daughtersIds().front() + 1
957+
<< " (expected " << N << " final)" << std::endl;
958+
// Printf("MC Gen: Rejected: too many direct daughters: %d (expected %ld final)", candidate.daughtersIds().back() - candidate.daughtersIds().front() + 1, N);
959+
}
942960
// Printf("MC Gen: Rejected: too many direct daughters: %d (expected %ld final)", candidate.daughtersIds().back() - candidate.daughtersIds().front() + 1, N);
943961
return false;
944962
}
945963
}
946964
// Get the list of actual final daughters.
947965
getDaughters<checkProcess>(candidate, &arrAllDaughtersIndex, arrPdgDaughters, depthMax);
966+
if (print) {
967+
std::cout << "MC Gen: Mother " << candidate.globalIndex() << " has " << arrAllDaughtersIndex.size() << " final daughters:";
968+
// Printf("MC Gen: Mother %ld has %ld final daughters:", candidate.globalIndex(), arrAllDaughtersIndex.size());
969+
for (auto i : arrAllDaughtersIndex) {
970+
std::cout << " " << i;
971+
// Printf(" %d", i);
972+
}
973+
std::cout << std::endl;
974+
// Printf("\n");
975+
}
948976
// printf("MC Gen: Mother %ld has %ld final daughters:", candidate.globalIndex(), arrAllDaughtersIndex.size());
949977
// for (auto i : arrAllDaughtersIndex) {
950978
// printf(" %d", i);
951979
// }
952980
// printf("\n");
953981
// Check whether the number of final daughters is equal to the required number.
954982
if (arrAllDaughtersIndex.size() != N) {
983+
if (print) {
984+
std::cout << "MC Gen: Rejected: incorrect number of final daughters: " << arrAllDaughtersIndex.size() << " (expected " << N << ")" << std::endl;
985+
// Printf("MC Gen: Rejected: incorrect number of final daughters: %ld (expected %ld)", arrAllDaughtersIndex.size(), N);
986+
}
955987
// Printf("MC Gen: Rejected: incorrect number of final daughters: %ld (expected %ld)", arrAllDaughtersIndex.size(), N);
956988
return false;
957989
}
@@ -969,6 +1001,10 @@ struct RecoDecay {
9691001
for (auto indexDaughterI : arrAllDaughtersIndex) { // o2-linter: disable=const-ref-in-for-loop (int elements)
9701002
auto candidateDaughterI = particlesMC.rawIteratorAt(indexDaughterI - particlesMC.offset()); // ith daughter particle
9711003
auto pdgCandidateDaughterI = candidateDaughterI.pdgCode(); // PDG code of the ith daughter
1004+
if (print) {
1005+
std::cout << "MC Gen: Daughter " << indexDaughterI << " PDG: " << pdgCandidateDaughterI << std::endl;
1006+
// Printf("MC Gen: Daughter %d PDG: %d", indexDaughterI, pdgCandidateDaughterI);
1007+
}
9721008
// Printf("MC Gen: Daughter %d PDG: %d", indexDaughterI, pdgCandidateDaughterI);
9731009
bool isPdgFound = false; // Is the PDG code of this daughter among the remaining expected PDG codes?
9741010
for (std::size_t iProngCp = 0; iProngCp < N; ++iProngCp) {
@@ -987,6 +1023,10 @@ struct RecoDecay {
9871023
*listIndexDaughters = arrAllDaughtersIndex;
9881024
}
9891025
}
1026+
if (print) {
1027+
std::cout << "MC Gen: Accepted: m: " << candidate.globalIndex() << std::endl;
1028+
// Printf("MC Gen: Accepted: m: %ld", candidate.globalIndex());
1029+
}
9901030
// Printf("MC Gen: Accepted: m: %d", candidate.globalIndex());
9911031
if (sign) {
9921032
*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)