analysis_INTT/src/pdgidfunc.h

0001 #ifndef PDGIDFUNC_H
0002 #define PDGIDFUNC_H
0003
0004 // Direct translate from Reference: https://github.com/scikit-hep/particle/blob/master/src/particle/pdgid/functions.py
0005
0006 #include <iostream>
0007 #include <fstream>
0008 #include <vector>
0009 #include <string>
0010 #include <sstream>
0011 #include <cmath>
0012 #include <stdio.h>
0013
0014 enum class Location
0015 {
0016     Nj = 1,
0017     Nq3 = 2,
0018     Nq2 = 3,
0019     Nq1 = 4,
0020     Nl = 5,
0021     Nr = 6,
0022     N = 7,
0023     N8 = 8,
0024     N9 = 9,
0025     N10 = 10
0026 };
0027
0028 int abspid(int pdgid) {
0029     return std::abs(int(pdgid));
0030 }
0031
0032 int _digit(int pdgid, int loc)
0033 {
0034     // Provides a convenient index into the PDGID number, whose format is in base 10.
0035     // Returns the digit at position 'loc' given that the right-most digit is at position 1.
0036     std::string sid = std::to_string(abspid(pdgid));
0037     return (loc <= sid.length()) ? std::stoi(sid.substr(sid.length() - loc, 1)) : 0;
0038 }
0039
0040 int _extra_bits(int pdgid)
0041 {
0042     return abspid(pdgid) / 10000000;
0043 }
0044
0045 int _fundamental_id(int pdgid)
0046 {
0047     if (_extra_bits(pdgid) > 0)
0048     {
0049         return 0;
0050     }
0051     int aid = abspid(pdgid);
0052     if (aid <= 100)
0053     {
0054         return aid;
0055     }
0056     if (_digit(pdgid, (int) Location::Nq2) == _digit(pdgid, (int) Location::Nq1) == 0)
0057     {
0058         return aid % 10000;
0059     }
0060     return 0;
0061 }
0062
0063 bool is_quark(int pdgid)
0064 {
0065     return (abspid(pdgid) >= 1 && abspid(pdgid) <= 8) ? true:false;
0066 }
0067
0068 bool is_meson(int pdgid)
0069 {
0070     if (_extra_bits(pdgid) > 0)
0071     {
0072         return false;
0073     }
0074     int aid = abspid(pdgid);
0075     if (aid <= 100)
0076     {
0077         return false;
0078     }
0079     if (0 < _fundamental_id(pdgid) && _fundamental_id(pdgid) <= 100)
0080     {
0081         return false;
0082     }
0083     // Special IDs - K(L)0, ???, K(S)0
0084     if (aid == 130 || aid == 210 || aid == 310)
0085     {
0086         return true;
0087     }
0088     // Special IDs - B(L)0, B(sL)0, B(H)0, B(sH)0
0089     if (aid == 150 || aid == 350 || aid == 510 || aid == 530)
0090     {
0091         return true;
0092     }
0093     // Special particles - reggeon, pomeron, odderon
0094     if (pdgid == 110 || pdgid == 990 || pdgid == 9990)
0095     {
0096         return true;
0097     }
0098     // Generator-specific "particles" for GEANT tracking purposes
0099     if (aid == 998 || aid == 999)
0100     {
0101         return false;
0102     }
0103     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)
0104     {
0105         // check for illegal antiparticles
0106         return (_digit(pdgid, (int) Location::Nq3) != _digit(pdgid, (int) Location::Nq2) || pdgid >= 0);
0107     }
0108     return false;
0109 }
0110
0111 bool is_baryon(int pdgid)
0112 {
0113     int aid = abspid(pdgid);
0114     if (aid <= 100)
0115     {
0116         return false;
0117     }
0118     // Special case of proton and neutron:
0119     // needs to be checked first since _extra_bits(pdgid) > 0 for nuclei
0120     if (aid == 1000000010 || aid == 1000010010)
0121     {
0122         return true;
0123     }
0124
0125     if (_extra_bits(pdgid) > 0)
0126     {
0127         return false;
0128     }
0129
0130     if (0 < _fundamental_id(pdgid) && _fundamental_id(pdgid) <= 100)
0131     {
0132         return false;
0133     }
0134
0135     // Old codes for diffractive p and n (MC usage)
0136     if (aid == 2110 || aid == 2210)
0137     {
0138         return true;
0139     }
0140
0141     if (aid == 9221132 || aid == 9331122)
0142     {
0143         return false;
0144     }
0145
0146     return (
0147         _digit(pdgid, (int) Location::Nj) > 0
0148         && _digit(pdgid, (int) Location::Nq3) > 0
0149         && _digit(pdgid, (int) Location::Nq2) > 0
0150         && _digit(pdgid, (int) Location::Nq1) > 0
0151     );
0152 }
0153
0154 bool is_SUSY(int pdgid) {
0155     if (_extra_bits(pdgid) > 0) {
0156         return false;
0157     }
0158     if (_digit(pdgid, (int) Location::N) != 1 && _digit(pdgid, (int) Location::N) != 2)
0159     {
0160         return false;
0161     }
0162     if (_digit(pdgid, (int) Location::Nr) != 0)
0163     {
0164         return false;
0165     }
0166     return _fundamental_id(pdgid) != 0;
0167 }
0168
0169 bool is_Rhadron(int pdgid) {
0170     if (_extra_bits(pdgid) > 0) {
0171         return false;
0172     }
0173     if (_digit(pdgid, (int) Location::N) != 1) {
0174         return false;
0175     }
0176     if (_digit(pdgid, (int) Location::Nr) != 0) {
0177         return false;
0178     }
0179     if (is_SUSY(pdgid)) {
0180         return false;
0181     }
0182     // All R-hadrons have at least 3 core digits
0183     return (
0184         _digit(pdgid, (int) Location::Nq2) != 0
0185         && _digit(pdgid, (int) Location::Nq3) != 0
0186         && _digit(pdgid, (int) Location::Nj) != 0
0187     );
0188 }
0189
0190 bool is_hadron(int pdgid) {
0191     // Special case of proton and neutron:
0192     // needs to be checked first since _extra_bits(pdgid) > 0 for nuclei
0193     if (abs(pdgid) == 1000000010 || abs(pdgid) == 1000010010) {
0194         return true;
0195     }
0196     if (_extra_bits(pdgid) > 0) {
0197         return false;
0198     }
0199     if (is_meson(pdgid)) {
0200         return true;
0201     }
0202     if (is_baryon(pdgid)) {
0203         return true;
0204     }
0205     // Irrelevant test since all pentaquarks are baryons!
0206     // if (is_pentaquark(pdgid)) return true;
0207     return is_Rhadron(pdgid);
0208 }
0209
0210 int A(int pdgid) {
0211     /*
0212     Returns the atomic mass number A if the PDG ID corresponds to a nucleus, else it returns None.
0213     The (atomic) mass number is also known as the nucleon number, the total number of protons and neutrons.
0214     The number of neutrons is hence N = A - Z.
0215     */
0216     // A proton can be a Hydrogen nucleus
0217     // A neutron can be considered as a nucleus when given the PDG ID 1000000010,
0218     // hence consistency demands that A(neutron) = 1
0219     if (abspid(pdgid) == 2112 || abspid(pdgid) == 2212)
0220     {
0221         return 1;
0222     }
0223     if (_digit(pdgid, (int) Location::N10) != 1 || _digit(pdgid, (int) Location::N9) != 0)
0224     {
0225         return -1;
0226     }
0227     return (abspid(pdgid) / 10) % 1000;
0228 }
0229
0230 int Z(int pdgid) {
0231     /*
0232     Returns the nuclear charge Z if the PDG ID corresponds to a nucleus, else it returns None.
0233     The nuclear charge number Z is also known as the atomic number.
0234     For ordinary nuclei Z is equal to the number of protons.
0235     */
0236     // A proton can be a Hydrogen nucleus
0237     if (abspid(pdgid) == 2212)
0238     {
0239         return int(pdgid) / 2212;
0240     }
0241     // A neutron can be considered as a nucleus when given the PDG ID 1000000010,
0242     // hence consistency demands that Z(neutron) = 0
0243     if (abspid(pdgid) == 2112)
0244     {
0245         return 0;
0246     }
0247     if (_digit(pdgid, (int) Location::N10) != 1 || _digit(pdgid, (int) Location::N9) != 0)
0248     {
0249         return -1;
0250     }
0251     return ((abspid(pdgid) / 10000) % 1000) * (int(pdgid) / abs(int(pdgid)));
0252 }
0253
0254
0255 bool is_nucleus(int pdgid) {
0256     /*
0257     Does this PDG ID correspond to a nucleus?
0258     Ion numbers are +/- 10LZZZAAAI.
0259     AAA is A - total baryon number
0260     ZZZ is Z - total charge
0261     L is the total number of strange quarks.
0262     I is the isomer number, with I=0 corresponding to the ground state.
0263     */
0264     // A proton can be a Hydrogen nucleus
0265     // A neutron can be considered as a nucleus when given the PDG ID 1000000010,
0266     // hence consistency demands that is_nucleus(neutron) is True
0267     if (abspid(pdgid) == 2112 || abspid(pdgid) == 2212) {
0268         return true;
0269     }
0270     if (_digit(pdgid, (int) Location::N10) == 1 && _digit(pdgid, (int) Location::N9) == 0) {
0271         // Charge should always be less than or equal to the baryon number
0272         int A_pdgid = A(pdgid);
0273         int Z_pdgid = Z(pdgid);
0274
0275         // At this point neither A_pdgid nor Z_pdgid can be None,
0276         // see the definitions of the A and Z functions
0277         if (A_pdgid >= abs(Z_pdgid)) {
0278             return true;
0279         }
0280     }
0281     return false;
0282 }
0283
0284 bool is_Qball(int pdgid) {
0285     /*
0286     Does this PDG ID correspond to a Q-ball or any exotic particle with electric charge beyond the qqq scheme?
0287     Ad-hoc numbering for such particles is +/- 100XXXY0, where XXX.Y is the charge.
0288     */
0289     if (_extra_bits(pdgid) != 1) {
0290         return false;
0291     }
0292     if (_digit(pdgid, (int)  Location::N) != 0) {
0293         return false;
0294     }
0295     if (_digit(pdgid, (int)  Location::Nr) != 0) {
0296         return false;
0297     }
0298     if ((abspid(pdgid) / 10) % 10000 == 0) {
0299         return false;
0300     }
0301     return _digit(pdgid, (int) Location::Nj) == 0;
0302 }
0303
0304 bool is_dyon(int pdgid)
0305 {
0306     /*
0307     Does this PDG ID correspond to a Dyon, a magnetic monopole?
0308     Magnetic monopoles and Dyons are assumed to have one unit of Dirac monopole charge
0309     and a variable integer number xyz units of electric charge,
0310     where xyz stands for Nq1 Nq2 Nq3.
0311     Codes 411xyz0 are used when the magnetic and electrical charge sign agree and 412xyz0 when they disagree,
0312     with the overall sign of the particle set by the magnetic charge.
0313     For now, no spin information is provided.
0314     */
0315     if (_extra_bits(pdgid) > 0) {
0316         return false;
0317     }
0318     if (_digit(pdgid, (int) Location::N) != 4) {
0319         return false;
0320     }
0321     if (_digit(pdgid, (int) Location::Nr) != 1) {
0322         return false;
0323     }
0324     if (_digit(pdgid, (int) Location::Nl) != 1 && _digit(pdgid, (int) Location::Nl) != 2) {
0325         return false;
0326     }
0327     if (_digit(pdgid, (int) Location::Nq3) == 0) {
0328         return false;
0329     }
0330     return _digit(pdgid, (int) Location::Nj) == 0;
0331 }
0332
0333 bool is_diquark(int pdgid) {
0334     if (_extra_bits(pdgid) > 0) {
0335         return false;
0336     }
0337     if (abs(abspid(pdgid)) <= 100) {
0338         return false;
0339     }
0340     if (0 < _fundamental_id(pdgid) && _fundamental_id(pdgid) <= 100) {
0341         return false;
0342     }
0343     return (
0344         _digit(pdgid, (int) Location::Nj) > 0
0345         && _digit(pdgid, (int) Location::Nq3) == 0
0346         && _digit(pdgid, (int) Location::Nq2) > 0
0347         && _digit(pdgid, (int) Location::Nq1) > 0
0348     );
0349 }
0350
0351 bool is_generator_specific(int pdgid) {
0352     int aid = abs(abspid(pdgid));
0353     if (81 <= aid && aid <= 100) {
0354         return true;
0355     }
0356     if (901 <= aid && aid <= 930) {
0357         return true;
0358     }
0359     if (1901 <= aid && aid <= 1930) {
0360         return true;
0361     }
0362     if (2901 <= aid && aid <= 2930) {
0363         return true;
0364     }
0365     if (3901 <= aid && aid <= 3930) {
0366         return true;
0367     }
0368     if (aid == 998 || aid == 999) {
0369         return true;
0370     }
0371     if (aid == 20022 || aid == 480000000) {  // Special cases of opticalphoton and geantino
0372         return true;
0373     }
0374     return false;
0375 }
0376
0377 bool is_technicolor(int pdgid)
0378 {
0379     if (_extra_bits(pdgid) > 0) {
0380         return false;
0381     }
0382     return _digit(pdgid, (int) Location::N) == 3;
0383 }
0384
0385 bool is_excited_quark_or_lepton(int pdgid)
0386 {
0387     if (_extra_bits(pdgid) > 0)
0388     {
0389         return false;
0390     }
0391     if (_fundamental_id(pdgid) == 0)
0392     {
0393         return false;
0394     }
0395
0396     return _digit(pdgid, (int) Location::N) == 4 && _digit(pdgid, (int) Location::Nr) == 0;
0397 }
0398
0399 bool is_gauge_boson_or_higgs(int pdgid) {
0400     /*
0401     Does this PDG ID correspond to a gauge boson or a Higgs?
0402     Codes 21-30 are reserved for the Standard Model gauge bosons and the Higgs.
0403     The graviton and the boson content of a two-Higgs-doublet scenario
0404     and of additional SU(2)xU(1) groups are found in the range 31-40.
0405     */
0406     return (abs(pdgid) >= 21 && abs(pdgid) <= 40) ? true:false ;
0407 }
0408
0409 bool is_pentaquark(int pdgid) {
0410     if (_extra_bits(pdgid) > 0) {
0411         return false;
0412     }
0413     if (_digit(pdgid, (int) Location::N) != 9) {
0414         return false;
0415     }
0416     if (_digit(pdgid, (int) Location::Nr) == 9 || _digit(pdgid, (int) Location::Nr) == 0) {
0417         return false;
0418     }
0419     if (_digit(pdgid, (int) Location::Nj) == 9 || _digit(pdgid, (int) Location::Nl) == 0) {
0420         return false;
0421     }
0422     if (_digit(pdgid, (int) Location::Nq1) == 0) {
0423         return false;
0424     }
0425     if (_digit(pdgid, (int) Location::Nq2) == 0) {
0426         return false;
0427     }
0428     if (_digit(pdgid, (int) Location::Nq3) == 0) {
0429         return false;
0430     }
0431     if (_digit(pdgid, (int) Location::Nj) == 0) {
0432         return false;
0433     }
0434     if (_digit(pdgid, (int) Location::Nq2) > _digit(pdgid, (int) Location::Nq1)) {
0435         return false;
0436     }
0437     if (_digit(pdgid, (int) Location::Nq1) > _digit(pdgid, (int) Location::Nl)) {
0438         return false;
0439     }
0440     return _digit(pdgid, (int) Location::Nl) <= _digit(pdgid, (int) Location::Nr);
0441 }
0442
0443 bool is_valid(int pdgid) {
0444     /*Is it a valid PDG ID?*/
0445     if (is_gauge_boson_or_higgs(pdgid)) { // test first since quickest check
0446         return true;
0447     }
0448     if (_fundamental_id(pdgid) != 0) { // function always returns a number >= 0
0449         return true;
0450     }
0451     if (is_meson(pdgid)) {
0452         return true;
0453     }
0454     if (is_baryon(pdgid)) {
0455         return true;
0456     }
0457     if (is_pentaquark(pdgid)) {
0458         return true;
0459     }
0460     if (is_SUSY(pdgid)) {
0461         return true;
0462     }
0463     if (is_Rhadron(pdgid)) {
0464         return true;
0465     }
0466     if (is_dyon(pdgid)) {
0467         return true;
0468     }
0469     if (is_diquark(pdgid)) {
0470         return true;
0471     }
0472     if (is_generator_specific(pdgid)) {
0473         return true;
0474     }
0475     if (is_technicolor(pdgid)) {
0476         return true;
0477     }
0478     if (is_excited_quark_or_lepton(pdgid)) {
0479         return true;
0480     }
0481     if (_extra_bits(pdgid) > 0) {
0482         return is_Qball(pdgid) || is_nucleus(pdgid);
0483     }
0484     return false;
0485 }
0486
0487 int three_charge(int pdgid)
0488 {
0489     if (!is_valid(pdgid))
0490     {
0491         // return NULL;
0492         return -999;
0493     }
0494
0495     int aid = abspid(pdgid);
0496     // int charge = NULL; // pylint: disable=redefined-outer-name
0497     int charge = -999;
0498     int ch100[] = {
0499         -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
0500     };
0501     int q1 = _digit(pdgid, (int) Location::Nq1);
0502     int q2 = _digit(pdgid, (int) Location::Nq2);
0503     int q3 = _digit(pdgid, (int) Location::Nq3);
0504     int sid = _fundamental_id(pdgid);
0505     // std::cout << "sid=" << sid << std::endl;
0506
0507     if (_extra_bits(pdgid) > 0)
0508     {
0509         if (is_nucleus(pdgid))
0510         { // ion
0511             int Z_pdgid = Z(pdgid);
0512             // return (Z_pdgid == NULL) ? NULL : 3 * Z_pdgid;
0513             return 3 * Z_pdgid;
0514         }
0515         if (is_Qball(pdgid))
0516         { // Qball
0517             charge = 3 * ((aid / 10) % 10000);
0518         }
0519         else
0520         { // this should never be reached in the present numbering scheme
0521             // return NULL; // since extra bits exist only for Q-balls and nuclei
0522             return -999;
0523         }
0524     }
0525     else if (is_dyon(pdgid))
0526     { // Dyon
0527         charge = 3 * ((aid / 10) % 1000);
0528         // this is half right
0529         // the charge sign will be changed below if pid < 0
0530         if (_digit(pdgid, (int) Location::Nl) == 2)
0531         {
0532             charge = -charge;
0533         }
0534     }
0535     else if (sid > 0 && sid <= 100)
0536     { // use table
0537         charge = ch100[sid - 1];
0538         if (aid == 1000017 || aid == 1000018 || aid == 1000034 || aid == 1000052 || aid == 1000053 || aid == 1000054)
0539         {
0540             charge = 0;
0541         }
0542         if (aid == 5100061 || aid == 5100062)
0543         {
0544             charge = 6;
0545         }
0546     }
0547     else if (_digit(pdgid, (int) Location::Nj) == 0)
0548     { // KL, KS, or undefined
0549         return 0;
0550     }
0551     else if (q1 == 0 || (is_Rhadron(pdgid) && q1 == 9))
0552     { // mesons
0553         if (q2 == 3 || q2 == 5)
0554         {
0555             charge = ch100[q3 - 1] - ch100[q2 - 1];
0556         }
0557         else
0558         {
0559             charge = ch100[q2 - 1] - ch100[q3 - 1];
0560         }
0561     }
0562     else if (q3 == 0)
0563     { // diquarks
0564         charge = ch100[q2 - 1] + ch100[q1 - 1];
0565     }
0566     else if (is_baryon(pdgid) || (is_Rhadron(pdgid) && _digit(pdgid, (int) Location::Nl) == 9))
0567     { // baryons
0568         charge = ch100[q3 - 1] + ch100[q2 - 1] + ch100[q1 - 1];
0569     }
0570     else if (is_pentaquark(pdgid))
0571     {
0572         if (aid == 9221132)
0573         {
0574             charge = 3;
0575         }
0576         if (aid == 9331122)
0577         {
0578             charge = -6;
0579         }
0580     }
0581
0582     if (charge != -999 && (int) pdgid < 0) {
0583         charge = -charge;
0584     }
0585     return charge;
0586 }
0587
0588 float charge(int pdgid) {
0589     float three_charge_pdgid = three_charge(pdgid);
0590     // if (three_charge_pdgid == NULL)
0591     if (three_charge_pdgid == -999)
0592     {
0593         // return NULL;
0594         return -999;
0595     }
0596     if (!is_Qball(pdgid))
0597     {
0598         return three_charge_pdgid / 3.0;
0599     }
0600     return three_charge_pdgid / 30.0;
0601 }
0602
0603 bool is_chargedHadron(int pdgid)
0604 {
0605     bool is_hadron = (is_meson(pdgid) || is_baryon(pdgid)) ? true : false;
0606     bool is_charged = (fabs(charge(pdgid)) > 0) ? true : false;
0607     return (is_hadron && is_charged) ? true : false;
0608 }
0609
0610 bool is_charged(int pdgid)
0611 {
0612     bool is_charged = (fabs(charge(pdgid)) > 0) ? true : false;
0613     return (is_charged) ? true : false;
0614 }
0615
0616
0617 #endif