26 #include <HepMC/GenEvent.h>
27 #include <HepMC/GenParticle.h>
28 #include <HepMC/GenVertex.h>
29 #include <HepMC/IteratorRange.h>
30 #include <HepMC/SimpleVector.h>
32 #include <TDatabasePDG.h>
42 int listOfResonantPIDs[] = {111, 113, 213, 333, 310, 311, 313, 323, 413, 423, 513, 523, 441, 443, 100443, 9000111, 9000211, 100111, 100211, 10111,
43 10211, 9010111, 9010211, 10113, 10213, 20113, 20213, 9000113, 9000213, 100113, 100213, 9010113, 9010213, 9020113, 9020213,
44 30113, 30213, 9030113, 9030213, 9040113, 9040213, 115, 215, 10115, 10215, 9000115, 9000215, 9010115, 9010215, 117, 217,
45 9000117, 9000217, 9010117, 9010217, 119, 219, 221, 331, 9000221, 9010221, 100221, 10221, 9020221, 100331, 9030221, 10331,
46 9040221, 9050221, 9060221, 9070221, 9080221, 223, 10223, 20223, 10333, 20333, 1000223, 9000223, 9010223, 30223, 100333, 225,
47 9000225, 335, 9010225, 9020225, 10225, 9030225, 10335, 9040225, 9050225, 9060225, 9070225, 9080225, 9090225, 227, 337, 229,
48 9000229, 9010229, 9000311, 9000321, 10311, 10321, 100311, 100321, 9010311, 9010321, 9020311, 9020321, 10313, 10323, 20313,
49 20323, 100313, 100323, 9000313, 9000323, 30313, 30323, 315, 325, 9000315, 9000325, 10315, 10325, 20315, 20325, 9010315,
50 9010325, 9020315, 9020325, 317, 327, 9010317, 9010327, 319, 329, 9000319, 9000329, 10411, 10421, 10413, 10423, 20413, 20423,
51 415, 425, 431, 10431, 433, 10433, 20433, 435, 10511, 10521, 10513, 10523, 20513, 20523, 515, 525, 10531, 533, 10533, 20533, 535,
52 10541, 543, 10543, 20543, 545, 10441, 100441, 10443, 20443, 30443, 9000443, 9010443, 9020443, 445, 100445, 551, 10551, 100551,
53 110551, 200551, 210551, 10553, 20553, 30553, 110553, 120553, 130553, 210553, 220553, 9000553, 9010553, 555, 10555, 20555, 100555,
54 110555, 120555, 200555, 557, 100557, 2224, 2214, 2114, 1114, 3212, 3224, 3214, 3114, 3324, 3314, 4222, 4212, 4112, 4224, 4214,
55 4114, 4232, 4132, 4322, 4312, 4324, 4314, 4332, 4334, 4412, 4422, 4414, 4424, 4432, 4434, 4444};
71 std::cout <<
"DecayFinder name: " <<
Name() << std::endl;
79 return canSearchDecay;
112 bool ddCanBeParsed =
true;
114 size_t daughterLocator;
117 std::string intermediate;
118 std::string daughter;
120 int mother_charge = 0;
121 std::vector<int> intermediates_charge;
122 std::vector<int> daughters_charge;
124 std::string decayArrow =
"->";
125 std::string chargeIndicator =
"^";
126 std::string startIntermediate =
"{";
127 std::string endIntermediate =
"}";
133 while ((pos = manipulateDecayDescriptor.find(
" ")) != std::string::npos) manipulateDecayDescriptor.replace(pos, 1,
"");
136 std::string checkForCC = manipulateDecayDescriptor.substr(0, 1) + manipulateDecayDescriptor.substr(manipulateDecayDescriptor.size() - 3, 3);
137 std::for_each(checkForCC.begin(), checkForCC.end(), [](
char&
c) {
c = ::toupper(
c); });
140 if (checkForCC ==
"[]CC")
142 manipulateDecayDescriptor = manipulateDecayDescriptor.substr(1, manipulateDecayDescriptor.size() - 4);
147 size_t findMotherEndPoint = manipulateDecayDescriptor.find(decayArrow);
148 mother = manipulateDecayDescriptor.substr(0, findMotherEndPoint);
152 ddCanBeParsed =
false;
153 manipulateDecayDescriptor.erase(0, findMotherEndPoint + decayArrow.length());
156 while ((pos = manipulateDecayDescriptor.find(startIntermediate)) != std::string::npos)
158 size_t findIntermediateStartPoint = manipulateDecayDescriptor.find(startIntermediate, pos);
159 size_t findIntermediateEndPoint = manipulateDecayDescriptor.find(endIntermediate, pos);
160 std::string intermediateDecay = manipulateDecayDescriptor.substr(pos + 1, findIntermediateEndPoint - (pos + 1));
162 intermediate = intermediateDecay.substr(0, intermediateDecay.find(decayArrow));
168 ddCanBeParsed =
false;
172 intermediateDecay.erase(0, intermediateDecay.find(decayArrow) + decayArrow.length());
173 while ((daughterLocator = intermediateDecay.find(chargeIndicator)) != std::string::npos)
175 daughter = intermediateDecay.substr(0, daughterLocator);
176 daughter += intermediateDecay.substr(daughterLocator + 1, 1);
180 daughters_charge.push_back(
get_charge(daughter));
183 ddCanBeParsed =
false;
184 intermediateDecay.erase(0, daughterLocator + 2);
187 manipulateDecayDescriptor.erase(findIntermediateStartPoint, findIntermediateEndPoint + 1 - findIntermediateStartPoint);
193 while ((daughterLocator = manipulateDecayDescriptor.find(chargeIndicator)) != std::string::npos)
195 daughter = manipulateDecayDescriptor.substr(0, daughterLocator);
196 daughter += manipulateDecayDescriptor.substr(daughterLocator + 1, 1);
200 daughters_charge.push_back(
get_charge(daughter));
203 ddCanBeParsed =
false;
204 manipulateDecayDescriptor.erase(0, daughterLocator + 2);
208 unsigned int trackStart = 0;
209 unsigned int trackEnd = 0;
212 trackStart = trackEnd;
217 for (
unsigned int j = trackStart; j < trackEnd; ++j)
219 vtxCharge += daughters_charge[j];
222 intermediates_charge.push_back(vtxCharge);
225 for (
unsigned int i = 0; i <
m_daughters_ID.size(); ++i) mother_charge += daughters_charge[i];
246 if (
Verbosity() >=
VERBOSITY_SOME) std::cout <<
"Your decay descriptor cannot be parsed, " <<
Name() <<
" will not be registered" << std::endl;
253 bool decayWasFound =
false;
254 bool aTrackFailedPT =
false;
255 bool aTrackFailedETA =
false;
261 m_geneventmap = findNode::getClass<PHHepMCGenEventMap>(topNode,
"PHHepMCGenEventMap");
264 std::cout <<
"DecayFinder: Missing node PHHepMCGenEventMap" << std::endl;
271 std::cout <<
"DecayFinder: Missing node PHHepMCGenEvent" << std::endl;
277 std::vector<int> positive_motherDecayProducts;
281 for (HepMC::GenEvent::particle_const_iterator
p = theEvent->particles_begin();
p != theEvent->particles_end(); ++
p)
287 bool breakOut =
false;
288 std::vector<int> correctMotherProducts;
289 decayChain.push_back(std::make_pair((*p)->barcode(), (*p)->pdg_id()));
291 for (HepMC::GenVertex::particle_iterator children = (*p)->end_vertex()->particles_begin(HepMC::children);
292 children != (*p)->end_vertex()->particles_end(HepMC::children); ++children)
294 if (
Verbosity() >=
VERBOSITY_MAX) std::cout <<
"--children->pdg_id(): " << (*children)->pdg_id() << std::endl;
301 if (!
m_allowPi0 && (*children)->pdg_id() == 111)
308 if (std::find(positive_motherDecayProducts.begin(), positive_motherDecayProducts.end(),
309 std::abs((*children)->pdg_id())) != positive_motherDecayProducts.end())
313 if (needThisParticle)
315 correctMotherProducts.push_back((*children)->pdg_id());
316 decayChain.push_back(std::make_pair((*children)->barcode(), (*children)->pdg_id()));
323 for (HepMC::GenVertex::particle_iterator grandchildren = (*children)->end_vertex()->particles_begin(HepMC::children);
324 grandchildren != (*children)->end_vertex()->particles_end(HepMC::children); ++grandchildren)
327 if (needThisParticle)
329 correctMotherProducts.push_back((*grandchildren)->pdg_id());
330 decayChain.push_back(std::make_pair((*grandchildren)->barcode(), (*grandchildren)->pdg_id()));
347 std::cout <<
"Printing required mother decay products: ";
349 std::cout << std::endl;
350 std::cout <<
"Printing actual mother decay products: ";
351 for (
unsigned int i = 0; i < correctMotherProducts.size(); ++i) std::cout << correctMotherProducts[i] <<
", ";
352 std::cout << std::endl;
354 std::cout <<
"*\n* These vectors match\n*\n"
357 std::cout <<
"*\n* These vectors DONT match\n*\n"
369 std::cout <<
"Checking CC state" << std::endl;
371 std::cout <<
"*\n* These vectors match\n*\n"
374 std::cout <<
"*\n* These vectors DONT match\n*\n"
383 if (aTrackFailedPT && !aTrackFailedETA)
385 else if (!aTrackFailedPT && aTrackFailedETA)
387 else if (aTrackFailedPT && aTrackFailedETA)
397 return decayWasFound;
402 bool acceptParticle =
false;
405 std::vector<int> positive_intermediates_ID;
409 if (std::find(positive_intermediates_ID.begin(), positive_intermediates_ID.end(),
410 abs(particle->pdg_id())) != positive_intermediates_ID.end())
412 std::vector<int> requiredIntermediateDecayProducts;
413 std::vector<int> actualIntermediateDecayProducts;
415 auto it = std::find(positive_intermediates_ID.begin(), positive_intermediates_ID.end(),
abs(particle->pdg_id()));
418 unsigned int trackStart = 0, trackStop = 0;
427 for (
unsigned int i = trackStart; i < trackStop; ++i) requiredIntermediateDecayProducts.push_back(
m_daughters_ID[i]);
429 for (HepMC::GenVertex::particle_iterator grandchildren = particle->end_vertex()->particles_begin(HepMC::children);
430 grandchildren != particle->end_vertex()->particles_end(HepMC::children); ++grandchildren)
432 if (
Verbosity() >=
VERBOSITY_MAX) std::cout <<
"----grandchildren->pdg_id(): " << (*grandchildren)->pdg_id() << std::endl;
437 for (HepMC::GenVertex::particle_iterator greatgrandchildren = (*grandchildren)->end_vertex()->particles_begin(HepMC::children);
438 greatgrandchildren != (*grandchildren)->end_vertex()->particles_end(HepMC::children); ++greatgrandchildren)
440 if (
Verbosity() >=
VERBOSITY_MAX) std::cout <<
"--------greatgrandchildren->pdg_id(): " << (*greatgrandchildren)->pdg_id() << std::endl;
444 else if (!
m_allowPhotons && (*greatgrandchildren)->pdg_id() == 22)
449 else if (
m_allowPi0 && (*greatgrandchildren)->pdg_id() == 111)
451 else if (!
m_allowPi0 && (*greatgrandchildren)->pdg_id() == 111)
458 actualIntermediateDecayProducts.push_back((*greatgrandchildren)->pdg_id());
459 HepMC::FourVector myFourVector = (*greatgrandchildren)->momentum();
460 if (myFourVector.perp() < 0.2) trackFailedPT =
true;
461 if (
std::abs(myFourVector.eta()) > 1.1) trackFailedETA =
true;
472 else if (
m_allowPi0 && (*grandchildren)->pdg_id() == 111)
474 else if (!
m_allowPi0 && (*grandchildren)->pdg_id() == 111)
481 actualIntermediateDecayProducts.push_back((*grandchildren)->pdg_id());
482 HepMC::FourVector myFourVector = (*grandchildren)->momentum();
483 if (myFourVector.perp() < 0.2) trackFailedPT =
true;
484 if (
std::abs(myFourVector.eta()) > 1.1) trackFailedETA =
true;
493 std::cout <<
"Printing required intermediate decay products: ";
494 for (
unsigned int i = 0; i < requiredIntermediateDecayProducts.size(); ++i) std::cout << requiredIntermediateDecayProducts[i] <<
", ";
495 std::cout << std::endl;
496 std::cout <<
"Printing actual intermediate decay products: ";
497 for (
unsigned int i = 0; i < actualIntermediateDecayProducts.size(); ++i) std::cout << actualIntermediateDecayProducts[i] <<
", ";
498 std::cout << std::endl;
499 if (requiredIntermediateDecayProducts == actualIntermediateDecayProducts)
500 std::cout <<
"*\n* These vectors match\n*\n"
503 std::cout <<
"*\n* These vectors DONT match\n*\n"
507 if (requiredIntermediateDecayProducts == actualIntermediateDecayProducts) acceptParticle =
true;
512 if (requiredIntermediateDecayProducts == actualIntermediateDecayProducts) acceptParticle =
true;
516 std::cout <<
"Checking CC state" << std::endl;
517 if (requiredIntermediateDecayProducts == actualIntermediateDecayProducts)
518 std::cout <<
"*\n* These vectors match\n*\n"
521 std::cout <<
"*\n* These vectors DONT match\n*\n"
526 else if (particle->pdg_id() == 22)
528 else if (particle->pdg_id() == 111)
533 HepMC::FourVector myFourVector = particle->momentum();
534 if (myFourVector.perp() < 0.2) trackFailedPT =
true;
535 if (
std::abs(myFourVector.eta()) > 1.1) trackFailedETA =
true;
536 acceptParticle =
true;
539 return acceptParticle;
547 for (i = 0; i <
n; i++)
557 for (
int j = i; j <
n; j++)
566 bool particleFound =
true;
567 if (!TDatabasePDG::Instance()->GetParticle(particle.c_str()))
571 std::cout <<
"The particle, " << particle <<
" is not in the TDatabasePDG particle list" << std::endl;
573 particleFound =
false;
576 return particleFound;
583 std::sort(v.begin(), v.end());
589 return TDatabasePDG::Instance()->GetParticle(name.c_str())->PdgCode();
597 return TDatabasePDG::Instance()->GetParticle(name.c_str())->Charge()/3;
612 std::cout <<
"DST node added" << std::endl;
615 std::string baseName;
624 std::string undrscr =
"_";
625 std::string nothing =
"";
626 std::map<std::string, std::string> forbiddenStrings;
627 forbiddenStrings[
"/"] = undrscr;
628 forbiddenStrings[
"("] = undrscr;
629 forbiddenStrings[
")"] = nothing;
630 forbiddenStrings[
"+"] =
"plus";
631 forbiddenStrings[
"-"] =
"minus";
632 forbiddenStrings[
"*"] =
"star";
633 forbiddenStrings[
" "] = nothing;
634 for (
auto const& [badString, goodString] : forbiddenStrings)
636 while ((pos = baseName.find(badString)) != std::string::npos) baseName.replace(pos, 1, goodString);
644 std::cout <<
m_nodeName <<
" node added" << std::endl;
657 std::cout <<
"\n---------------DecayFinder information---------------" << std::endl;
658 std::cout <<
"Module name: " <<
Name() << std::endl;
660 std::cout <<
"Number of generated decays: " <<
m_counter << std::endl;
661 std::cout <<
" Number of decays that failed pT requirement: " <<
m_nCandFail_pT << std::endl;
662 std::cout <<
" Number of decays that failed eta requirement: " <<
m_nCandFail_eta << std::endl;
663 std::cout <<
" Number of decays that failed pT and eta requirements: " <<
m_nCandFail_pT_and_eta << std::endl;
665 std::cout <<
"-----------------------------------------------------\n"
671 std::cout <<
"----------------";
672 std::cout <<
" DecayFinderNode: " <<
m_nodeName <<
" information ";
673 std::cout <<
"----------------" << std::endl;
678 for (
unsigned int i = 0; i < decay.size(); ++i)
680 std::cout <<
"Particle Barcode: " << decay[i].first <<
", Particle PDG ID: " << decay[i].second << std::endl;
683 std::cout <<
"--------------------------------------------------------------------------------------------------" << std::endl;