Skip to content

Commit bc22810

Browse files
FDUEnrichFDUEnrich
authored
Add GeneratorHF_D2H_ccbar_and_bbbar_PbPb_replaced teat
Test if the replacements and the decay model are correct.
1 parent 7eb8a8c commit bc22810

File tree

1 file changed

+215
-0
lines changed

1 file changed

+215
-0
lines changed

‎MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_PbPb_replaced.C‎

Lines changed: 215 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,215 @@
1+
int External() {
2+
std::string path{"o2sim_Kine.root"};
3+
//std::string path{"tf1/sgn_Kine.root"};
4+
5+
// Particle replacement configuration: {original, replacement}, frequencies
6+
std::array<std::array<int, 2>, 3> pdgReplParticles = {
7+
std::array{423, 4132}, // D*0 -> Xic0
8+
std::array{423, 4232}, // D*0 -> Xic+
9+
std::array{4212, 4332} // Sigmac+ -> Omegac0
10+
};
11+
std::array<std::array<int, 2>, 3> pdgReplPartCounters = {
12+
std::array{0, 0},
13+
std::array{0, 0},
14+
std::array{0, 0}
15+
};
16+
std::array<float, 3> freqRepl = {0.5, 0.5, 1.0};
17+
std::map<int, int> sumOrigReplacedParticles = {{423, 0}, {4212, 0}};
18+
19+
// Signal hadrons to check (only final charm baryons after replacement)
20+
std::array<int, 4> checkPdgHadron{4122, 4132, 4232, 4332};
21+
22+
// Expected decay channels - will be sorted automatically
23+
// Define both particle and antiparticle versions
24+
std::map<int, std::vector<std::vector<int>>> checkHadronDecaysRaw{
25+
// Λc+ decays (from cfg: 4122:addChannel + resonance decays)
26+
{4122, {
27+
{2212, 311}, {-2212, -311}, // p K0s
28+
{2212, -321, 211}, {-2212, 321, -211}, // p K- π+ (non-resonant)
29+
{2212, 313}, {-2212, -313}, // p K*0 (not decayed)
30+
{2212, 321, 211}, {-2212, -321, -211}, // p K*0 -> p (K- π+) [K*0 decayed]
31+
{2224, -321}, {-2224, 321}, // Delta++ K- (not decayed)
32+
{2212, 211, -321}, {-2212, -211, 321}, // Delta++ K- -> (p π+) K- [Delta decayed]
33+
{102134, 211}, {-102134, -211}, // Lambda(1520) π+ (not decayed)
34+
{2212, 321, 211}, {-2212, -321, -211}, // Lambda(1520) π+ -> (p K-) π+ [Lambda* decayed]
35+
{2212, -321, 211, 111}, {-2212, 321, -211, 111}, // p K- π+ π0
36+
{2212, -211, 211}, {-2212, 211, -211}, // p π- π+ (cfg line 61: 2212 -211 211)
37+
{2212, 333}, {-2212, 333}, // p φ (not decayed)
38+
{2212, 321, -321}, {-2212, -321, 321} // p φ -> p (K+ K-) [φ decayed]
39+
}},
40+
// Ξc0 decays (from cfg: 4132:onIfMatch)
41+
{4132, {
42+
{3312, 211}, {-3312, -211}, // Ξ- π+
43+
{3334, 321}, {-3334, -321} // Ω- K+
44+
}},
45+
// Ξc+ decays (from cfg: 4232:onIfMatch + resonance decays)
46+
{4232, {
47+
{2212, -321, 211}, {-2212, 321, -211}, // p K- π+
48+
{2212, -313}, {-2212, 313}, // p K̄*0 (not decayed)
49+
{2212, -321, 211}, {-2212, 321, -211}, // p K̄*0 -> p (K+ π-) [K*0 decayed]
50+
{2212, 333}, {-2212, 333}, // p φ (not decayed)
51+
{2212, 321, -321}, {-2212, -321, 321}, // p φ -> p (K+ K-) [φ decayed]
52+
{3222, -211, 211}, {-3222, 211, -211}, // Σ+ π- π+
53+
{3324, 211}, {-3324, -211}, // Ξ*0 π+
54+
{3312, 211, 211}, {-3312, -211, -211} // Ξ- π+ π+
55+
}},
56+
// Ωc0 decays (from cfg: 4332:onIfMatch)
57+
{4332, {
58+
{3334, 211}, {-3334, -211}, // Ω- π+
59+
{3312, 211}, {-3312, -211}, // Ξ- π+
60+
{3334, 321}, {-3334, -321} // Ω- K+
61+
}}
62+
};
63+
64+
// Sort all decay channels
65+
std::map<int, std::vector<std::vector<int>>> checkHadronDecays;
66+
for (auto &[pdg, decays] : checkHadronDecaysRaw) {
67+
for (auto decay : decays) {
68+
std::sort(decay.begin(), decay.end());
69+
checkHadronDecays[pdg].push_back(decay);
70+
}
71+
}
72+
73+
TFile file(path.c_str(), "READ");
74+
if (file.IsZombie())
75+
{
76+
std::cerr << "Cannot open ROOT file " << path << "\n";
77+
return 1;
78+
}
79+
80+
auto tree = (TTree *)file.Get("o2sim");
81+
std::vector<o2::MCTrack> *tracks{};
82+
tree->SetBranchAddress("MCTrack", &tracks);
83+
84+
int nSignals{}, nSignalGoodDecay{};
85+
std::map<int, int> failedDecayCount{{4122, 0}, {4132, 0}, {4232, 0}, {4332, 0}};
86+
std::map<int, std::set<std::vector<int>>> unknownDecays;
87+
auto nEvents = tree->GetEntries();
88+
89+
for (int i = 0; i < nEvents; i++)
90+
{
91+
tree->GetEntry(i);
92+
for (auto &track : *tracks)
93+
{
94+
auto pdg = track.GetPdgCode();
95+
auto absPdg = std::abs(pdg);
96+
97+
if (std::find(checkPdgHadron.begin(), checkPdgHadron.end(), absPdg) != checkPdgHadron.end())
98+
{
99+
nSignals++; // count signal PDG
100+
101+
// Count replacement particles (single-match per track)
102+
int matchedIdx = -1;
103+
bool matchedIsReplacement = false;
104+
for (int iRepl{0}; iRepl < 3; ++iRepl) {
105+
if (absPdg == pdgReplParticles[iRepl][0]) {
106+
matchedIdx = iRepl;
107+
matchedIsReplacement = false;
108+
break;
109+
}
110+
if (absPdg == pdgReplParticles[iRepl][1]) {
111+
matchedIdx = iRepl;
112+
matchedIsReplacement = true;
113+
break;
114+
}
115+
}
116+
if (matchedIdx >= 0) {
117+
if (matchedIsReplacement) {
118+
pdgReplPartCounters[matchedIdx][1]++;
119+
} else {
120+
pdgReplPartCounters[matchedIdx][0]++;
121+
}
122+
// Count the original-particle population once for this matched group
123+
sumOrigReplacedParticles[pdgReplParticles[matchedIdx][0]]++;
124+
}
125+
126+
// Collect decay products
127+
std::vector<int> pdgsDecay{};
128+
if (track.getFirstDaughterTrackId() >= 0 && track.getLastDaughterTrackId() >= 0) {
129+
for (int j{track.getFirstDaughterTrackId()}; j <= track.getLastDaughterTrackId(); ++j) {
130+
auto pdgDau = tracks->at(j).GetPdgCode();
131+
pdgsDecay.push_back(pdgDau);
132+
}
133+
}
134+
135+
std::sort(pdgsDecay.begin(), pdgsDecay.end());
136+
137+
// Check if decay matches expected channels
138+
bool foundMatch = false;
139+
for (auto &decay : checkHadronDecays[absPdg]) {
140+
if (pdgsDecay == decay) {
141+
nSignalGoodDecay++;
142+
foundMatch = true;
143+
break;
144+
}
145+
}
146+
147+
// Record failed decays for debugging
148+
if (!foundMatch && pdgsDecay.size() > 0) {
149+
failedDecayCount[absPdg]++;
150+
unknownDecays[absPdg].insert(pdgsDecay);
151+
}
152+
}
153+
}
154+
}
155+
156+
std::cout << "--------------------------------\n";
157+
std::cout << "# Events: " << nEvents << "\n";
158+
std::cout << "# signal charm baryons: " << nSignals << "\n";
159+
std::cout << "# signal charm baryons decaying in the correct channel: " << nSignalGoodDecay << "\n";
160+
161+
// Print failed decay statistics
162+
std::cout << "\nFailed decay counts:\n";
163+
for (auto &[pdg, count] : failedDecayCount) {
164+
if (count > 0) {
165+
std::cout << "PDG " << pdg << ": " << count << " failed decays\n";
166+
std::cout << " Unknown decay channels (first 5):\n";
167+
int printed = 0;
168+
for (auto &decay : unknownDecays[pdg]) {
169+
if (printed++ >= 5) break;
170+
std::cout << " [";
171+
for (size_t i = 0; i < decay.size(); ++i) {
172+
std::cout << decay[i];
173+
if (i < decay.size()-1) std::cout << ", ";
174+
}
175+
std::cout << "]\n";
176+
}
177+
}
178+
}
179+
std::cout << "\n";
180+
181+
std::cout << "# D*0 (original): " << pdgReplPartCounters[0][0] << "\n";
182+
std::cout << "# Xic0 (replaced from D*0): " << pdgReplPartCounters[0][1] << "\n";
183+
std::cout << "# Xic+ (replaced from D*0): " << pdgReplPartCounters[1][1] << "\n";
184+
std::cout << "# Sigmac+ (original): " << pdgReplPartCounters[2][0] << "\n";
185+
std::cout << "# Omegac0 (replaced from Sigmac+): " << pdgReplPartCounters[2][1] << "\n";
186+
187+
// Check forced decay fraction
188+
float fracForcedDecays = nSignals ? float(nSignalGoodDecay) / nSignals : 0.0f;
189+
std::cout << "# fraction of signals decaying into the correct channel: " << fracForcedDecays
190+
<< " (" << fracForcedDecays * 100.0f << "%)\n";
191+
if (fracForcedDecays < 0.9f) { // 90% threshold with tolerance
192+
std::cerr << "Fraction of signals decaying into the correct channel " << fracForcedDecays << " lower than expected\n";
193+
return 1;
194+
}
195+
196+
// Check particle replacement ratios (2-sigma statistical compatibility)
197+
for (int iRepl{0}; iRepl < 3; ++iRepl) {
198+
if (sumOrigReplacedParticles[pdgReplParticles[iRepl][0]] == 0) {
199+
continue; // Skip if no original particles found
200+
}
201+
202+
float expectedCount = freqRepl[iRepl] * sumOrigReplacedParticles[pdgReplParticles[iRepl][0]];
203+
float sigma = std::sqrt(freqRepl[iRepl] * sumOrigReplacedParticles[pdgReplParticles[iRepl][0]]);
204+
205+
if (std::abs(pdgReplPartCounters[iRepl][1] - expectedCount) > 2 * sigma) {
206+
float fracMeas = float(pdgReplPartCounters[iRepl][1]) / sumOrigReplacedParticles[pdgReplParticles[iRepl][0]];
207+
std::cerr << "Fraction of replaced " << pdgReplParticles[iRepl][0]
208+
<< " into " << pdgReplParticles[iRepl][1]
209+
<< " is " << fracMeas << " (expected " << freqRepl[iRepl] << ")\n";
210+
return 1;
211+
}
212+
}
213+
214+
return 0;
215+
}

0 commit comments

Comments
 (0)