13 #include <analysis/utility/MCMatching.h>
14 #include <analysis/utility/AnalysisConfiguration.h>
16 #include <analysis/dataobjects/Particle.h>
17 #include <mdst/dataobjects/MCParticle.h>
19 #include <framework/datastore/StoreArray.h>
20 #include <framework/logging/Logger.h>
22 #include <unordered_set>
31 if (flags == c_Correct)
39 case c_MissFSR : s +=
"c_MissFSR";
break;
40 case c_MissingResonance : s +=
"c_MissingResonance";
break;
41 case c_DecayInFlight : s +=
"c_DecayInFlight";
break;
42 case c_MissNeutrino : s +=
"c_MissNeutrino";
break;
43 case c_MissGamma : s +=
"c_MissGamma";
break;
44 case c_MissMassiveParticle : s +=
"c_MissMassiveParticle";
break;
45 case c_MissKlong : s +=
"c_MissKlong";
break;
46 case c_MisID : s +=
"c_MisID";
break;
47 case c_AddedWrongParticle : s +=
"c_AddedWrongParticle";
break;
48 case c_InternalError : s +=
"c_InternalError";
break;
49 case c_MissPHOTOS : s +=
"c_MissPHOTOS";
break;
50 case c_AddedRecoBremsPhoton : s +=
"c_AddedRecoBremsPhoton";
break;
53 B2ERROR(
"MCMatching::explainFlags() doesn't know about flag " << f <<
", please update it.");
67 genMCPMothers.push_back(mcP->
getIndex());
77 for (
unsigned int i = lastMother; i < firstMothers.size(); i++) {
78 if (firstMothers[i] == idx)
93 if (particle->hasExtraInfo(c_extraInfoMCErrors))
101 int nChildren = particle->getNDaughters();
102 if (nChildren == 0) {
108 bool daugMCTruth =
true;
109 for (
int i = 0; i < nChildren; ++i) {
110 const Particle* daugP = particle->getDaughter(i);
112 daugMCTruth &= setMCTruth(daugP);
118 if (nChildren == 1) {
120 const Particle* daugP = particle->getDaughter(0);
134 vector<int> firstDaugMothers;
137 for (
int i = 0; i < nChildren; ++i) {
138 const Particle* daugP = particle->getDaughter(i);
142 fillGenMothers(daugMCP, firstDaugMothers);
144 lastMother = findCommonMother(daugMCP, firstDaugMothers, lastMother);
145 if (lastMother == -1)
150 motherIndex = firstDaugMothers[lastMother];
163 B2ERROR(
"setMCTruth(): sanity check failed!");
167 const MCParticle* mcMatch = mcParticles[motherIndex - 1];
168 particle->addRelationTo(mcMatch);
177 void appendParticles(
const Particle* p, unordered_set<const MCParticle*>& mcMatchedParticles)
179 for (
unsigned i = 0; i < p->getNDaughters(); ++i) {
180 const Particle* daug = p->getDaughter(i);
185 mcMatchedParticles.insert(mcParticle);
187 if (daug->getNDaughters() != 0) {
189 appendParticles(daug, mcMatchedParticles);
200 mcMatchedParticles.insert(mcParticle);
213 void appendParticles(
const MCParticle* gen, vector<const MCParticle*>& children)
218 const vector<MCParticle*>& genDaughters = gen->
getDaughters();
219 for (
auto daug : genDaughters) {
220 children.push_back(daug);
221 appendParticles(daug, children);
228 const int property = part->getProperty();
230 const int absPDG = abs(mcDaug->
getPDG());
234 if (property & Particle::PropertyFlags::c_IsIgnoreIntermediate)
239 if (property & Particle::PropertyFlags::c_IsIgnoreRadiatedPhotons)
242 if (property & Particle::PropertyFlags::c_IsIgnoreGamma)
245 }
else if (absPDG == 12 or absPDG == 14 or absPDG == 16) {
246 if (property & Particle::PropertyFlags::c_IsIgnoreNeutrino)
249 if (property & Particle::PropertyFlags::c_IsIgnoreMassive)
257 void appendAcceptedMissingDaughters(
const Particle* p, unordered_set<const MCParticle*>& acceptedParticles)
261 vector<const MCParticle*> genDaughters;
262 appendParticles(mcParticle, genDaughters);
264 for (
auto mcDaug : genDaughters) {
265 if (isDaughterAccepted(mcDaug, p))
266 acceptedParticles.insert(mcDaug);
268 acceptedParticles.erase(mcDaug);
272 for (
unsigned i = 0; i < p->getNDaughters(); ++i) {
273 const Particle* daug = p->getDaughter(i);
274 appendAcceptedMissingDaughters(daug, acceptedParticles);
284 if (particle->hasExtraInfo(c_extraInfoMCErrors)) {
285 return particle->getExtraInfo(c_extraInfoMCErrors);
288 mcParticle = particle->getRelatedTo<
MCParticle>();
289 return setMCErrorsExtraInfo(
const_cast<Particle*
>(particle), mcParticle);
295 auto setStatus = [](
Particle * part,
int s) ->
int {
296 part->addExtraInfo(c_extraInfoMCErrors, s);
301 return setStatus(particle, MCErrorFlags::c_InternalError);
303 if (particle->getNDaughters() == 0) {
305 return setStatus(particle, getFlagsOfFSP(particle, mcParticle));;
318 if (mother and particle->getPDGCode() == mother->getPDG() and getNumberOfDaughtersWithoutNeutrinos(mother) == 1
319 and !particle->hasExtraInfo(
"bremsCorrected")) {
320 if (abs(mother->getPDG()) != 15 and abs(mcParticle->
getPDG()) != 15) {
321 B2WARNING(
"Special treatment in MCMatching for tau is called for a non-tau particle. Check if you discovered another special case here, or if we have a bug! "
322 << mother->getPDG() <<
" " << particle->getPDGCode() <<
" " << mcParticle->
getPDG());
325 if (not(particle->getProperty() & Particle::PropertyFlags::c_IsIgnoreIntermediate))
326 status |= MCErrorFlags::c_MissingResonance;
327 }
else if (!(particle->getProperty() & Particle::PropertyFlags::c_IsUnspecified)) {
329 status |= MCErrorFlags::c_AddedWrongParticle;
333 status |= getFlagsOfDaughters(particle, mcParticle);
335 status |= getMissingParticleFlags(particle, mcParticle);
338 status &= ~(getFlagsIgnoredByProperty(particle));
340 return setStatus(particle, status);
345 if (particle->getPDGCode() == mcParticle->
getPDG())
346 return MCErrorFlags::c_Correct;
351 return MCErrorFlags::c_MisID;
354 int status = MCErrorFlags::c_DecayInFlight;
362 status |= MCErrorFlags::c_InternalError;
363 }
else if (particle->getPDGCode() != primary->
getPDG()) {
365 status |= MCErrorFlags::c_MisID;
373 unsigned nChildren = particle->getNDaughters();
376 vector<const MCParticle*> genParts;
378 if (particle->hasExtraInfo(
"bremsCorrected") && nChildren > 1) {
379 if (mcParticle && mcParticle->
getMother())
380 appendParticles(mcParticle->
getMother(), genParts);
383 int daughterStatuses = 0;
384 vector<int> daughterProperties = particle->getDaughterProperties();
385 for (
unsigned i = 0; i < nChildren; ++i) {
386 const Particle* daughter = particle->getDaughter(i);
387 int daughterStatus = 0;
390 daughterStatus |= getFlagsOfBremsPhotonDaughter(daughter, mcParticle, genParts);
392 daughterStatus |= getMCErrors(daughter);
395 int daughterStatusAcceptMask = (~c_Correct);
396 if (i < daughterProperties.size())
397 daughterStatusAcceptMask = makeDaughterAcceptMask(daughterProperties[i]);
399 daughterStatuses |= (daughterStatus & daughterStatusAcceptMask);
403 const int daughterStatusesAcceptMask = c_MisID | c_AddedWrongParticle | c_DecayInFlight | c_InternalError | c_AddedRecoBremsPhoton;
404 return (daughterStatuses & daughterStatusesAcceptMask);
410 const vector<const MCParticle*>& genParts)
413 int daughterStatus = getMCErrors(daughter);
419 daughterStatus &= (~c_InternalError);
420 daughterStatus |= c_AddedRecoBremsPhoton;
423 else if (mcDaughter == mcParticle) {
424 daughterStatus &= (~c_MisID);
427 else if (std::find(genParts.begin(), genParts.end(), mcDaughter) == genParts.end()) {
428 daughterStatus |= c_AddedRecoBremsPhoton;
432 return daughterStatus;
438 unsigned int number_of_neutrinos = 0;
439 for (
auto& p : daughters) {
440 auto pdg = abs(p->getPDG());
441 if (pdg == 12 || pdg == 14 || pdg == 16) {
442 number_of_neutrinos++;
445 return daughters.size() - number_of_neutrinos;
482 B2ERROR(
"The hell? why would this gamma not have a mother?");
486 int ndaug = mother->getNDaughters();
506 vector<const MCParticle*> genParts;
507 appendParticles(mcParticle, genParts);
509 unordered_set<const MCParticle*> mcMatchedParticles;
510 appendParticles(particle, mcMatchedParticles);
511 unordered_set<const MCParticle*> acceptedParticles;
512 appendAcceptedMissingDaughters(particle, acceptedParticles);
513 for (
auto part : acceptedParticles) {
514 mcMatchedParticles.insert(part);
521 if (mcMatchedParticles.find(genPart) != mcMatchedParticles.end())
525 const int generatedPDG = genPart->getPDG();
526 const int absGeneratedPDG = abs(generatedPDG);
528 if (!
isFSP(generatedPDG)) {
529 flags |= c_MissingResonance;
532 if (flags & (c_MissFSR | c_MissGamma)) {
535 if (isFSRLegacy(genPart)) {
537 flags |= c_MissPHOTOS;
539 flags |= c_MissGamma;
546 flags |= c_MissPHOTOS;
548 flags |= c_MissGamma;
550 }
else if (absGeneratedPDG == 12 || absGeneratedPDG == 14 || absGeneratedPDG == 16) {
551 flags |= c_MissNeutrino;
553 flags |= c_MissMassiveParticle;
554 if (absGeneratedPDG == 130)
555 flags |= c_MissKlong;
563 unordered_set<const MCParticle*> mcMatchedParticles;
566 appendParticles(particle, mcMatchedParticles);
567 vector<const MCParticle*> genParts;
568 appendParticles(mcParticle, genParts);
570 int nMissingDaughter = 0;
573 const bool missing = (mcMatchedParticles.find(genPart) == mcMatchedParticles.end());
576 const int generatedPDG = genPart->getPDG();
577 const int absGeneratedPDG = abs(generatedPDG);
579 auto result = find(daughterPDG.begin(), daughterPDG.end(), absGeneratedPDG);
580 if (result != daughterPDG.end())
585 return nMissingDaughter;
601 if (daughterProperty & Particle::PropertyFlags::c_IsIgnoreMisID) flags |= (
MCMatching::c_MisID);