Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
pdgidfunc.h
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pdgidfunc.h
1 #ifndef PDGIDFUNC_H
2 #define PDGIDFUNC_H
3 
4 // Direct translate from Reference: https://github.com/scikit-hep/particle/blob/master/src/particle/pdgid/functions.py
5 
6 #include <iostream>
7 #include <fstream>
8 #include <vector>
9 #include <string>
10 #include <sstream>
11 #include <cmath>
12 #include <stdio.h>
13 
14 enum class Location
15 {
16  Nj = 1,
17  Nq3 = 2,
18  Nq2 = 3,
19  Nq1 = 4,
20  Nl = 5,
21  Nr = 6,
22  N = 7,
23  N8 = 8,
24  N9 = 9,
25  N10 = 10
26 };
27 
28 int abspid(int pdgid) {
29  return std::abs(int(pdgid));
30 }
31 
32 int _digit(int pdgid, int loc)
33 {
34  // Provides a convenient index into the PDGID number, whose format is in base 10.
35  // Returns the digit at position 'loc' given that the right-most digit is at position 1.
37  return (loc <= sid.length()) ? std::stoi(sid.substr(sid.length() - loc, 1)) : 0;
38 }
39 
40 int _extra_bits(int pdgid)
41 {
42  return abspid(pdgid) / 10000000;
43 }
44 
45 int _fundamental_id(int pdgid)
46 {
47  if (_extra_bits(pdgid) > 0)
48  {
49  return 0;
50  }
51  int aid = abspid(pdgid);
52  if (aid <= 100)
53  {
54  return aid;
55  }
56  if (_digit(pdgid, (int) Location::Nq2) == _digit(pdgid, (int) Location::Nq1) == 0)
57  {
58  return aid % 10000;
59  }
60  return 0;
61 }
62 
63 bool is_quark(int pdgid)
64 {
65  return (abspid(pdgid) >= 1 && abspid(pdgid) <= 8) ? true:false;
66 }
67 
68 bool is_meson(int pdgid)
69 {
70  if (_extra_bits(pdgid) > 0)
71  {
72  return false;
73  }
74  int aid = abspid(pdgid);
75  if (aid <= 100)
76  {
77  return false;
78  }
79  if (0 < _fundamental_id(pdgid) && _fundamental_id(pdgid) <= 100)
80  {
81  return false;
82  }
83  // Special IDs - K(L)0, ???, K(S)0
84  if (aid == 130 || aid == 210 || aid == 310)
85  {
86  return true;
87  }
88  // Special IDs - B(L)0, B(sL)0, B(H)0, B(sH)0
89  if (aid == 150 || aid == 350 || aid == 510 || aid == 530)
90  {
91  return true;
92  }
93  // Special particles - reggeon, pomeron, odderon
94  if (pdgid == 110 || pdgid == 990 || pdgid == 9990)
95  {
96  return true;
97  }
98  // Generator-specific "particles" for GEANT tracking purposes
99  if (aid == 998 || aid == 999)
100  {
101  return false;
102  }
103  if (_digit(pdgid, (int) Location::Nj) > 0 && _digit(pdgid, (int) Location::Nq3) > 0 && _digit(pdgid, (int) Location::Nq2) > 0 && _digit(pdgid, (int) Location::Nq1) == 0)
104  {
105  // check for illegal antiparticles
106  return (_digit(pdgid, (int) Location::Nq3) != _digit(pdgid, (int) Location::Nq2) || pdgid >= 0);
107  }
108  return false;
109 }
110 
111 bool is_baryon(int pdgid)
112 {
113  int aid = abspid(pdgid);
114  if (aid <= 100)
115  {
116  return false;
117  }
118  // Special case of proton and neutron:
119  // needs to be checked first since _extra_bits(pdgid) > 0 for nuclei
120  if (aid == 1000000010 || aid == 1000010010)
121  {
122  return true;
123  }
124 
125  if (_extra_bits(pdgid) > 0)
126  {
127  return false;
128  }
129 
130  if (0 < _fundamental_id(pdgid) && _fundamental_id(pdgid) <= 100)
131  {
132  return false;
133  }
134 
135  // Old codes for diffractive p and n (MC usage)
136  if (aid == 2110 || aid == 2210)
137  {
138  return true;
139  }
140 
141  if (aid == 9221132 || aid == 9331122)
142  {
143  return false;
144  }
145 
146  return (
147  _digit(pdgid, (int) Location::Nj) > 0
148  && _digit(pdgid, (int) Location::Nq3) > 0
149  && _digit(pdgid, (int) Location::Nq2) > 0
150  && _digit(pdgid, (int) Location::Nq1) > 0
151  );
152 }
153 
154 bool is_SUSY(int pdgid) {
155  if (_extra_bits(pdgid) > 0) {
156  return false;
157  }
158  if (_digit(pdgid, (int) Location::N) != 1 && _digit(pdgid, (int) Location::N) != 2)
159  {
160  return false;
161  }
162  if (_digit(pdgid, (int) Location::Nr) != 0)
163  {
164  return false;
165  }
166  return _fundamental_id(pdgid) != 0;
167 }
168 
169 bool is_Rhadron(int pdgid) {
170  if (_extra_bits(pdgid) > 0) {
171  return false;
172  }
173  if (_digit(pdgid, (int) Location::N) != 1) {
174  return false;
175  }
176  if (_digit(pdgid, (int) Location::Nr) != 0) {
177  return false;
178  }
179  if (is_SUSY(pdgid)) {
180  return false;
181  }
182  // All R-hadrons have at least 3 core digits
183  return (
184  _digit(pdgid, (int) Location::Nq2) != 0
185  && _digit(pdgid, (int) Location::Nq3) != 0
186  && _digit(pdgid, (int) Location::Nj) != 0
187  );
188 }
189 
190 bool is_hadron(int pdgid) {
191  // Special case of proton and neutron:
192  // needs to be checked first since _extra_bits(pdgid) > 0 for nuclei
193  if (abs(pdgid) == 1000000010 || abs(pdgid) == 1000010010) {
194  return true;
195  }
196  if (_extra_bits(pdgid) > 0) {
197  return false;
198  }
199  if (is_meson(pdgid)) {
200  return true;
201  }
202  if (is_baryon(pdgid)) {
203  return true;
204  }
205  // Irrelevant test since all pentaquarks are baryons!
206  // if (is_pentaquark(pdgid)) return true;
207  return is_Rhadron(pdgid);
208 }
209 
210 int A(int pdgid) {
211  /*
212  Returns the atomic mass number A if the PDG ID corresponds to a nucleus, else it returns None.
213  The (atomic) mass number is also known as the nucleon number, the total number of protons and neutrons.
214  The number of neutrons is hence N = A - Z.
215  */
216  // A proton can be a Hydrogen nucleus
217  // A neutron can be considered as a nucleus when given the PDG ID 1000000010,
218  // hence consistency demands that A(neutron) = 1
219  if (abspid(pdgid) == 2112 || abspid(pdgid) == 2212)
220  {
221  return 1;
222  }
223  if (_digit(pdgid, (int) Location::N10) != 1 || _digit(pdgid, (int) Location::N9) != 0)
224  {
225  return -1;
226  }
227  return (abspid(pdgid) / 10) % 1000;
228 }
229 
230 int Z(int pdgid) {
231  /*
232  Returns the nuclear charge Z if the PDG ID corresponds to a nucleus, else it returns None.
233  The nuclear charge number Z is also known as the atomic number.
234  For ordinary nuclei Z is equal to the number of protons.
235  */
236  // A proton can be a Hydrogen nucleus
237  if (abspid(pdgid) == 2212)
238  {
239  return int(pdgid) / 2212;
240  }
241  // A neutron can be considered as a nucleus when given the PDG ID 1000000010,
242  // hence consistency demands that Z(neutron) = 0
243  if (abspid(pdgid) == 2112)
244  {
245  return 0;
246  }
247  if (_digit(pdgid, (int) Location::N10) != 1 || _digit(pdgid, (int) Location::N9) != 0)
248  {
249  return -1;
250  }
251  return ((abspid(pdgid) / 10000) % 1000) * (int(pdgid) / abs(int(pdgid)));
252 }
253 
254 
255 bool is_nucleus(int pdgid) {
256  /*
257  Does this PDG ID correspond to a nucleus?
258  Ion numbers are +/- 10LZZZAAAI.
259  AAA is A - total baryon number
260  ZZZ is Z - total charge
261  L is the total number of strange quarks.
262  I is the isomer number, with I=0 corresponding to the ground state.
263  */
264  // A proton can be a Hydrogen nucleus
265  // A neutron can be considered as a nucleus when given the PDG ID 1000000010,
266  // hence consistency demands that is_nucleus(neutron) is True
267  if (abspid(pdgid) == 2112 || abspid(pdgid) == 2212) {
268  return true;
269  }
270  if (_digit(pdgid, (int) Location::N10) == 1 && _digit(pdgid, (int) Location::N9) == 0) {
271  // Charge should always be less than or equal to the baryon number
272  int A_pdgid = A(pdgid);
273  int Z_pdgid = Z(pdgid);
274 
275  // At this point neither A_pdgid nor Z_pdgid can be None,
276  // see the definitions of the A and Z functions
277  if (A_pdgid >= abs(Z_pdgid)) {
278  return true;
279  }
280  }
281  return false;
282 }
283 
284 bool is_Qball(int pdgid) {
285  /*
286  Does this PDG ID correspond to a Q-ball or any exotic particle with electric charge beyond the qqq scheme?
287  Ad-hoc numbering for such particles is +/- 100XXXY0, where XXX.Y is the charge.
288  */
289  if (_extra_bits(pdgid) != 1) {
290  return false;
291  }
292  if (_digit(pdgid, (int) Location::N) != 0) {
293  return false;
294  }
295  if (_digit(pdgid, (int) Location::Nr) != 0) {
296  return false;
297  }
298  if ((abspid(pdgid) / 10) % 10000 == 0) {
299  return false;
300  }
301  return _digit(pdgid, (int) Location::Nj) == 0;
302 }
303 
304 bool is_dyon(int pdgid)
305 {
306  /*
307  Does this PDG ID correspond to a Dyon, a magnetic monopole?
308  Magnetic monopoles and Dyons are assumed to have one unit of Dirac monopole charge
309  and a variable integer number xyz units of electric charge,
310  where xyz stands for Nq1 Nq2 Nq3.
311  Codes 411xyz0 are used when the magnetic and electrical charge sign agree and 412xyz0 when they disagree,
312  with the overall sign of the particle set by the magnetic charge.
313  For now, no spin information is provided.
314  */
315  if (_extra_bits(pdgid) > 0) {
316  return false;
317  }
318  if (_digit(pdgid, (int) Location::N) != 4) {
319  return false;
320  }
321  if (_digit(pdgid, (int) Location::Nr) != 1) {
322  return false;
323  }
324  if (_digit(pdgid, (int) Location::Nl) != 1 && _digit(pdgid, (int) Location::Nl) != 2) {
325  return false;
326  }
327  if (_digit(pdgid, (int) Location::Nq3) == 0) {
328  return false;
329  }
330  return _digit(pdgid, (int) Location::Nj) == 0;
331 }
332 
333 bool is_diquark(int pdgid) {
334  if (_extra_bits(pdgid) > 0) {
335  return false;
336  }
337  if (abs(abspid(pdgid)) <= 100) {
338  return false;
339  }
340  if (0 < _fundamental_id(pdgid) && _fundamental_id(pdgid) <= 100) {
341  return false;
342  }
343  return (
344  _digit(pdgid, (int) Location::Nj) > 0
345  && _digit(pdgid, (int) Location::Nq3) == 0
346  && _digit(pdgid, (int) Location::Nq2) > 0
347  && _digit(pdgid, (int) Location::Nq1) > 0
348  );
349 }
350 
351 bool is_generator_specific(int pdgid) {
352  int aid = abs(abspid(pdgid));
353  if (81 <= aid && aid <= 100) {
354  return true;
355  }
356  if (901 <= aid && aid <= 930) {
357  return true;
358  }
359  if (1901 <= aid && aid <= 1930) {
360  return true;
361  }
362  if (2901 <= aid && aid <= 2930) {
363  return true;
364  }
365  if (3901 <= aid && aid <= 3930) {
366  return true;
367  }
368  if (aid == 998 || aid == 999) {
369  return true;
370  }
371  if (aid == 20022 || aid == 480000000) { // Special cases of opticalphoton and geantino
372  return true;
373  }
374  return false;
375 }
376 
377 bool is_technicolor(int pdgid)
378 {
379  if (_extra_bits(pdgid) > 0) {
380  return false;
381  }
382  return _digit(pdgid, (int) Location::N) == 3;
383 }
384 
386 {
387  if (_extra_bits(pdgid) > 0)
388  {
389  return false;
390  }
391  if (_fundamental_id(pdgid) == 0)
392  {
393  return false;
394  }
395 
396  return _digit(pdgid, (int) Location::N) == 4 && _digit(pdgid, (int) Location::Nr) == 0;
397 }
398 
399 bool is_gauge_boson_or_higgs(int pdgid) {
400  /*
401  Does this PDG ID correspond to a gauge boson or a Higgs?
402  Codes 21-30 are reserved for the Standard Model gauge bosons and the Higgs.
403  The graviton and the boson content of a two-Higgs-doublet scenario
404  and of additional SU(2)xU(1) groups are found in the range 31-40.
405  */
406  return (abs(pdgid) >= 21 && abs(pdgid) <= 40) ? true:false ;
407 }
408 
409 bool is_pentaquark(int pdgid) {
410  if (_extra_bits(pdgid) > 0) {
411  return false;
412  }
413  if (_digit(pdgid, (int) Location::N) != 9) {
414  return false;
415  }
416  if (_digit(pdgid, (int) Location::Nr) == 9 || _digit(pdgid, (int) Location::Nr) == 0) {
417  return false;
418  }
419  if (_digit(pdgid, (int) Location::Nj) == 9 || _digit(pdgid, (int) Location::Nl) == 0) {
420  return false;
421  }
422  if (_digit(pdgid, (int) Location::Nq1) == 0) {
423  return false;
424  }
425  if (_digit(pdgid, (int) Location::Nq2) == 0) {
426  return false;
427  }
428  if (_digit(pdgid, (int) Location::Nq3) == 0) {
429  return false;
430  }
431  if (_digit(pdgid, (int) Location::Nj) == 0) {
432  return false;
433  }
434  if (_digit(pdgid, (int) Location::Nq2) > _digit(pdgid, (int) Location::Nq1)) {
435  return false;
436  }
437  if (_digit(pdgid, (int) Location::Nq1) > _digit(pdgid, (int) Location::Nl)) {
438  return false;
439  }
440  return _digit(pdgid, (int) Location::Nl) <= _digit(pdgid, (int) Location::Nr);
441 }
442 
443 bool is_valid(int pdgid) {
444  /*Is it a valid PDG ID?*/
445  if (is_gauge_boson_or_higgs(pdgid)) { // test first since quickest check
446  return true;
447  }
448  if (_fundamental_id(pdgid) != 0) { // function always returns a number >= 0
449  return true;
450  }
451  if (is_meson(pdgid)) {
452  return true;
453  }
454  if (is_baryon(pdgid)) {
455  return true;
456  }
457  if (is_pentaquark(pdgid)) {
458  return true;
459  }
460  if (is_SUSY(pdgid)) {
461  return true;
462  }
463  if (is_Rhadron(pdgid)) {
464  return true;
465  }
466  if (is_dyon(pdgid)) {
467  return true;
468  }
469  if (is_diquark(pdgid)) {
470  return true;
471  }
472  if (is_generator_specific(pdgid)) {
473  return true;
474  }
475  if (is_technicolor(pdgid)) {
476  return true;
477  }
478  if (is_excited_quark_or_lepton(pdgid)) {
479  return true;
480  }
481  if (_extra_bits(pdgid) > 0) {
482  return is_Qball(pdgid) || is_nucleus(pdgid);
483  }
484  return false;
485 }
486 
487 int three_charge(int pdgid)
488 {
489  if (!is_valid(pdgid))
490  {
491  // return NULL;
492  return -999;
493  }
494 
495  int aid = abspid(pdgid);
496  // int charge = NULL; // pylint: disable=redefined-outer-name
497  int charge = -999;
498  int ch100[] = {
499  -1, 2, -1, 2, -1, 2, -1, 2, 0, 0, -3, 0, -3, 0, -3, 0, -3, 0, 0, 0, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, 0, 0, 3, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 3, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
500  };
501  int q1 = _digit(pdgid, (int) Location::Nq1);
502  int q2 = _digit(pdgid, (int) Location::Nq2);
503  int q3 = _digit(pdgid, (int) Location::Nq3);
504  int sid = _fundamental_id(pdgid);
505  // std::cout << "sid=" << sid << std::endl;
506 
507  if (_extra_bits(pdgid) > 0)
508  {
509  if (is_nucleus(pdgid))
510  { // ion
511  int Z_pdgid = Z(pdgid);
512  // return (Z_pdgid == NULL) ? NULL : 3 * Z_pdgid;
513  return 3 * Z_pdgid;
514  }
515  if (is_Qball(pdgid))
516  { // Qball
517  charge = 3 * ((aid / 10) % 10000);
518  }
519  else
520  { // this should never be reached in the present numbering scheme
521  // return NULL; // since extra bits exist only for Q-balls and nuclei
522  return -999;
523  }
524  }
525  else if (is_dyon(pdgid))
526  { // Dyon
527  charge = 3 * ((aid / 10) % 1000);
528  // this is half right
529  // the charge sign will be changed below if pid < 0
530  if (_digit(pdgid, (int) Location::Nl) == 2)
531  {
532  charge = -charge;
533  }
534  }
535  else if (sid > 0 && sid <= 100)
536  { // use table
537  charge = ch100[sid - 1];
538  if (aid == 1000017 || aid == 1000018 || aid == 1000034 || aid == 1000052 || aid == 1000053 || aid == 1000054)
539  {
540  charge = 0;
541  }
542  if (aid == 5100061 || aid == 5100062)
543  {
544  charge = 6;
545  }
546  }
547  else if (_digit(pdgid, (int) Location::Nj) == 0)
548  { // KL, KS, or undefined
549  return 0;
550  }
551  else if (q1 == 0 || (is_Rhadron(pdgid) && q1 == 9))
552  { // mesons
553  if (q2 == 3 || q2 == 5)
554  {
555  charge = ch100[q3 - 1] - ch100[q2 - 1];
556  }
557  else
558  {
559  charge = ch100[q2 - 1] - ch100[q3 - 1];
560  }
561  }
562  else if (q3 == 0)
563  { // diquarks
564  charge = ch100[q2 - 1] + ch100[q1 - 1];
565  }
566  else if (is_baryon(pdgid) || (is_Rhadron(pdgid) && _digit(pdgid, (int) Location::Nl) == 9))
567  { // baryons
568  charge = ch100[q3 - 1] + ch100[q2 - 1] + ch100[q1 - 1];
569  }
570  else if (is_pentaquark(pdgid))
571  {
572  if (aid == 9221132)
573  {
574  charge = 3;
575  }
576  if (aid == 9331122)
577  {
578  charge = -6;
579  }
580  }
581 
582  if (charge != -999 && (int) pdgid < 0) {
583  charge = -charge;
584  }
585  return charge;
586 }
587 
588 float charge(int pdgid) {
589  float three_charge_pdgid = three_charge(pdgid);
590  // if (three_charge_pdgid == NULL)
591  if (three_charge_pdgid == -999)
592  {
593  // return NULL;
594  return -999;
595  }
596  if (!is_Qball(pdgid))
597  {
598  return three_charge_pdgid / 3.0;
599  }
600  return three_charge_pdgid / 30.0;
601 }
602 
603 bool is_chargedHadron(int pdgid)
604 {
605  bool is_hadron = (is_meson(pdgid) || is_baryon(pdgid)) ? true : false;
606  bool is_charged = (fabs(charge(pdgid)) > 0) ? true : false;
607  return (is_hadron && is_charged) ? true : false;
608 }
609 
610 
611 #endif