Hammer  1.0.0
Helicity Amplitude Module for Matrix Element Reweighting
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Pdg.cc
Go to the documentation of this file.
1 ///
2 /// @file Pdg.cc
3 /// @brief Hammer particle data class
4 ///
5 
6 //**** This file is a part of the HAMMER library
7 //**** Copyright (C) 2016 - 2020 The HAMMER Collaboration
8 //**** HAMMER is licensed under version 3 of the GPL; see COPYING for details
9 //**** Please note the MCnet academic guidelines; see GUIDELINES for details
10 
11 // -*- C++ -*-
12 #include <cmath>
13 #include <regex>
14 
15 #include "Hammer/Tools/Pdg.hh"
16 #include "Hammer/Tools/ParticleUtils.hh"
17 #include "Hammer/Tools/Logging.hh"
18 
19 #include <iostream>
20 
21 using namespace std;
22 
23 namespace Hammer {
24 
25  PID* PID::_thePID = nullptr;
26 
27  PID::PID() {
28  }
29 
30  PID::~PID() {
31  _masses.clear();
32  _widths.clear();
33  _names.clear();
34  _brs.clear();
35  if (_thePID != nullptr) {
36  PID* pTemp = _thePID;
37  _thePID = nullptr;
38  delete pTemp;
39  }
40  }
41 
42  PdgId PID::toPdgCode(const string& name) const {
43  auto it = _names.find(name);
44  if(it == _names.end() || it->second.size() != 1) {
45  return 0;
46  }
47  else {
48  return it->second.front();
49  }
50  }
51 
52  vector<PdgId> PID::toPdgList(const string& name) const {
53  auto it = _names.find(name);
54  if(it == _names.end()) {
55  return {};
56  }
57  else {
58  return it->second;
59  }
60  }
61 
62  double PID::getMass(PdgId id) const {
63  auto it = _masses.find(abspid(id));
64  if(it != _masses.end()) {
65  return it->second;
66  }
67  else {
68  return -1.;
69  }
70  }
71 
72  double PID::getMass(const string& name) const {
73  return getMass(toPdgCode(name));
74  }
75 
76 
77  void PID::setMass(PdgId id, double value) {
78  auto it = _masses.find(abspid(id));
79  if(it != _masses.end()) {
80  it->second = value;
81  }
82  else {
83  _masses.insert(make_pair(abspid(id), value));
84  }
85  }
86 
87  void PID::setMass(const string& name, double value) {
88  setMass(toPdgCode(name), value);
89  }
90 
91  double PID::getWidth(PdgId id) const {
92  auto it = _widths.find(abspid(id));
93  if(it != _widths.end()) {
94  return it->second;
95  }
96  else {
97  return -1.;
98  }
99  }
100 
101  double PID::getWidth(const string& name) const {
102  return getWidth(toPdgCode(name));
103  }
104 
105  void PID::setWidth(PdgId id, double value) {
106  auto it = _widths.find(abspid(id));
107  if(it != _widths.end()) {
108  it->second = value;
109  }
110  else {
111  _widths.insert(make_pair(abspid(id), value));
112  }
113  }
114 
115  void PID::setWidth(const string& name, double value) {
116  setWidth(toPdgCode(name), value);
117  }
118 
119  size_t PID::getSpinMultiplicity(PdgId id) const {
120  if(abspid(id) <= 18 || id == PHOTON) {
121  return 2.;
122  }
123  if(isMeson(id) || isBaryon(id)) {
124  return digit(nj, id);
125  }
126  return 1.;
127  }
128 
129  size_t PID::getSpinMultiplicity(const string& name) const {
130  return getSpinMultiplicity(toPdgCode(name));
131  }
132 
133  size_t PID::getSpinMultiplicities(const vector<PdgId>& ids) const {
134  size_t result = 1ul;
135  for(auto elem: ids) {
136  result *= getSpinMultiplicity(elem);
137  }
138  return result;
139  }
140 
141  size_t PID::getSpinMultiplicities(const vector<string>& names) const {
142  size_t result = 1ul;
143  for(auto elem: names) {
144  result *= getSpinMultiplicity(elem);
145  }
146  return result;
147  }
148 
149  int PID::getThreeCharge(PdgId pid) const {
150  int charge = 0;
151  int ida, sid;
152  unsigned short q1, q2, q3;
153  static int ch100[100] = { -1, 2, -1, 2, -1, 2, -1, 2, 0, 0,
154  -3, 0, -3, 0, -3, 0, -3, 0, 0, 0,
155  0, 0, 0, 3, 0, 0, 0, 0, 0, 0,
156  0, 0, 0, 3, 0, 0, 3, 0, 0, 0,
157  0, -1, 0, 0, 0, 0, 0, 0, 0, 0,
158  0, 6, 3, 6, 0, 0, 0, 0, 0, 0,
159  0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
160  0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
161  0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
162  0, 0, 0, 0, 0, 0, 0, 0, 0, 0
163  };
164  q1 = digit(nq1, pid);
165  q2 = digit(nq2, pid);
166  q3 = digit(nq3, pid);
167  ida = abspid(pid);
168  sid = fundamentalID(pid);
169 
170  if( ida == 0 || extraBits(pid) > 0 ) { // ion or illegal
171  return 0;
172  }
173  else if( sid > 0 && sid <= 100 ) { // use table
174  charge = ch100[sid - 1];
175  }
176  else if( digit(nj, pid) == 0 ) { // KL, Ks, or undefined
177  return 0;
178  }
179  else if( isMeson(pid) ) { // mesons
180  if( q2 == 3 || q2 == 5 ) {
181  charge = ch100[q3 - 1] - ch100[q2 - 1];
182  }
183  else {
184  charge = ch100[q2 - 1] - ch100[q3 - 1];
185  }
186  }
187  else if( isDiQuark(pid) ) { // diquarks
188  charge = ch100[q2 - 1] + ch100[q1 - 1];
189  }
190  else if( isBaryon(pid) ) { // baryons
191  charge = ch100[q3 - 1] + ch100[q2 - 1] + ch100[q1 - 1];
192  }
193  else { // unknown
194  return 0;
195  }
196 
197  if( charge == 0 ) {
198  return 0;
199  }
200  else if( pid < 0 ) {
201  charge = -charge;
202  }
203 
204  return charge;
205  }
206 
207  int PID::getThreeCharge(const vector<PdgId>& ids) const {
208  int result = 0;
209  for(auto elem: ids) {
210  result += getThreeCharge(elem);
211  }
212  return result;
213  }
214 
215  int PID::getLeptonNumber(PdgId id) const {
216  if(abspid(id) < 11 || abspid(id) > 18) {
217  return 0;
218  }
219  else {
220  return (id > 0) ? 1 : -1;
221  }
222  }
223 
224  int PID::getLeptonNumber(const vector<PdgId>& ids) const {
225  int result = 0;
226  for(auto elem: ids) {
227  result += getLeptonNumber(elem);
228  }
229  return result;
230  }
231 
232  int PID::getBaryonNumber(PdgId id) const {
233  if(isBaryon(id)) {
234  return (id > 0) ? 1 : -1;
235  }
236  else {
237  return 0;
238  }
239  }
240 
241  int PID::getBaryonNumber(const vector<PdgId>& ids) const {
242  int result = 0;
243  for(auto elem: ids) {
244  result += getBaryonNumber(elem);
245  }
246  return result;
247  }
248 
249  tuple<int, int, int> PID::getLeptonFlavorNumber(PdgId id) const {
250  if(!isLepton(id)) {
251  return make_tuple(0,0,0);
252  }
253  else if(abspid(id)-11 < 2) {
254  return make_tuple((id > 0) ? 1 : -1,0,0);
255  }
256  else if(abspid(id)-13 < 2) {
257  return make_tuple(0, (id > 0) ? 1 : -1,0);
258  }
259  else {
260  return make_tuple(0,0, (id > 0) ? 1 : -1);
261  }
262  }
263 
264  tuple<int, int, int> PID::getLeptonFlavorNumber(const vector<PdgId>& ids) const {
265  tuple<int, int, int> result = make_tuple(0,0,0);
266  for(auto elem: ids) {
267  tuple<int, int, int> tmp = getLeptonFlavorNumber(elem);
268  get<0>(result) += get<0>(tmp);
269  get<1>(result) += get<1>(tmp);
270  get<2>(result) += get<2>(tmp);
271  }
272  return result;
273  }
274 
275  PID* PID::getPIDInstance() {
276  if(_thePID == nullptr) {
277  _thePID = new PID();
278  _thePID->init();
279  }
280  return _thePID;
281  }
282 
283  PID& PID::instance() {
284  return *PID::getPIDInstance();
285  }
286 
287  void PID::addBR(PdgId parent, const vector<PdgId>& daughters, const double br){
288  auto tmp = combineDaughters(daughters, {});
289  _brs[processID(parent, tmp)] = make_pair(br, parent);
290  }
291 
292  map<HashId, double> PID::getPartialWidths(){
293  map<HashId, double> res;
294  for(auto& elem: _brs){
295  auto it = _widths.find(abspid(elem.second.second));
296  if(it != _widths.end()){
297  res.insert({elem.first, (elem.second.first)*(it->second)});
298  }
299  }
300  return res;
301  }
302 
303  void PID::init() {
304  _masses[ELECTRON] = 5.11e-4;
305  _masses[MUON] = 0.1056583745;
306  _masses[TAU] = 1.77686;
307  _masses[NU_E] = 0.;
308  _masses[NU_MU] = 0.;
309  _masses[NU_TAU] = 0.;
310  _masses[PHOTON] = 0.;
311  _masses[PROTON] = 0.938272081;
312  _masses[NEUTRON] = 0.939565413;
313  _masses[PI0] = 0.1349770;
314  _masses[PIPLUS] = 0.13957061;
315  _masses[RHO0] = 0.77526;
316  _masses[RHOPLUS] = 0.77511;
317  _masses[K0L] = 0.497611;
318  _masses[K0S] = 0.497611;
319  _masses[KPLUS] = 0.493677;
320  _masses[D0] = 1.86483;
321  _masses[DPLUS] = 1.86959;
322  _masses[DSTAR] = 2.00685;
323  _masses[DSTARPLUS] = 2.01026;
324  _masses[DSSD0STAR] = 2.300;
325  _masses[DSSD0STARPLUS] = 2.349;
326  _masses[DSSD1STAR] = 2.427;
327  _masses[DSSD1STARPLUS] = 2.427;
328  _masses[DSSD1] = 2.421;
329  _masses[DSSD1PLUS] = 2.423;
330  _masses[DSSD2STAR] = 2.461;
331  _masses[DSSD2STARPLUS] = 2.465;
332  _masses[DSPLUS] = 1.96828;
333  _masses[DSSTARPLUS] = 2.1122;
334  _masses[DSSDS0STARPLUS] = 2.3178;
335  _masses[DSSDS1STARPLUS] = 2.4595;
336  _masses[DSSDS1PLUS] = 2.5351;
337  _masses[DSSDS2STARPLUS] = 2.5691;
338  _masses[BZERO] = 5.27963;
339  _masses[BPLUS] = 5.27932;
340  _masses[BS] = 5.36689;
341  _masses[LAMBDAB] = 5.620;
342  _masses[LAMBDACPLUS] = 2.286;
343 
344  const double hbarGevs = 6.582e-25;
345  _widths[TAU] = hbarGevs/2.903e-13;
346  _widths[RHO0] = 0.147;
347  _widths[RHOPLUS] = 0.149;
348  _widths[D0] = hbarGevs/4.101e-13;
349  _widths[DPLUS] = hbarGevs/1.040e-12;
350  _widths[DSTAR] = 2.1e-3;
351  _widths[DSTARPLUS] = 83.4e-6;
352  _widths[DSSD0STAR] = 0.274;
353  _widths[DSSD0STARPLUS] = 0.221;
354  _widths[DSSD1STAR] = 0.384;
355  _widths[DSSD1STARPLUS] = 0.384;
356  _widths[DSSD1] = 31.7e-3;
357  _widths[DSSD1PLUS] = 25e-3;
358  _widths[DSSD2STAR] = 47.5e-3;
359  _widths[DSSD2STARPLUS] = 46.7e-3;
360  _widths[DSPLUS] = hbarGevs/5.04e-13;
361  _widths[DSSTARPLUS] = 1.9e-3;
362  _widths[DSSDS0STARPLUS] = 3.8e-3;
363  _widths[DSSDS1STARPLUS] = 3.5e-3;
364  _widths[DSSDS1PLUS] = 0.92e-3;
365  _widths[DSSDS2STARPLUS] = 16.9e-3;
366  _widths[BZERO] = hbarGevs/1.52e-12;
367  _widths[BPLUS] = hbarGevs/1.638e-12;
368  _widths[BS] = hbarGevs/1.527e-12;
369  _widths[LAMBDAB] = hbarGevs/1.470e-12;
370  _widths[LAMBDACPLUS] = hbarGevs/2.e-13;
371 
372  addBR(ANTITAU, {NU_TAUBAR, NU_MU, ANTIMUON}, 0.1739);
373  addBR(ANTITAU, {NU_TAUBAR, NU_E, POSITRON}, 0.1782);
374  addBR(ANTITAU, {PIPLUS, NU_TAUBAR}, 0.1082);
375  addBR(RHO0, {PIPLUS, PIMINUS}, 1.0);
376  addBR(RHOPLUS, {PIPLUS, PI0}, 1.0);
377  addBR(DSTAR, {D0, PI0}, 0.647);
378  addBR(DSTAR, {D0, GAMMA}, 0.353);
379  addBR(DSTARMINUS, {DMINUS, PI0}, 0.307);
380  addBR(DSTARMINUS, {-D0, PIMINUS}, 0.677);
381  addBR(DSTARMINUS, {DMINUS, GAMMA}, 0.016);
382  addBR(DSSTARMINUS, {DSMINUS, PI0}, 0.058);
383  addBR(DSSTARMINUS, {DSMINUS, GAMMA}, 0.935);
384 
385  _names["Tau"] = {TAU, ANTITAU};
386  _names["Nu"] = {NU_E, NU_EBAR, NU_MU, NU_MUBAR, NU_TAU, NU_TAUBAR};
387 
388  _names["D"] = {DPLUS, DMINUS, D0, -D0};
389  _names["D*"] = {DSTAR, DSTARMINUS, DSTARPLUS, -DSTAR};
390  _names["D**0*"] = {DSSD0STAR, DSSD0STARMINUS, DSSD0STARPLUS, -DSSD0STAR};
391  _names["D**1*"] = {DSSD1STAR, DSSD1STARMINUS, DSSD1STARPLUS, -DSSD1STAR};
392  _names["D**1"] = {DSSD1, DSSD1MINUS, DSSD1PLUS, -DSSD1};
393  _names["D**2*"] = {DSSD2STAR, DSSD2STARMINUS, DSSD2STARPLUS, -DSSD2STAR};
394  _names["D**"] = {DSSD0STAR, DSSD0STARMINUS, DSSD0STARPLUS, -DSSD0STAR,
395  DSSD1STAR, DSSD1STARMINUS, DSSD1STARPLUS, -DSSD1STAR,
396  DSSD1, DSSD1MINUS, DSSD1PLUS, -DSSD1,
397  DSSD2STAR, DSSD2STARMINUS, DSSD2STARPLUS, -DSSD2STAR};
398 
399  _names["Ds"] = {DSPLUS, DSMINUS};
400  _names["Ds*"] = {DSSTARMINUS, DSSTARPLUS};
401  _names["Ds**0*"] = {DSSDS0STARMINUS, DSSDS0STARPLUS};
402  _names["Ds**1*"] = {DSSDS1STARMINUS, DSSDS1STARPLUS};
403  _names["Ds**1"] = {DSSDS1MINUS, DSSDS1PLUS};
404  _names["Ds**2*"] = {DSSDS2STARMINUS, DSSDS2STARPLUS};
405  _names["Ds**"] = {DSSDS0STARMINUS, DSSDS0STARPLUS,
406  DSSDS1STARMINUS, DSSDS1STARPLUS,
407  DSSDS1MINUS, DSSDS1PLUS,
408  DSSDS2STARMINUS, DSSDS2STARPLUS};
409 
410  _names["B"] = {BZERO, BMINUS, BPLUS, -BZERO};
411  _names["Bs"] = {BS, -BS};
412  _names["Lb"] = {LAMBDAB,-LAMBDAB};
413  _names["Lc"] = {LAMBDACPLUS,LAMBDACMINUS};
414  _names["K"] = {KPLUS, KMINUS, K0L, K0S};
415  _names["Pi"] = {PI0, PIPLUS, PIMINUS};
416  _names["Rho"] = {RHO0, RHOPLUS, RHOMINUS};
417 
418  _names["Ell"] = {MUON, ANTIMUON, ELECTRON, POSITRON};
419  _names["E"] = {POSITRON, ELECTRON};
420  _names["Mu"] = {MUON, ANTIMUON};
421  _names["Gamma"] = {GAMMA};
422  _names["Tau+"] = {ANTITAU};
423  _names["Tau-"] = {TAU};
424  _names["E+"] = {POSITRON};
425  _names["E-"] = {ELECTRON};
426  _names["Mu+"] = {ANTIMUON};
427  _names["Mu-"] = {MUON};
428 
429  _names["Rho0"] = {RHO0};
430  _names["Rho+"] = {RHOPLUS};
431  _names["Rho-"] = {RHOMINUS};
432  _names["Rho+-"] = {RHOPLUS, RHOMINUS};
433 
434  _names["K+"] = {KPLUS};
435  _names["K-"] = {KMINUS};
436  _names["K+-"] = {KPLUS, KMINUS};
437  _names["K0S"] = {K0S};
438  _names["K0L"] = {K0L};
439  _names["Kabar"] = {K0S, K0L};
440 
441  _names["D0"] = {D0};
442  _names["D+"] = {DPLUS};
443  _names["D-"] = {DMINUS};
444  _names["D0bar"] = {-D0};
445  _names["D+-"] = {DPLUS, DMINUS};
446  _names["Dabar"] = {D0, -D0};
447  _names["D*0"] = {DSTAR};
448  _names["D*+"] = {DSTARPLUS};
449  _names["D*-"] = {DSTARMINUS};
450  _names["D*0bar"] = {-DSTAR};
451  _names["D*+-"] = {DSTARPLUS, DSTARMINUS};
452  _names["D*abar"] = {DSTAR, -DSTAR};
453 
454  _names["D**0*0"] = {DSSD0STAR};
455  _names["D**0*+"] = {DSSD0STARPLUS};
456  _names["D**0*-"] = {DSSD0STARMINUS};
457  _names["D**0*0bar"] = {-DSSD0STAR};
458  _names["D**0*+-"] = {DSSD0STARPLUS, DSSD0STARMINUS};
459  _names["D**0*abar"] = {DSSD0STAR, -DSSD0STAR};
460  _names["D**1*0"] = {DSSD1STAR};
461  _names["D**1*+"] = {DSSD1STARPLUS};
462  _names["D**1*-"] = {DSSD1STARMINUS};
463  _names["D**1*0bar"] = {-DSSD1STAR};
464  _names["D**1*+-"] = {DSSD1STARPLUS, DSSD1STARMINUS};
465  _names["D**1*abar"] = {DSSD1STAR, -DSSD1STAR};
466  _names["D**10"] = {DSSD1};
467  _names["D**1+"] = {DSSD1PLUS};
468  _names["D**1-"] = {DSSD1MINUS};
469  _names["D**10bar"] = {-DSSD1};
470  _names["D**1+-"] = {DSSD1PLUS, DSSD1MINUS};
471  _names["D**1abar"] = {DSSD1, -DSSD1};
472  _names["D**2*0"] = {DSSD2STAR};
473  _names["D**2*+"] = {DSSD2STARPLUS};
474  _names["D**2*-"] = {DSSD2STARMINUS};
475  _names["D**2*0bar"] = {-DSSD2STAR};
476  _names["D**2*+-"] = {DSSD2STARPLUS, DSSD2STARMINUS};
477  _names["D**2*abar"] = {DSSD2STAR, -DSSD2STAR};
478 
479  _names["Lc+"] = {LAMBDACPLUS};
480  _names["Lc-"] = {LAMBDACMINUS};
481 
482  _names["Ds+"] = {DSPLUS};
483  _names["Ds-"] = {DSMINUS};
484  _names["Ds+-"] = {DSPLUS, DSMINUS};
485  _names["Ds*+"] = {DSSTARPLUS};
486  _names["Ds*-"] = {DSSTARMINUS};
487  _names["Ds*+-"] = {DSSTARPLUS, DSSTARMINUS};
488 
489  _names["Ds**0*+"] = {DSSDS0STARPLUS};
490  _names["Ds**0*-"] = {DSSDS0STARMINUS};
491  _names["Ds**0*+-"] = {DSSDS0STARPLUS, DSSDS0STARMINUS};
492  _names["Ds**1*+"] = {DSSDS1STARPLUS};
493  _names["Ds**1*-"] = {DSSDS1STARMINUS};
494  _names["Ds**1*+-"] = {DSSDS1STARPLUS, DSSDS1STARMINUS};
495  _names["Ds**1+"] = {DSSDS1PLUS};
496  _names["Ds**1-"] = {DSSDS1MINUS};
497  _names["Ds**1+-"] = {DSSDS1PLUS, DSSDS1MINUS};
498  _names["Ds**2*+"] = {DSSDS2STARPLUS};
499  _names["Ds**2*-"] = {DSSDS2STARMINUS};
500  _names["Ds**2*+-"] = {DSSDS2STARPLUS, DSSDS2STARMINUS};
501 
502  _names["B0"] = {BZERO};
503  _names["B+"] = {BPLUS};
504  _names["B-"] = {BMINUS};
505  _names["B0bar"] = {-BZERO};
506  _names["B+-"] = {BPLUS, BMINUS};
507  _names["Babar"] = {BZERO, -BZERO};
508  _names["Bs0"] = {BS};
509  _names["Bs0bar"] = {-BS};
510  _names["Lb0"] = {LAMBDAB};
511  _names["Lb0bar"] = {-LAMBDAB};
512  _names["Pi0"] = {PI0};
513  _names["Pi+"] = {PIPLUS};
514  _names["Pi-"] = {PIMINUS};
515  _names["Pi+-"] = {PIPLUS, PIMINUS};
516  _names["Nut"] = {NU_TAU};
517  _names["Nutbar"] = {NU_TAUBAR};
518  _names["Num"] = {NU_MU};
519  _names["Numbar"] = {NU_MUBAR};
520  _names["Nue"] = {NU_E};
521  _names["Nuebar"] = {NU_EBAR};
522  _names["W+"] = {WBOSON};
523  _names["W-"] = {-WBOSON};
524  _names["W"] = {WBOSON, -WBOSON};
525  }
526 
527  unsigned short PID::digit( location loc, PdgId pid ) const {
528  // PID digits (base 10) are: n nr nl nq1 nq2 nq3 nj
529  // the location enum provides a convenient index into the PID
530  int numerator = static_cast<int>(pow(10.0, (loc - 1)));
531  return static_cast<unsigned short>((abspid(pid) / numerator) % 10);
532  }
533 
534  PdgId PID::abspid( PdgId pid ) const {
535  return (pid < 0) ? -pid : pid;
536  }
537 
538  PdgId PID::extraBits( PdgId pid ) const {
539  return abspid(pid) / 10000000;
540  }
541 
542  PdgId PID::fundamentalID( PdgId pid ) const {
543  if( extraBits(pid) > 0 ) {
544  return 0;
545  }
546 
547  if( digit(nq2, pid) == 0 && digit(nq1, pid) == 0) {
548  return abspid(pid) % 10000;
549  }
550  else if( abspid(pid) <= 100 ) {
551  return abspid(pid);
552  }
553  else {
554  return 0;
555  }
556  }
557 
558  bool PID::isMeson( PdgId pid ) const {
559  if( extraBits(pid) > 0 ) {
560  return false;
561  }
562 
563  if( abspid(pid) <= 100 ) {
564  return false;
565  }
566 
567  if( fundamentalID(pid) <= 100 && fundamentalID(pid) > 0 ) {
568  return false;
569  }
570 
571  int aid = abspid(pid);
572 
573  if( aid == 130 || aid == 310 || aid == 111 || aid == 113) {
574  if (pid < 0) {
575  return false;
576  } else {
577  return true;
578  }
579  }
580 
581  // EvtGen uses some odd numbers
582  if( aid == 150 || aid == 350 || aid == 510 || aid == 530 ) {
583  return true;
584  }
585 
586  // pomeron, etc.
587  if( pid == 110 || pid == 990 || pid == 9990 ) {
588  return true;
589  }
590 
591  if( digit(nj, pid) > 0 && digit(nq3, pid) > 0
592  && digit(nq2, pid) > 0 && digit(nq1, pid) == 0 ) {
593  // check for illegal antiparticles
594  if( digit(nq3, pid) == digit(nq2, pid) && pid < 0 ) {
595  return false;
596  }
597  else {
598  return true;
599  }
600  }
601 
602  return false;
603  }
604 
605  // check to see if this is a valid baryon
606  bool PID::isBaryon( PdgId pid ) const {
607  if( extraBits(pid) > 0 ) {
608  return false;
609  }
610 
611  if( abspid(pid) <= 100 ) {
612  return false;
613  }
614 
615  if( fundamentalID(pid) <= 100 && fundamentalID(pid) > 0 ) {
616  return false;
617  }
618 
619  if( abspid(pid) == 2110 || abspid(pid) == 2210 ) {
620  return true;
621  }
622 
623  if( digit(nj, pid) > 0 && digit(nq3, pid) > 0
624  && digit(nq2, pid) > 0 && digit(nq1, pid) > 0 ) {
625  return true;
626  }
627 
628  return false;
629  }
630 
631  // check to see if this is a valid diquark
632  bool PID::isDiQuark( PdgId pid ) const {
633  if( extraBits(pid) > 0 ) {
634  return false;
635  }
636 
637  if( abspid(pid) <= 100 ) {
638  return false;
639  }
640 
641  if( fundamentalID(pid) <= 100 && fundamentalID(pid) > 0 ) {
642  return false;
643  }
644 
645  if( digit(nj, pid) > 0 && digit(nq3, pid) == 0
646  && digit(nq2, pid) > 0 && digit(nq1, pid) > 0 ) { // diquark signature
647  // EvtGen uses the diquarks for quark pairs, so, for instance,
648  // 5501 is a valid "diquark" for EvtGen
649  //if( digit(nj) == 1 && digit(nq2) == digit(nq1) ) { // illegal
650  // return false;
651  //} else {
652  return true;
653  //}
654  }
655 
656  return false;
657  }
658 
659  // is this a valid hadron ID?
660  bool PID::isHadron( PdgId pid ) const {
661  if( extraBits(pid) > 0 ) {
662  return false;
663  }
664 
665  if( isMeson(pid) ) {
666  return true;
667  }
668 
669  if( isBaryon(pid) ) {
670  return true;
671  }
672 
673  if( isPentaquark(pid) ) {
674  return true;
675  }
676 
677  return false;
678  }
679  // is this a valid lepton ID?
680  bool PID::isLepton( PdgId pid ) const {
681  if( extraBits(pid) > 0 ) {
682  return false;
683  }
684 
685  if( fundamentalID(pid) >= 11 && fundamentalID(pid) <= 18 ) {
686  return true;
687  }
688 
689  return false;
690  }
691  // is this a neutrino ID?
692  bool PID::isNeutrino( PdgId pid ) const {
693  if( extraBits(pid) > 0 ) {
694  return false;
695  }
696 
697  if( fundamentalID(pid) == 12 || fundamentalID(pid) == 14 || fundamentalID(pid) == 16 ) {
698  return true;
699  }
700 
701  return false;
702  }
703 
704  // check to see if this is a valid pentaquark
705  bool PID::isPentaquark( PdgId pid ) const {
706  // a pentaquark is of the form 9abcdej,
707  // where j is the spin and a, b, c, d, and e are quarks
708  if( extraBits(pid) > 0 ) {
709  return false;
710  }
711 
712  if( digit(n, pid) != 9 ) {
713  return false;
714  }
715 
716  if( digit(nr, pid) == 9 || digit(nr, pid) == 0 ) {
717  return false;
718  }
719 
720  if( digit(nj, pid) == 9 || digit(nl, pid) == 0 ) {
721  return false;
722  }
723 
724  if( digit(nq1, pid) == 0 ) {
725  return false;
726  }
727 
728  if( digit(nq2, pid) == 0 ) {
729  return false;
730  }
731 
732  if( digit(nq3, pid) == 0 ) {
733  return false;
734  }
735 
736  if( digit(nj, pid) == 0 ) {
737  return false;
738  }
739 
740  // check ordering
741  if( digit(nq2, pid) > digit(nq1, pid) ) {
742  return false;
743  }
744 
745  if( digit(nq1, pid) > digit(nl, pid) ) {
746  return false;
747  }
748 
749  if( digit(nl, pid) > digit(nr, pid) ) {
750  return false;
751  }
752 
753  return true;
754  }
755 
756  vector<pair<PdgId, vector<PdgId>>> PID::expandToValidVertices(const string& name) const {
757  regex regexpression ("[A-Z][^A-Z]*");
758  regex_iterator<string::const_iterator> rit (name.begin(), name.end(), regexpression);
759  regex_iterator<string::const_iterator> rend;
760  vector<string> names{};
761  while(rit!=rend) {
762  names.push_back(rit->str());
763  ++rit;
764  }
765  vector<pair<PdgId, vector<PdgId>>> results;
766  vector<vector<PdgId>> tempRes;
767  for(auto elem: names) {
768  auto res = toPdgList(elem);
769  if(res.size() == 0) {
770  return {};
771  }
772  if(tempRes.size() == 0) {
773  for(auto elem2: res) {
774  tempRes.push_back({elem2});
775  }
776  }
777  else {
778  vector<vector<PdgId>> tmpMat(tempRes.size() * res.size());
779  for(size_t i = 0; i< tempRes.size(); ++i) {
780  for(size_t j = 0; j< res.size(); ++j) {
781  tmpMat[i * res.size() + j] = tempRes[i];
782  tmpMat[i * res.size() + j].push_back(res[j]);
783  }
784  }
785  tempRes.swap(tmpMat);
786  }
787  }
788  for(auto& elemres: tempRes) {
789  PdgId parent = elemres[0];
790  if(getLeptonNumber(elemres)-2*getLeptonNumber(parent) != 0) {
791  continue;
792  }
793  if(getBaryonNumber(elemres)-2*getBaryonNumber(parent) != 0) {
794  continue;
795  }
796  if(getThreeCharge(elemres)-2*getThreeCharge(parent) != 0) {
797  continue;
798  }
799  auto left = getLeptonFlavorNumber(elemres);
800  auto right = getLeptonFlavorNumber(parent);
801  if(get<0>(left) -2*get<0>(right) !=0 ||
802  get<1>(left) -2*get<1>(right) !=0 ||
803  get<2>(left) -2*get<2>(right) !=0) {
804  continue;
805  }
806  elemres.erase(elemres.begin());
807  auto daughters = combineDaughters(elemres);
808  bool found = false;
809  for(auto elem: results) {
810  if(elem.first != parent || elem.second.size() != daughters.size()) {
811  continue;
812  }
813  bool equal = true;
814  for(size_t j = 0; j < daughters.size(); ++j) {
815  if(elem.second[j] != daughters[j]) {
816  equal = false;
817  break;
818  }
819  }
820  if(equal) {
821  found = true;
822  break;
823  }
824  }
825  if(!found) {
826  results.push_back(make_pair(parent, daughters));
827 // results.push_back(make_pair(flipSign(parent), flipSigns(daughters)));
828  }
829  }
830  return results;
831  }
832 
833  vector<HashId> PID::expandToValidVertexUIDs(const string& name, const bool& hadonly) const {
834  auto tmp = expandToValidVertices(name);
835  vector<HashId> result;
836  for(auto& elem: tmp) {
837  result.push_back(processID(elem.first,elem.second));
838  }
839  if (hadonly) {
840  auto tmpW = expandToValidVertices(name + "W");
841  auto tmpN = expandToValidVertices(name + "Nu");
842  tmpW.insert(tmpW.end(), tmpN.begin(), tmpN.end());
843  for(auto& elem: tmpW) {
844  vector<PdgId> hadDaughters = elem.second;
845  hadDaughters.erase(remove_if(hadDaughters.begin(), hadDaughters.end(), [](PdgId idx){ return (abs(idx) < 100); }),
846  hadDaughters.end());
847 // elem.second.erase(remove(elem.second.begin(), elem.second.end(), 24), elem.second.end());
848 // elem.second.erase(remove(elem.second.begin(), elem.second.end(), -24), elem.second.end());
849  result.push_back(processID(elem.first,hadDaughters));
850  }
851  }
852  return result;
853  }
854 
855  const PdgId PID::ELECTRON;
856  const PdgId PID::POSITRON;
857  const PdgId PID::EMINUS;
858  const PdgId PID::EPLUS;
859  const PdgId PID::MUON;
860  const PdgId PID::ANTIMUON;
861  const PdgId PID::TAU;
862  const PdgId PID::ANTITAU;
863  const PdgId PID::NU_E;
864  const PdgId PID::NU_EBAR;
865  const PdgId PID::NU_MU;
866  const PdgId PID::NU_MUBAR;
867  const PdgId PID::NU_TAU;
868  const PdgId PID::NU_TAUBAR;
869  const PdgId PID::PHOTON;
870  const PdgId PID::GAMMA;
871  const PdgId PID::PROTON;
872  const PdgId PID::WBOSON;
873  const PdgId PID::ANTIPROTON;
874  const PdgId PID::PBAR;
875  const PdgId PID::NEUTRON;
876  const PdgId PID::ANTINEUTRON;
877  const PdgId PID::PI0;
878  const PdgId PID::PIPLUS;
879  const PdgId PID::PIMINUS;
880  const PdgId PID::RHO0;
881  const PdgId PID::RHOPLUS;
882  const PdgId PID::RHOMINUS;
883  const PdgId PID::K0L;
884  const PdgId PID::K0S;
885  const PdgId PID::KPLUS;
886  const PdgId PID::KMINUS;
887  const PdgId PID::ETA;
888  const PdgId PID::ETAPRIME;
889  const PdgId PID::PHI;
890  const PdgId PID::OMEGA;
891  const PdgId PID::ETAC;
892  const PdgId PID::JPSI;
893  const PdgId PID::PSI2S;
894  const PdgId PID::D0;
895  const PdgId PID::DPLUS;
896  const PdgId PID::DMINUS;
897  const PdgId PID::DSTAR;
898  const PdgId PID::DSTARPLUS;
899  const PdgId PID::DSTARMINUS;
900  const PdgId PID::DSSD0STAR;
901  const PdgId PID::DSSD0STARPLUS;
902  const PdgId PID::DSSD0STARMINUS;
903  const PdgId PID::DSSD1STAR;
904  const PdgId PID::DSSD1STARPLUS;
905  const PdgId PID::DSSD1STARMINUS;
906  const PdgId PID::DSSD1;
907  const PdgId PID::DSSD1PLUS;
908  const PdgId PID::DSSD1MINUS;
909  const PdgId PID::DSSD2STAR;
910  const PdgId PID::DSSD2STARPLUS;
911  const PdgId PID::DSSD2STARMINUS;
912  const PdgId PID::DSPLUS;
913  const PdgId PID::DSMINUS;
914  const PdgId PID::DSSTARPLUS;
915  const PdgId PID::DSSTARMINUS;
916  const PdgId PID::DSSDS0STARPLUS;
917  const PdgId PID::DSSDS0STARMINUS;
918  const PdgId PID::DSSDS1STARPLUS;
919  const PdgId PID::DSSDS1STARMINUS;
920  const PdgId PID::DSSDS1PLUS;
921  const PdgId PID::DSSDS1MINUS;
922  const PdgId PID::DSSDS2STARPLUS;
923  const PdgId PID::DSSDS2STARMINUS;
924  const PdgId PID::ETAB;
925  const PdgId PID::UPSILON1S;
926  const PdgId PID::UPSILON2S;
927  const PdgId PID::UPSILON3S;
928  const PdgId PID::UPSILON4S;
929  const PdgId PID::BZERO;
930  const PdgId PID::BPLUS;
931  const PdgId PID::BMINUS;
932  const PdgId PID::BS;
933  const PdgId PID::BCPLUS;
934  const PdgId PID::BCMINUS;
935  const PdgId PID::LAMBDA;
936  const PdgId PID::SIGMA0;
937  const PdgId PID::SIGMAPLUS;
938  const PdgId PID::SIGMAMINUS;
939  const PdgId PID::LAMBDACPLUS;
940  const PdgId PID::LAMBDACMINUS;
941  const PdgId PID::LAMBDAB;
942  const PdgId PID::XI0;
943  const PdgId PID::XIMINUS;
944  const PdgId PID::XIPLUS;
945  const PdgId PID::OMEGAMINUS;
946  const PdgId PID::OMEGAPLUS;
947 
948 } // namespace Hammer
ParticleUtils.hh
PDG codes to UID functions.
Hammer::combineDaughters
std::vector< PdgId > combineDaughters(const std::vector< PdgId > &daughters, const std::vector< PdgId > &subDaughters)
combine list of codes of daughters and grandaughters (for processes which parameterise two subsequent...
Definition: ParticleUtils.cc:44
Logging.hh
Message logging routines.
Hammer::PID
Hammer class for dealing with particle data.
Definition: Pdg.hh:32
Pdg.hh
Hammer particle data class.
Hammer::processID
HashId processID(PdgId parent, const std::vector< PdgId > &allDaughters)
compute a unique ID for a given process based on the PDG codes of the parent particle and the ordered...
Definition: ParticleUtils.cc:55
Hammer::PdgId
int PdgId
Definition: Pdg.fhh:17