Skip to content

Commit 2f79661

Browse files
iarseneIonut Cristian Arsene
iarsene
and
Ionut Cristian Arsene
authored
[PWGDQ] Changes to MCProng/MCSignal and dq-effieciency for generator level signal counting (AliceO2Group#10510)
Co-authored-by: Ionut Cristian Arsene <iarsene@cern.ch>
1 parent 6440e7d commit 2f79661

File tree

7 files changed

+198
-61
lines changed

7 files changed

+198
-61
lines changed

‎PWGDQ/Core/MCProng.cxx‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -21,7 +21,8 @@ std::map<TString, int> MCProng::fgSourceNames = {
2121
{"kPhysicalPrimary", MCProng::kPhysicalPrimary},
2222
{"kProducedInTransport", MCProng::kProducedInTransport},
2323
{"kProducedByGenerator", MCProng::kProducedByGenerator},
24-
{"kFromBackgroundEvent", MCProng::kFromBackgroundEvent}};
24+
{"kFromBackgroundEvent", MCProng::kFromBackgroundEvent},
25+
{"kHEPMCFinalState", MCProng::kHEPMCFinalState}};
2526

2627
//________________________________________________________________________________________________________________
2728
MCProng::MCProng() : fNGenerations(0),

‎PWGDQ/Core/MCProng.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -72,6 +72,7 @@ class MCProng
7272
kProducedInTransport, // Produced during transport through the detector (e.g. GEANT)
7373
kProducedByGenerator, // Produced by generator (if not, then produced by GEANT)
7474
kFromBackgroundEvent, // Produced in the underlying event
75+
kHEPMCFinalState, // HEPMC code 11
7576
kNSources
7677
};
7778

‎PWGDQ/Core/MCSignal.h‎

