ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DecayFinder.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file DecayFinder.cc
1 /*
2  * Find decay topologies in HepMC
3  * Cameron Dean
4  * 04/06/2021
5  */
6 
7 //This is an array of PID's which decay faster than we can see in the detector
8 //This means if we find them in the HepMC record we need to skip over them
9 //If your decay descriptor contains one of these IDs then it will be removed from the list before starting the search
10 
11 #include "DecayFinder.h"
12 
13 #include "DecayFinderContainerBase.h" // for DecayFinderContainerBase::Iter
14 #include "DecayFinderContainer_v1.h" // for DecayFinderContainer_v1
15 
18 
20 
21 #include <phool/PHCompositeNode.h>
22 #include <phool/PHIODataNode.h> // for PHIODataNode
23 #include <phool/PHNodeIterator.h> // for PHNodeIterator
24 #include <phool/getClass.h>
25 
26 #include <HepMC/GenEvent.h>
27 #include <HepMC/GenParticle.h>
28 #include <HepMC/GenVertex.h> // for GenVertex::particle_iterator
29 #include <HepMC/IteratorRange.h>
30 #include <HepMC/SimpleVector.h>
31 
32 #include <TDatabasePDG.h>
33 
34 #include <algorithm> // for find, for_each, sort, transform
35 #include <ctype.h> // for toupper
36 #include <cstdlib> // for abs, size_t
37 #include <iostream> // for operator<<, endl, basic_ostream
38 #include <iterator> // for end, begin
39 #include <map> // for map, map<>::mapped_type, _Rb...
40 #include <memory> // for allocator_traits<>::value_type
41 
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};
56 
58  : m_save_dst(false)
59 {
60 }
61 
62 DecayFinder::DecayFinder(const std::string& /*name*/)
63  : m_save_dst(false)
64 {
65 }
66 
68 {
69  if (Verbosity() >= VERBOSITY_SOME)
70  {
71  std::cout << "DecayFinder name: " << Name() << std::endl;
72  std::cout << "Decay descriptor: " << m_decayDescriptor << std::endl;
73  }
74 
75  int canSearchDecay = parseDecayDescriptor();
76 
77  if (m_save_dst) createDecayNode(topNode);
78 
79  return canSearchDecay;
80 }
81 
83 {
84  bool decayFound = findDecay(topNode);
85 
86  if (decayFound)
87  {
88  m_counter += 1;
89  if (Verbosity() >= VERBOSITY_MORE) printNode(topNode);
90  }
91 
92  if (m_triggerOnDecay && !decayFound)
93  {
94  if (Verbosity() >= VERBOSITY_MORE) std::cout << "The decay, " << m_decayDescriptor << " was not found in this event, skipping" << std::endl;
96  }
97  else
98  {
100  }
101 }
102 
104 {
105  if (Verbosity() >= VERBOSITY_SOME) printInfo();
106 
107  return 0;
108 }
109 
111 {
112  bool ddCanBeParsed = true;
113 
114  size_t daughterLocator;
115 
116  std::string mother;
117  std::string intermediate;
118  std::string daughter;
119 
120  int mother_charge = 0;
121  std::vector<int> intermediates_charge;
122  std::vector<int> daughters_charge;
123 
124  std::string decayArrow = "->";
125  std::string chargeIndicator = "^";
126  std::string startIntermediate = "{";
127  std::string endIntermediate = "}";
128 
129  std::string manipulateDecayDescriptor = m_decayDescriptor;
130 
131  //Remove all white space before we begin
132  size_t pos;
133  while ((pos = manipulateDecayDescriptor.find(" ")) != std::string::npos) manipulateDecayDescriptor.replace(pos, 1, "");
134 
135  //Check for charge conjugate requirement
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); });
138 
139  //Remove the CC check if needed
140  if (checkForCC == "[]CC")
141  {
142  manipulateDecayDescriptor = manipulateDecayDescriptor.substr(1, manipulateDecayDescriptor.size() - 4);
143  m_getChargeConjugate = true;
144  }
145 
146  //Try and find the initial particle
147  size_t findMotherEndPoint = manipulateDecayDescriptor.find(decayArrow);
148  mother = manipulateDecayDescriptor.substr(0, findMotherEndPoint);
149  if (findParticle(mother))
150  m_mother_ID = abs(get_pdgcode(mother));
151  else
152  ddCanBeParsed = false;
153  manipulateDecayDescriptor.erase(0, findMotherEndPoint + decayArrow.length());
154 
155  //Try and find the intermediates
156  while ((pos = manipulateDecayDescriptor.find(startIntermediate)) != std::string::npos)
157  {
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));
161 
162  intermediate = intermediateDecay.substr(0, intermediateDecay.find(decayArrow));
163  if (findParticle(intermediate))
164  {
165  m_intermediates_ID.push_back(abs(get_pdgcode(intermediate)));
166  }
167  else
168  ddCanBeParsed = false;
169 
170  //Now find the daughters associated to this intermediate
171  int nDaughters = 0;
172  intermediateDecay.erase(0, intermediateDecay.find(decayArrow) + decayArrow.length());
173  while ((daughterLocator = intermediateDecay.find(chargeIndicator)) != std::string::npos)
174  {
175  daughter = intermediateDecay.substr(0, daughterLocator);
176  daughter += intermediateDecay.substr(daughterLocator + 1, 1);
177  if (findParticle(daughter))
178  {
179  m_daughters_ID.push_back(abs(get_pdgcode(daughter)));
180  daughters_charge.push_back(get_charge(daughter));
181  }
182  else
183  ddCanBeParsed = false;
184  intermediateDecay.erase(0, daughterLocator + 2);
185  ++nDaughters;
186  }
187  manipulateDecayDescriptor.erase(findIntermediateStartPoint, findIntermediateEndPoint + 1 - findIntermediateStartPoint);
188  m_nTracksFromIntermediates.push_back(nDaughters);
189  m_nTracksFromMother += 1;
190  }
191 
192  //Now find any remaining reconstructable tracks from the mother
193  while ((daughterLocator = manipulateDecayDescriptor.find(chargeIndicator)) != std::string::npos)
194  {
195  daughter = manipulateDecayDescriptor.substr(0, daughterLocator);
196  daughter += manipulateDecayDescriptor.substr(daughterLocator + 1, 1);
197  if (findParticle(daughter))
198  {
199  m_daughters_ID.push_back(abs(get_pdgcode(daughter)));
200  daughters_charge.push_back(get_charge(daughter));
201  }
202  else
203  ddCanBeParsed = false;
204  manipulateDecayDescriptor.erase(0, daughterLocator + 2);
205  m_nTracksFromMother += 1;
206  }
207 
208  unsigned int trackStart = 0;
209  unsigned int trackEnd = 0;
210  for (unsigned int i = 0; i < m_intermediates_ID.size(); ++i)
211  {
212  trackStart = trackEnd;
213  trackEnd = m_nTracksFromIntermediates[i] + trackStart;
214 
215  int vtxCharge = 0;
216 
217  for (unsigned int j = trackStart; j < trackEnd; ++j)
218  {
219  vtxCharge += daughters_charge[j];
220  }
221 
222  intermediates_charge.push_back(vtxCharge);
223  }
224 
225  for (unsigned int i = 0; i < m_daughters_ID.size(); ++i) mother_charge += daughters_charge[i];
226 
227  m_mother_ID = mother_charge == 0 ? m_mother_ID : mother_charge * m_mother_ID;
228  for (unsigned int i = 0; i < m_intermediates_ID.size(); ++i)
229  {
230  m_intermediates_ID[i] = intermediates_charge[i] == 0 ? m_intermediates_ID[i] : intermediates_charge[i] * m_intermediates_ID[i];
231  m_motherDecayProducts.push_back(m_intermediates_ID[i]);
232  }
233  for (unsigned int i = 0; i < m_daughters_ID.size(); ++i)
234  {
235  m_daughters_ID[i] = daughters_charge[i] == 0 ? m_daughters_ID[i] : daughters_charge[i] * m_daughters_ID[i];
236  if (i >= trackEnd) m_motherDecayProducts.push_back(m_daughters_ID[i]);
237  }
238 
239  if (ddCanBeParsed)
240  {
241  if (Verbosity() >= VERBOSITY_MORE) std::cout << "Your decay descriptor can be parsed" << std::endl;
242  return 0;
243  }
244  else
245  {
246  if (Verbosity() >= VERBOSITY_SOME) std::cout << "Your decay descriptor cannot be parsed, " << Name() << " will not be registered" << std::endl;
248  }
249 }
250 
252 {
253  bool decayWasFound = false;
254  bool aTrackFailedPT = false;
255  bool aTrackFailedETA = false;
256 
257  int n = sizeof(listOfResonantPIDs) / sizeof(listOfResonantPIDs[0]);
258  for (unsigned int i = 0; i < m_intermediates_ID.size(); ++i)
260 
261  m_geneventmap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
262  if (!m_geneventmap)
263  {
264  std::cout << "DecayFinder: Missing node PHHepMCGenEventMap" << std::endl;
265  return 0;
266  }
267 
268  m_genevt = m_geneventmap->get(1);
269  if (!m_genevt)
270  {
271  std::cout << "DecayFinder: Missing node PHHepMCGenEvent" << std::endl;
272  return 0;
273  }
274 
275  HepMC::GenEvent* theEvent = m_genevt->getEvent();
276 
277  std::vector<int> positive_motherDecayProducts;
278  for (unsigned int i = 0; i < m_motherDecayProducts.size(); ++i)
279  positive_motherDecayProducts.push_back(std::abs(m_motherDecayProducts[i]));
280 
281  for (HepMC::GenEvent::particle_const_iterator p = theEvent->particles_begin(); p != theEvent->particles_end(); ++p)
282  {
283  if ((*p)->pdg_id() == m_mother_ID)
284  {
285  if (Verbosity() >= VERBOSITY_MAX) std::cout << "parent->pdg_id(): " << (*p)->pdg_id() << std::endl;
286 
287  bool breakOut = false;
288  std::vector<int> correctMotherProducts;
289  decayChain.push_back(std::make_pair((*p)->barcode(), (*p)->pdg_id()));
290 
291  for (HepMC::GenVertex::particle_iterator children = (*p)->end_vertex()->particles_begin(HepMC::children);
292  children != (*p)->end_vertex()->particles_end(HepMC::children); ++children)
293  {
294  if (Verbosity() >= VERBOSITY_MAX) std::cout << "--children->pdg_id(): " << (*children)->pdg_id() << std::endl;
295 
296  if (!m_allowPhotons && (*children)->pdg_id() == 22)
297  {
298  breakOut = true;
299  break;
300  }
301  if (!m_allowPi0 && (*children)->pdg_id() == 111)
302  {
303  breakOut = true;
304  break;
305  }
306 
307  //This is one of the children we are looking for
308  if (std::find(positive_motherDecayProducts.begin(), positive_motherDecayProducts.end(),
309  std::abs((*children)->pdg_id())) != positive_motherDecayProducts.end())
310  {
311  if (Verbosity() >= VERBOSITY_MAX) std::cout << "This is a child you were looking for" << std::endl;
312  int needThisParticle = checkIfCorrectParticle((*children), aTrackFailedPT, aTrackFailedETA);
313  if (needThisParticle)
314  {
315  correctMotherProducts.push_back((*children)->pdg_id());
316  decayChain.push_back(std::make_pair((*children)->barcode(), (*children)->pdg_id()));
317  }
318  } //Now check if it's part of the other resonance list
319  else if (std::find(std::begin(listOfResonantPIDs), std::end(listOfResonantPIDs),
320  std::abs((*children)->pdg_id())) != std::end(listOfResonantPIDs))
321  {
322  if (Verbosity() >= VERBOSITY_MAX) std::cout << "This is a resonance to investigate further" << std::endl;
323  for (HepMC::GenVertex::particle_iterator grandchildren = (*children)->end_vertex()->particles_begin(HepMC::children);
324  grandchildren != (*children)->end_vertex()->particles_end(HepMC::children); ++grandchildren)
325  {
326  int needThisParticle = checkIfCorrectParticle((*grandchildren), aTrackFailedPT, aTrackFailedETA);
327  if (needThisParticle)
328  {
329  correctMotherProducts.push_back((*grandchildren)->pdg_id());
330  decayChain.push_back(std::make_pair((*grandchildren)->barcode(), (*grandchildren)->pdg_id()));
331  }
332  }
333  }
334  else
335  breakOut = true; //This particle is not in the decay descriptor, stop
336 
337  if (breakOut) break;
338  }
339 
340  if (breakOut) break;
341 
343  multiplyVectorByScalarAndSort(correctMotherProducts, +1);
344 
345  if (Verbosity() >= VERBOSITY_MAX)
346  {
347  std::cout << "Printing required mother decay products: ";
348  for (unsigned int i = 0; i < m_motherDecayProducts.size(); ++i) std::cout << m_motherDecayProducts[i] << ", ";
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;
353  if (m_motherDecayProducts == correctMotherProducts)
354  std::cout << "*\n* These vectors match\n*\n"
355  << std::endl;
356  else
357  std::cout << "*\n* These vectors DONT match\n*\n"
358  << std::endl;
359  }
360 
361  if (m_motherDecayProducts == correctMotherProducts) decayWasFound = true;
362  if (m_getChargeConjugate && !decayWasFound)
363  {
365  if (m_motherDecayProducts == correctMotherProducts) decayWasFound = true;
366 
367  if (Verbosity() >= VERBOSITY_MAX)
368  {
369  std::cout << "Checking CC state" << std::endl;
370  if (m_motherDecayProducts == correctMotherProducts)
371  std::cout << "*\n* These vectors match\n*\n"
372  << std::endl;
373  else
374  std::cout << "*\n* These vectors DONT match\n*\n"
375  << std::endl;
376  }
377  }
378  }
379  }
380 
381  if (decayWasFound)
382  {
383  if (aTrackFailedPT && !aTrackFailedETA)
384  m_nCandFail_pT += 1;
385  else if (!aTrackFailedPT && aTrackFailedETA)
386  m_nCandFail_eta += 1;
387  else if (aTrackFailedPT && aTrackFailedETA)
389  else
391 
392  if (m_save_dst) fillDecayNode(topNode, decayChain);
393  }
394 
395  if (decayChain.size() != 0) decayChain.clear();
396 
397  return decayWasFound;
398 }
399 
400 int DecayFinder::checkIfCorrectParticle(HepMC::GenParticle* particle, bool& trackFailedPT, bool& trackFailedETA)
401 {
402  bool acceptParticle = false;
403 // bool breakOut = false;
404 
405  std::vector<int> positive_intermediates_ID;
406  for (unsigned int i = 0; i < m_intermediates_ID.size(); ++i) positive_intermediates_ID.push_back(abs(m_intermediates_ID[i]));
407 
408  //Check if it is an intermediate or a final track
409  if (std::find(positive_intermediates_ID.begin(), positive_intermediates_ID.end(),
410  abs(particle->pdg_id())) != positive_intermediates_ID.end())
411  {
412  std::vector<int> requiredIntermediateDecayProducts;
413  std::vector<int> actualIntermediateDecayProducts;
414  //Which intermediate decay list to we need
415  auto it = std::find(positive_intermediates_ID.begin(), positive_intermediates_ID.end(), abs(particle->pdg_id()));
416  int index = it - m_intermediates_ID.begin();
417 
418  unsigned int trackStart = 0, trackStop = 0;
419  if (index == 0)
420  trackStop = m_nTracksFromIntermediates[0];
421  else
422  {
423  for (int i = 0; i < index; ++i) trackStart += m_nTracksFromIntermediates[i];
424  trackStop = trackStart + m_nTracksFromIntermediates[index];
425  }
426 
427  for (unsigned int i = trackStart; i < trackStop; ++i) requiredIntermediateDecayProducts.push_back(m_daughters_ID[i]);
428 
429  for (HepMC::GenVertex::particle_iterator grandchildren = particle->end_vertex()->particles_begin(HepMC::children);
430  grandchildren != particle->end_vertex()->particles_end(HepMC::children); ++grandchildren)
431  {
432  if (Verbosity() >= VERBOSITY_MAX) std::cout << "----grandchildren->pdg_id(): " << (*grandchildren)->pdg_id() << std::endl;
433 
434  if (std::find(std::begin(listOfResonantPIDs), std::end(listOfResonantPIDs),
435  std::abs((*grandchildren)->pdg_id())) != std::end(listOfResonantPIDs))
436  {
437  for (HepMC::GenVertex::particle_iterator greatgrandchildren = (*grandchildren)->end_vertex()->particles_begin(HepMC::children);
438  greatgrandchildren != (*grandchildren)->end_vertex()->particles_end(HepMC::children); ++greatgrandchildren)
439  {
440  if (Verbosity() >= VERBOSITY_MAX) std::cout << "--------greatgrandchildren->pdg_id(): " << (*greatgrandchildren)->pdg_id() << std::endl;
441 
442  if (m_allowPhotons && (*greatgrandchildren)->pdg_id() == 22)
443  continue;
444  else if (!m_allowPhotons && (*greatgrandchildren)->pdg_id() == 22)
445  {
446 // breakOut = true;
447  break;
448  }
449  else if (m_allowPi0 && (*greatgrandchildren)->pdg_id() == 111)
450  continue;
451  else if (!m_allowPi0 && (*greatgrandchildren)->pdg_id() == 111)
452  {
453 // breakOut = true;
454  break;
455  }
456  else
457  {
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;
462  }
463  }
464  }
465  else if (m_allowPhotons && (*grandchildren)->pdg_id() == 22)
466  continue;
467  else if (!m_allowPhotons && (*grandchildren)->pdg_id() == 22)
468  {
469 // breakOut = true;
470  break;
471  }
472  else if (m_allowPi0 && (*grandchildren)->pdg_id() == 111)
473  continue;
474  else if (!m_allowPi0 && (*grandchildren)->pdg_id() == 111)
475  {
476 // breakOut = true;
477  break;
478  }
479  else
480  {
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;
485  }
486  }
487 
488  multiplyVectorByScalarAndSort(requiredIntermediateDecayProducts, +1);
489  multiplyVectorByScalarAndSort(actualIntermediateDecayProducts, +1);
490 
491  if (Verbosity() >= VERBOSITY_MAX)
492  {
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"
501  << std::endl;
502  else
503  std::cout << "*\n* These vectors DONT match\n*\n"
504  << std::endl;
505  }
506 
507  if (requiredIntermediateDecayProducts == actualIntermediateDecayProducts) acceptParticle = true;
508 
509  if (m_getChargeConjugate && !acceptParticle)
510  {
511  multiplyVectorByScalarAndSort(requiredIntermediateDecayProducts, -1);
512  if (requiredIntermediateDecayProducts == actualIntermediateDecayProducts) acceptParticle = true;
513 
514  if (Verbosity() >= VERBOSITY_MAX)
515  {
516  std::cout << "Checking CC state" << std::endl;
517  if (requiredIntermediateDecayProducts == actualIntermediateDecayProducts)
518  std::cout << "*\n* These vectors match\n*\n"
519  << std::endl;
520  else
521  std::cout << "*\n* These vectors DONT match\n*\n"
522  << std::endl;
523  }
524  }
525  }
526  else if (particle->pdg_id() == 22)
527  return 0;
528  else if (particle->pdg_id() == 111)
529  return 0;
530  else
531  {
532  if (Verbosity() >= VERBOSITY_MAX) std::cout << "This is a final state track" << std::endl;
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;
537  }
538 
539  return acceptParticle;
540 }
541 
542 int DecayFinder::deleteElement(int arr[], int n, int x)
543 {
544  // https://www.geeksforgeeks.org/delete-an-element-from-array-using-two-traversals-and-one-traversal/
545  // Search x in array
546  int i;
547  for (i = 0; i < n; i++)
548  if (arr[i] == x)
549  break;
550 
551  // If x found in array
552  if (i < n)
553  {
554  // reduce size of array and move all
555  // elements on space ahead
556  n = n - 1;
557  for (int j = i; j < n; j++)
558  arr[j] = arr[j + 1];
559  }
560 
561  return n;
562 }
563 
565 {
566  bool particleFound = true;
567  if (!TDatabasePDG::Instance()->GetParticle(particle.c_str()))
568  {
569  if (Verbosity() >= VERBOSITY_SOME)
570  {
571  std::cout << "The particle, " << particle << " is not in the TDatabasePDG particle list" << std::endl;
572  }
573  particleFound = false;
574  }
575 
576  return particleFound;
577 }
578 
579 void DecayFinder::multiplyVectorByScalarAndSort(std::vector<int>& v, int k)
580 {
581  //https://slaystudy.com/c-multiply-vector-by-scalar/
582  std::transform(v.begin(), v.end(), v.begin(), [k](int& c) { return c * k; });
583  std::sort(v.begin(), v.end());
584 }
585 
587 {
588  if (findParticle(name))
589  return TDatabasePDG::Instance()->GetParticle(name.c_str())->PdgCode();
590  else
591  return 0;
592 }
593 
595 {
596  if (findParticle(name))
597  return TDatabasePDG::Instance()->GetParticle(name.c_str())->Charge()/3;
598  else
599  return -99;
600 }
601 
602 
604 {
605  PHNodeIterator iter(topNode);
606 
607  PHCompositeNode* lowerNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
608  if (!lowerNode)
609  {
610  lowerNode = new PHCompositeNode("DST");
611  topNode->addNode(lowerNode);
612  std::cout << "DST node added" << std::endl;
613  }
614 
615  std::string baseName;
616 
617  if (m_container_name.empty())
618  baseName = "decay";
619  else
620  baseName = m_container_name;
621 
622  //Cant have forward slashes in DST or else you make a subdirectory on save!!!
623  size_t pos;
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)
635  {
636  while ((pos = baseName.find(badString)) != std::string::npos) baseName.replace(pos, 1, goodString);
637  }
638 
639  m_nodeName = baseName + "_DecayMap";
640 
642  PHIODataNode<PHObject>* decayNode = new PHIODataNode<PHObject>(m_decayMap, m_nodeName.c_str(), "PHObject");
643  lowerNode->addNode(decayNode);
644  std::cout << m_nodeName << " node added" << std::endl;
645 
647 }
648 
650 {
651  m_decayMap = findNode::getClass<DecayFinderContainer_v1>(topNode, m_nodeName.c_str());
652  m_decayMap->insert(decay);
653 }
654 
656 {
657  std::cout << "\n---------------DecayFinder information---------------" << std::endl;
658  std::cout << "Module name: " << Name() << std::endl;
659  std::cout << "Decay descriptor: " << m_decayDescriptor << 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;
664  std::cout << "Number of decays that could be reconstructed: " << m_nCandReconstructable << std::endl;
665  std::cout << "-----------------------------------------------------\n"
666  << std::endl;
667 }
668 
670 {
671  std::cout << "----------------";
672  std::cout << " DecayFinderNode: " << m_nodeName << " information ";
673  std::cout << "----------------" << std::endl;
674  DecayFinderContainer_v1* map = findNode::getClass<DecayFinderContainer_v1>(topNode, m_nodeName.c_str());
675  for (DecayFinderContainer_v1::Iter iter = map->begin(); iter != map->end(); ++iter)
676  {
677  Decay decay = iter->second;
678  for (unsigned int i = 0; i < decay.size(); ++i)
679  {
680  std::cout << "Particle Barcode: " << decay[i].first << ", Particle PDG ID: " << decay[i].second << std::endl;
681  }
682  }
683  std::cout << "--------------------------------------------------------------------------------------------------" << std::endl;
684 }