Lines changed: 66 additions & 56 deletions
Original file line numberDiff line numberDiff line change
@@ -216,41 +216,50 @@ bool MCSignal::CheckProng(int i, bool checkSources, const T& track)
216216
currentMCParticle = track;
217217
for (int j = 0; j < fProngs[i].fNGenerations; j++) {
218218
// check whether sources are required for this generation
219-
if (!fProngs[i].fSourceBits[j]) {
220-
continue;
219+
bool hasSources = false;
220+
if (fProngs[i].fSourceBits[j]) {
221+
hasSources = true;
221222
}
222223
// check each source
223224
uint64_t sourcesDecision = 0;
224-
// Check kPhysicalPrimary
225-
if (fProngs[i].fSourceBits[j] & (static_cast<uint64_t>(1) << MCProng::kPhysicalPrimary)) {
226-
if ((fProngs[i].fExcludeSource[j] & (static_cast<uint64_t>(1) << MCProng::kPhysicalPrimary)) != currentMCParticle.isPhysicalPrimary()) {
227-
sourcesDecision |= (static_cast<uint64_t>(1) << MCProng::kPhysicalPrimary);
225+
if (hasSources) {
226+
// Check kPhysicalPrimary
227+
if (fProngs[i].fSourceBits[j] & (static_cast<uint64_t>(1) << MCProng::kPhysicalPrimary)) {
228+
if ((fProngs[i].fExcludeSource[j] & (static_cast<uint64_t>(1) << MCProng::kPhysicalPrimary)) != currentMCParticle.isPhysicalPrimary()) {
229+
sourcesDecision |= (static_cast<uint64_t>(1) << MCProng::kPhysicalPrimary);
230+
}
228231
}
229-
}
230-
// Check kProducedInTransport
231-
if (fProngs[i].fSourceBits[j] & (static_cast<uint64_t>(1) << MCProng::kProducedInTransport)) {
232-
if ((fProngs[i].fExcludeSource[j] & (static_cast<uint64_t>(1) << MCProng::kProducedInTransport)) != (!currentMCParticle.producedByGenerator())) {
233-
sourcesDecision |= (static_cast<uint64_t>(1) << MCProng::kProducedInTransport);
232+
// Check kProducedInTransport
233+
if (fProngs[i].fSourceBits[j] & (static_cast<uint64_t>(1) << MCProng::kProducedInTransport)) {
234+
if ((fProngs[i].fExcludeSource[j] & (static_cast<uint64_t>(1) << MCProng::kProducedInTransport)) != (!currentMCParticle.producedByGenerator())) {
235+
sourcesDecision |= (static_cast<uint64_t>(1) << MCProng::kProducedInTransport);
236+
}
234237
}
235-
}
236-
// Check kProducedByGenerator
237-
if (fProngs[i].fSourceBits[j] & (static_cast<uint64_t>(1) << MCProng::kProducedByGenerator)) {
238-
if ((fProngs[i].fExcludeSource[j] & (static_cast<uint64_t>(1) << MCProng::kProducedByGenerator)) != currentMCParticle.producedByGenerator()) {
239-
sourcesDecision |= (static_cast<uint64_t>(1) << MCProng::kProducedByGenerator);
238+
// Check kProducedByGenerator
239+
if (fProngs[i].fSourceBits[j] & (static_cast<uint64_t>(1) << MCProng::kProducedByGenerator)) {
240+
if ((fProngs[i].fExcludeSource[j] & (static_cast<uint64_t>(1) << MCProng::kProducedByGenerator)) != currentMCParticle.producedByGenerator()) {
241+
sourcesDecision |= (static_cast<uint64_t>(1) << MCProng::kProducedByGenerator);
242+
}
240243
}
241-
}
242-
// Check kFromBackgroundEvent
243-
if (fProngs[i].fSourceBits[j] & (static_cast<uint64_t>(1) << MCProng::kFromBackgroundEvent)) {
244-
if ((fProngs[i].fExcludeSource[j] & (static_cast<uint64_t>(1) << MCProng::kFromBackgroundEvent)) != currentMCParticle.fromBackgroundEvent()) {
245-
sourcesDecision |= (static_cast<uint64_t>(1) << MCProng::kFromBackgroundEvent);
244+
// Check kFromBackgroundEvent
245+
if (fProngs[i].fSourceBits[j] & (static_cast<uint64_t>(1) << MCProng::kFromBackgroundEvent)) {
246+
if ((fProngs[i].fExcludeSource[j] & (static_cast<uint64_t>(1) << MCProng::kFromBackgroundEvent)) != currentMCParticle.fromBackgroundEvent()) {
247+
sourcesDecision |= (static_cast<uint64_t>(1) << MCProng::kFromBackgroundEvent);
248+
}
246249
}
247-
}
250+
// Check HEPMC code 11 (final state)
251+
if (fProngs[i].fSourceBits[j] & (static_cast<uint64_t>(1) << MCProng::kHEPMCFinalState)) {
252+
if ((fProngs[i].fExcludeSource[j] & (static_cast<uint64_t>(1) << MCProng::kHEPMCFinalState)) != (currentMCParticle.getHepMCStatusCode() == 11)) {
253+
sourcesDecision |= (static_cast<uint64_t>(1) << MCProng::kHEPMCFinalState);
254+
}
255+
}
256+
} // end if(hasSources)
248257
// no source bit is fulfilled
249-
if (!sourcesDecision) {
258+
if (hasSources && !sourcesDecision) {
250259
return false;
251260
}
252261
// if fUseANDonSourceBitMap is on, request all bits
253-
if (fProngs[i].fUseANDonSourceBitMap[j] && (sourcesDecision != fProngs[i].fSourceBits[j])) {
262+
if (hasSources && (fProngs[i].fUseANDonSourceBitMap[j] && (sourcesDecision != fProngs[i].fSourceBits[j]))) {
254263
return false;
255264
}
256265

@@ -280,8 +289,8 @@ bool MCSignal::CheckProng(int i, bool checkSources, const T& track)
280289
}
281290
}
282291
}
283-
}
284-
}
292+
} // end loop over generations
293+
} // end if(checkSources)
285294

286295
if (fProngs[i].fPDGInHistory.size() == 0) {
287296
return true;
@@ -299,40 +308,41 @@ bool MCSignal::CheckProng(int i, bool checkSources, const T& track)
299308
// Note: Currently no need to check generation InTime, so disable if case and always check BackInTime (direction of mothers)
300309
// The option to check for daughter in decay chain is still implemented but commented out.
301310

302-
// if (!fProngs[i].fCheckGenerationsInTime) { // check generation back in time
303-
while (currentMCParticle.has_mothers()) {
304-
auto mother = currentMCParticle.template mothers_first_as<P>();
305-
if (!fProngs[i].fExcludePDGInHistory[k] && fProngs[i].ComparePDG(mother.pdgCode(), fProngs[i].fPDGInHistory[k], true, fProngs[i].fExcludePDGInHistory[k])) {
306-
pdgInHistory.emplace_back(mother.pdgCode());
307-
break;
311+
if (!fProngs[i].fCheckGenerationsInTime) { // check generation back in time
312+
while (currentMCParticle.has_mothers()) {
313+
auto mother = currentMCParticle.template mothers_first_as<P>();
314+
if (!fProngs[i].fExcludePDGInHistory[k] && fProngs[i].ComparePDG(mother.pdgCode(), fProngs[i].fPDGInHistory[k], true, fProngs[i].fExcludePDGInHistory[k])) {
315+
pdgInHistory.emplace_back(mother.pdgCode());
316+
break;
317+
}
318+
if (fProngs[i].fExcludePDGInHistory[k] && !fProngs[i].ComparePDG(mother.pdgCode(), fProngs[i].fPDGInHistory[k], true, fProngs[i].fExcludePDGInHistory[k])) {
319+
return false;
320+
}
321+
ith++;
322+
currentMCParticle = mother;
323+
if (ith > 10) { // need error message. Given pdg code was not found within 10 generations of the particles decay chain.
324+
break;
325+
}
308326
}
309-
if (fProngs[i].fExcludePDGInHistory[k] && !fProngs[i].ComparePDG(mother.pdgCode(), fProngs[i].fPDGInHistory[k], true, fProngs[i].fExcludePDGInHistory[k])) {
327+
} /*else { // check generation in time
328+
if (!currentMCParticle.has_daughters()) {
310329
return false;
311330
}
312-
ith++;
313-
currentMCParticle = mother;
314-
if (ith > 10) { // need error message. Given pdg code was not found within 10 generations of the particles decay chain.
315-
break;
331+
const auto& daughtersSlice = currentMCParticle.template daughters_as<P>();
332+
for (auto& d : daughtersSlice) {
333+
if (!fProngs[i].fExcludePDGInHistory[k] && fProngs[i].ComparePDG(d.pdgCode(), fProngs[i].fPDGInHistory[k], true, fProngs[i].fExcludePDGInHistory[k])) {
334+
pdgInHistory.emplace_back(d.pdgCode());
335+
break;
336+
}
337+
if (fProngs[i].fExcludePDGInHistory[k] && !fProngs[i].ComparePDG(d.pdgCode(), fProngs[i].fPDGInHistory[k], true, fProngs[i].fExcludePDGInHistory[k])) {
338+
return false;
339+
}
340+
ith++;
341+
if (ith > 10) { // need error message. Given pdg code was not found within 10 generations of the particles decay chain.
342+
break;
343+
}
316344
}
317-
}
318-
// } else { // check generation in time
319-
// if (!currentMCParticle.has_daughters())
320-
// return false;
321-
// const auto& daughtersSlice = currentMCParticle.template daughters_as<P>();
322-
// for (auto& d : daughtersSlice) {
323-
// if (!fProngs[i].fExcludePDGInHistory[k] && fProngs[i].ComparePDG(d.pdgCode(), fProngs[i].fPDGInHistory[k], true, fProngs[i].fExcludePDGInHistory[k])) {
324-
// pdgInHistory.emplace_back(d.pdgCode());
325-
// break;
326-
// }
327-
// if (fProngs[i].fExcludePDGInHistory[k] && !fProngs[i].ComparePDG(d.pdgCode(), fProngs[i].fPDGInHistory[k], true, fProngs[i].fExcludePDGInHistory[k])) {
328-
// return false;
329-
// }
330-
// ith++;
331-
// if (ith > 10) { // need error message. Given pdg code was not found within 10 generations of the particles decay chain.
332-
// break;
333-
// }
334-
// }
335-
// }
345+
}*/
336346
}
337347
if (pdgInHistory.size() != nIncludedPDG) { // vector has as many entries as mothers (daughters) defined for prong
338348
return false;

‎PWGDQ/Core/MCSignalLibrary.cxx‎

Lines changed: 31 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -205,11 +205,23 @@ MCSignal* o2::aod::dqmcsignals::GetMCSignal(const char* name)
205205
signal = new MCSignal(name, "All beauty hadrons", {prong}, {-1});
206206
return signal;
207207
}
208+
if (!nameStr.compare("allBeautyHadronsFS")) {
209+
MCProng prong(1, {503}, {true}, {false}, {0}, {0}, {false});
210+
prong.SetSourceBit(0, MCProng::kHEPMCFinalState);
211+
signal = new MCSignal(name, "All beauty hadrons", {prong}, {-1});
212+
return signal;
213+
}
208214
if (!nameStr.compare("allOpenBeautyHadrons")) {
209215
MCProng prong(1, {502}, {true}, {false}, {0}, {0}, {false});
210216
signal = new MCSignal(name, "All open beauty hadrons", {prong}, {-1});
211217
return signal;
212218
}
219+
if (!nameStr.compare("allOpenBeautyHadronsFS")) {
220+
MCProng prong(1, {502}, {true}, {false}, {0}, {0}, {false});
221+
prong.SetSourceBit(0, MCProng::kHEPMCFinalState);
222+
signal = new MCSignal(name, "All open beauty hadrons", {prong}, {-1});
223+
return signal;
224+
}
213225
if (!nameStr.compare("Bc")) {
214226
MCProng prong(1, {541}, {true}, {false}, {0}, {0}, {false});
215227
signal = new MCSignal(name, "Bc", {prong}, {-1});
@@ -236,11 +248,23 @@ MCSignal* o2::aod::dqmcsignals::GetMCSignal(const char* name)
236248
signal = new MCSignal(name, "Everything from beauty", {prong}, {-1});
237249
return signal;
238250
}
251+
if (!nameStr.compare("everythingFromBeautyFS")) {
252+
MCProng prong(2, {0, 503}, {true, true}, {false, false}, {0, 0}, {0, 0}, {false, false});
253+
prong.SetSourceBit(1, MCProng::kHEPMCFinalState);
254+
signal = new MCSignal(name, "Everything from beauty", {prong}, {-1});
255+
return signal;
256+
}
239257
if (!nameStr.compare("everythingFromEverythingFromBeauty")) {
240258
MCProng prong(3, {0, 0, 503}, {true, true, true}, {false, false, false}, {0, 0, 0}, {0, 0, 0}, {false, false, false});
241259
signal = new MCSignal(name, "Everything from everything from beauty", {prong}, {-1});
242260
return signal;
243261
}
262+
if (!nameStr.compare("everythingFromEverythingFromBeautyFS")) {
263+
MCProng prong(3, {0, 0, 503}, {true, true, true}, {false, false, false}, {0, 0, 0}, {0, 0, 0}, {false, false, false});
264+
prong.SetSourceBit(2, MCProng::kHEPMCFinalState);
265+
signal = new MCSignal(name, "Everything from everything from beauty", {prong}, {-1});
266+
return signal;
267+
}
244268
if (!nameStr.compare("allCharmHadrons")) {
245269
MCProng prong(1, {403}, {true}, {false}, {0}, {0}, {false});
246270
signal = new MCSignal(name, "All charm hadrons", {prong}, {-1});
@@ -1249,6 +1273,13 @@ MCSignal* o2::aod::dqmcsignals::GetMCSignal(const char* name)
12491273
return signal;
12501274
}
12511275

1276+
if (!nameStr.compare("BplusFS")) {
1277+
MCProng prong(1, {521}, {true}, {false}, {0}, {0}, {false});
1278+
prong.SetSourceBit(0, MCProng::kHEPMCFinalState);
1279+
signal = new MCSignal(name, "B+", {prong}, {-1});
1280+
return signal;
1281+
}
1282+
12521283
if (!nameStr.compare("beautyPairs")) {
12531284
MCProng prong(1, {503}, {true}, {false}, {0}, {0}, {false});
12541285
signal = new MCSignal("beautyPairs", "Beauty hadron pair", {prong, prong}, {-1, -1});

‎PWGDQ/DataModel/ReducedInfoTables.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -408,6 +408,7 @@ DECLARE_SOA_TABLE(ReducedMCTracks, "AOD", "REDUCEDMCTRACK", //! MC track inform
408408
mcparticle::FromBackgroundEvent<mcparticle::Flags>,
409409
mcparticle::GetGenStatusCode<mcparticle::Flags, mcparticle::StatusCode>,
410410
mcparticle::GetProcess<mcparticle::Flags, mcparticle::StatusCode>,
411+
mcparticle::GetHepMCStatusCode<mcparticle::Flags, mcparticle::StatusCode>,
411412
mcparticle::IsPhysicalPrimary<mcparticle::Flags>);
412413

413414
using ReducedMCTrack = ReducedMCTracks::iterator;

‎PWGDQ/TableProducer/tableMakerMC_withAssoc.cxx‎

Lines changed: 38 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -106,6 +106,14 @@ constexpr static uint32_t gkMuonFillMapWithCov = VarManager::ObjTypes::Muon | Va
106106
// constexpr static uint32_t gkTrackFillMapWithAmbi = VarManager::ObjTypes::Track | VarManager::ObjTypes::AmbiTrack;
107107
constexpr static uint32_t gkMFTFillMap = VarManager::ObjTypes::TrackMFT;
108108

109+
template <typename TMap>
110+
void PrintBitMap(TMap map, int nbits)
111+
{
112+
for (int i = 0; i < nbits; i++) {
113+
cout << ((map & (TMap(1) << i)) > 0 ? "1" : "0");
114+
}
115+
}
116+
109117
struct TableMakerMC {
110118

111119
Produces<ReducedMCEvents> eventMC;
@@ -457,6 +465,36 @@ struct TableMakerMC {
457465
}
458466
i++;
459467
}
468+
469+
/*if ((std::abs(mctrack.pdgCode())>400 && std::abs(mctrack.pdgCode())<599) ||
470+
(std::abs(mctrack.pdgCode())>4000 && std::abs(mctrack.pdgCode())<5999) ||
471+
(mcflags > 0)) {
472+
cout << ">>>>>>>>>>>>>>>>>>>>>>> track idx / pdg / process / status code / HEPMC status / primary : "
473+
<< mctrack.globalIndex() << " / " << mctrack.pdgCode() << " / "
474+
<< mctrack.getProcess() << " / " << mctrack.getGenStatusCode() << " / " << mctrack.getHepMCStatusCode() << " / " << mctrack.isPhysicalPrimary() << endl;
475+
cout << ">>>>>>>>>>>>>>>>>>>>>>> track bitmap: ";
476+
PrintBitMap(mcflags, 16);
477+
cout << endl;
478+
if (mctrack.has_mothers()) {
479+
for (auto& m : mctrack.mothersIds()) {
480+
if (m < mcTracks.size()) { // protect against bad mother indices
481+
auto aMother = mcTracks.rawIteratorAt(m);
482+
cout << "<<<<<< mother idx / pdg: " << m << " / " << aMother.pdgCode() << endl;
483+
}
484+
}
485+
}
486+
487+
if (mctrack.has_daughters()) {
488+
for (int d = mctrack.daughtersIds()[0]; d <= mctrack.daughtersIds()[1]; ++d) {
489+
490+
if (d < mcTracks.size()) { // protect against bad daughter indices
491+
auto aDaughter = mcTracks.rawIteratorAt(d);
492+
cout << "<<<<<< daughter idx / pdg: " << d << " / " << aDaughter.pdgCode() << endl;
493+
}
494+
}
495+
}
496+
}*/
497+
460498
// if no MC signals were matched, continue
461499
if (mcflags == 0) {
462500
continue;

0 commit comments

Comments
 (0)