Line data Source code
1 : // ----------------------------------------------------------------------
2 : //
3 : // ParticleIDMethods.cc
4 : //
5 : // ----------------------------------------------------------------------
6 :
7 : #include <cmath> // for pow()
8 :
9 : #include "HepPID/ParticleIDMethods.hh"
10 : #include "HepPID/ParticleName.hh"
11 :
12 : namespace HepPID {
13 :
14 : namespace {
15 :
16 : // internal function used by hasXXX methods
17 0 : bool findQ( const int & pid, const int & q )
18 : {
19 0 : if( isDyon(pid) ) { return false; }
20 0 : if( isRhadron(pid) ) {
21 : int iz = 7;
22 0 : for( int i=6; i > 1; --i ) {
23 0 : if( digit(location(i),pid) == 0 ) {
24 : iz = i;
25 0 : } else if ( i == iz-1 ) {
26 : // ignore squark or gluino
27 : } else {
28 0 : if( digit(location(i),pid) == q ) { return true; }
29 : }
30 : }
31 : return false;
32 : }
33 0 : if( digit(nq3,pid) == q || digit(nq2,pid) == q || digit(nq1,pid) == q ) { return true; }
34 0 : if( isPentaquark(pid) ) {
35 0 : if( digit(nl,pid) == q || digit(nr,pid) == q ) { return true; }
36 : }
37 : return false;
38 : }
39 :
40 : }
41 :
42 : // absolute value
43 300656416 : int abspid( const int & pid )
44 : {
45 300656416 : return (pid < 0) ? -pid : pid;
46 : }
47 :
48 : // returns everything beyond the 7th digit (e.g. outside the numbering scheme)
49 50372267 : int extraBits( const int & pid )
50 : {
51 50372267 : return abspid(pid)/10000000;
52 : }
53 :
54 : // split the PID into constituent integers
55 166363547 : unsigned short digit( location loc, const int & pid )
56 : {
57 : // PID digits (base 10) are: n nr nl nq1 nq2 nq3 nj
58 : // the location enum provides a convenient index into the PID
59 : //
60 : // Modified for CRPropa: use precalculated values isntead of pow for
61 : // performance
62 : static unsigned int p10[] = { 1, 10, 100, 1000, 10000, 100000, 1000000,
63 : 10000000, 100000000, 1000000000};
64 166363547 : return (abspid(pid)/ p10[loc-1])%10;
65 : // int numerator = (int) std::pow(10.0,(loc-1));
66 : // return (abspid(pid)/numerator)%10;
67 : }
68 :
69 : // return the first two digits if this is a "fundamental" particle
70 : // ID = 100 is a special case (internal generator ID's are 81-100)
71 20191854 : int fundamentalID( const int & pid )
72 : {
73 20191854 : if( extraBits(pid) > 0 ) return 0;
74 20191855 : if( digit(nq2,pid) == 0 && digit(nq1,pid) == 0) {
75 20191856 : return abspid(pid)%10000;
76 2 : } else if( abspid(pid) <= 100 ) {
77 0 : return abspid(pid);
78 : } else {
79 : return 0;
80 : }
81 : }
82 :
83 : // Ion numbers are +/- 10LZZZAAAI.
84 299901 : int Z( const int & pid )
85 : {
86 : // a proton can also be a Hydrogen nucleus
87 299901 : if( abspid(pid) == 2212 ) { return 1; }
88 299901 : if( isNucleus(pid) ) return (abspid(pid)/10000)%1000;
89 : return 0;
90 : }
91 :
92 : // Ion numbers are +/- 10LZZZAAAI.
93 152688 : int A( const int & pid )
94 : {
95 : // a proton can also be a Hydrogen nucleus
96 152688 : if( abspid(pid) == 2212 ) { return 1; }
97 152688 : if( isNucleus(pid) ) return (abspid(pid)/10)%1000;
98 : return 0;
99 : }
100 :
101 : // if this is a nucleus (ion), get nLambda
102 : // Ion numbers are +/- 10LZZZAAAI.
103 0 : int lambda( const int & pid )
104 : {
105 : // a proton can also be a Hydrogen nucleus
106 0 : if( abspid(pid) == 2212 ) { return 0; }
107 0 : if( isNucleus(pid) ) return digit(n8,pid);
108 : return 0;
109 : }
110 :
111 :
112 : // --- boolean methods:
113 : //
114 :
115 : // check to see if this is a valid PID
116 0 : bool isValid( const int & pid )
117 : {
118 0 : if( extraBits(pid) > 0 ) {
119 0 : if( isNucleus(pid) ) { return true; }
120 0 : if( isQBall(pid) ) { return true; }
121 : return false;
122 : }
123 0 : if( isSUSY(pid) ) { return true; }
124 0 : if( isRhadron(pid) ) { return true; }
125 0 : if( isDyon(pid) ) { return true; }
126 : // Meson signature
127 0 : if( isMeson(pid) ) { return true; }
128 : // Baryon signature
129 0 : if( isBaryon(pid) ) { return true; }
130 : // DiQuark signature
131 0 : if( isDiQuark(pid) ) { return true; }
132 : // fundamental particle
133 0 : if( fundamentalID(pid) > 0 ) {
134 0 : if(pid > 0 ) {
135 : return true;
136 : } else {
137 0 : if( hasFundamentalAnti(pid) ) { return true; }
138 : return false;
139 : }
140 : }
141 : // pentaquark
142 0 : if( isPentaquark(pid) ) { return true; }
143 : // don't recognize this number
144 : return false;
145 : }
146 :
147 : // if this is a fundamental particle, does it have a valid antiparticle?
148 0 : bool hasFundamentalAnti( const int & pid )
149 : {
150 : // these are defined by the generator and therefore are always valid
151 0 : if( fundamentalID(pid) <= 100 && fundamentalID(pid) >= 80 ) { return true; }
152 : // check id's from 1 to 79
153 0 : if( fundamentalID(pid) > 0 && fundamentalID(pid) < 80 ) {
154 0 : if( validParticleName(-pid) ) { return true; }
155 : }
156 : return false;
157 : }
158 :
159 : // check to see if this is a valid meson
160 0 : bool isMeson( const int & pid )
161 : {
162 0 : if( extraBits(pid) > 0 ) { return false; }
163 0 : if( abspid(pid) <= 100 ) { return false; }
164 0 : if( fundamentalID(pid) <= 100 && fundamentalID(pid) > 0 ) { return false; }
165 0 : if( isRhadron(pid) ) { return false; }
166 0 : int aid = abspid(pid);
167 0 : if( aid == 130 || aid == 310 || aid == 210 ) { return true; }
168 : // EvtGen uses some odd numbers
169 0 : if( aid == 150 || aid == 350 || aid == 510 || aid == 530 ) { return true; }
170 : // pomeron, etc.
171 0 : if( pid == 110 || pid == 990 || pid == 9990 ) { return true; }
172 0 : if( digit(nj,pid) > 0 && digit(nq3,pid) > 0
173 0 : && digit(nq2,pid) > 0 && digit(nq1,pid) == 0 ) {
174 : // check for illegal antiparticles
175 0 : if( digit(nq3,pid) == digit(nq2,pid) && pid < 0 ) {
176 : return false;
177 : } else {
178 : return true;
179 : }
180 : }
181 : return false;
182 : }
183 :
184 : // check to see if this is a valid baryon
185 0 : bool isBaryon( const int & pid )
186 : {
187 0 : if( extraBits(pid) > 0 ) { return false; }
188 0 : if( abspid(pid) <= 100 ) { return false; }
189 0 : if( fundamentalID(pid) <= 100 && fundamentalID(pid) > 0 ) { return false; }
190 0 : if( isRhadron(pid) ) { return false; }
191 0 : if( isPentaquark(pid) ) { return false; }
192 0 : if( abspid(pid) == 2110 || abspid(pid) == 2210 ) { return true; }
193 0 : if( digit(nj,pid) > 0 && digit(nq3,pid) > 0
194 0 : && digit(nq2,pid) > 0 && digit(nq1,pid) > 0 ) { return true; }
195 : return false;
196 : }
197 :
198 : // check to see if this is a valid diquark
199 0 : bool isDiQuark( const int & pid )
200 : {
201 0 : if( extraBits(pid) > 0 ) { return false; }
202 0 : if( abspid(pid) <= 100 ) { return false; }
203 0 : if( fundamentalID(pid) <= 100 && fundamentalID(pid) > 0 ) { return false; }
204 0 : if( digit(nj,pid) > 0 && digit(nq3,pid) == 0
205 0 : && digit(nq2,pid) > 0 && digit(nq1,pid) > 0 ) { // diquark signature
206 : // EvtGen uses the diquarks for quark pairs, so, for instance,
207 : // 5501 is a valid "diquark" for EvtGen
208 : //if( digit(nj) == 1 && digit(nq2) == digit(nq1) ) { // illegal
209 : // return false;
210 : //} else {
211 : return true;
212 : //}
213 : }
214 : return false;
215 : }
216 :
217 : // is this a valid hadron ID?
218 0 : bool isHadron( const int & pid )
219 : {
220 0 : if( extraBits(pid) > 0 ) { return false; }
221 0 : if( isMeson(pid) ) { return true; }
222 0 : if( isBaryon(pid) ) { return true; }
223 0 : if( isPentaquark(pid) ) { return true; }
224 0 : if( isRhadron(pid) ) { return true; }
225 : return false;
226 : }
227 : // is this a valid lepton ID?
228 0 : bool isLepton( const int & pid )
229 : {
230 0 : if( extraBits(pid) > 0 ) { return false; }
231 0 : if( fundamentalID(pid) >= 11 && fundamentalID(pid) <= 18 ) { return true; }
232 : return false;
233 : }
234 :
235 : //
236 : // This implements the 2006 Monte Carlo nuclear code scheme.
237 : // Ion numbers are +/- 10LZZZAAAI.
238 : // AAA is A - total baryon number
239 : // ZZZ is Z - total charge
240 : // L is the total number of strange quarks.
241 : // I is the isomer number, with I=0 corresponding to the ground state.
242 41133615 : bool isNucleus( const int & pid )
243 : {
244 : // a proton can also be a Hydrogen nucleus
245 41133615 : if( abspid(pid) == 2212 ) { return true; }
246 : // new standard: +/- 10LZZZAAAI
247 41133604 : if( ( digit(n10,pid) == 1 ) && ( digit(n9,pid) == 0 ) ) {
248 : // charge should always be less than or equal to baryon number
249 : // the following line is A >= Z
250 749979 : if( (abspid(pid)/10)%1000 >= (abspid(pid)/10000)%1000 ) { return true; }
251 : }
252 : return false;
253 : }
254 :
255 : // check to see if this is a valid pentaquark
256 0 : bool isPentaquark( const int & pid )
257 : {
258 : // a pentaquark is of the form 9abcdej,
259 : // where j is the spin and a, b, c, d, and e are quarks
260 0 : if( extraBits(pid) > 0 ) { return false; }
261 0 : if( digit(n,pid) != 9 ) { return false; }
262 0 : if( digit(nr,pid) == 9 || digit(nr,pid) == 0 ) { return false; }
263 0 : if( digit(nj,pid) == 9 || digit(nl,pid) == 0 ) { return false; }
264 0 : if( digit(nq1,pid) == 0 ) { return false; }
265 0 : if( digit(nq2,pid) == 0 ) { return false; }
266 0 : if( digit(nq3,pid) == 0 ) { return false; }
267 0 : if( digit(nj,pid) == 0 ) { return false; }
268 : // check ordering
269 0 : if( digit(nq2,pid) > digit(nq1,pid) ) { return false; }
270 0 : if( digit(nq1,pid) > digit(nl,pid) ) { return false; }
271 0 : if( digit(nl,pid) > digit(nr,pid) ) { return false; }
272 : return true;
273 : }
274 :
275 : // is this a SUSY?
276 0 : bool isSUSY( const int & pid )
277 : {
278 : // fundamental SUSY particles have n = 1 or 2
279 0 : if( extraBits(pid) > 0 ) { return false; }
280 0 : if( digit(n,pid) != 1 && digit(n,pid) != 2 ) { return false; }
281 0 : if( digit(nr,pid) != 0 ) { return false; }
282 : // check fundamental part
283 0 : if( fundamentalID(pid) == 0 ) { return false; }
284 : return true;
285 : }
286 :
287 : // is this an R-hadron?
288 0 : bool isRhadron( const int & pid )
289 : {
290 : // an R-hadron is of the form 10abcdj, 100abcj, or 1000abj
291 : // where j is the spin, b, c, and d are quarks or gluons,
292 : // and a (the digit following the zero's) is a SUSY particle
293 0 : if( extraBits(pid) > 0 ) { return false; }
294 0 : if( digit(n,pid) != 1 ) { return false; }
295 0 : if( digit(nr,pid) != 0 ) { return false; }
296 : // make sure this isn't a SUSY particle
297 0 : if( isSUSY(pid) ) { return false; }
298 : // All R-hadrons have at least 3 core digits
299 0 : if( digit(nq2,pid) == 0 ) { return false; }
300 0 : if( digit(nq3,pid) == 0 ) { return false; }
301 0 : if( digit(nj,pid) == 0 ) { return false; }
302 : return true;
303 : }
304 :
305 : // is this a Dyon (magnetic monopole)?
306 3329523 : bool isDyon( const int & pid )
307 : {
308 : ///Magnetic monopoles and Dyons are assumed to have one unit of
309 : ///Dirac monopole charge and a variable integer number xyz units
310 : ///of electric charge.
311 : ///
312 : ///Codes 411xyz0 are then used when the magnetic and electrical
313 : ///charge sign agree and 412xyz0 when they disagree,
314 : ///with the overall sign of the particle set by the magnetic charge.
315 : ///For now no spin information is provided.
316 : ///
317 3329523 : if( extraBits(pid) > 0 ) { return false; }
318 3329523 : if( digit(n,pid) != 4 ) { return false; }
319 0 : if( digit(nr,pid) != 1 ) { return false; }
320 0 : if( (digit(nl,pid) != 1) && (digit(nl,pid) != 2) ) { return false; }
321 : // All Dyons have at least 1 core digit
322 0 : if( digit(nq3,pid) == 0 ) { return false; }
323 : // Dyons have spin zero for now
324 0 : if( digit(nj,pid) != 0 ) { return false; }
325 : return true;
326 : }
327 :
328 : // Check for QBalls
329 : // Ad-hoc numbering for such particles is 100xxxx0,
330 : // where xxxx is the charge in tenths.
331 23521378 : bool isQBall( const int & pid )
332 : {
333 23521378 : if( extraBits(pid) != 1 ) { return false; }
334 0 : if( digit(n,pid) != 0 ) { return false; }
335 0 : if( digit(nr,pid) != 0 ) { return false; }
336 : // check the core number
337 0 : if( (abspid(pid)/10)%10000 == 0 ) { return false; }
338 : // these particles have spin zero for now
339 0 : if( digit(nj,pid) != 0 ) { return false; }
340 : return true;
341 : }
342 :
343 : // does this particle contain an up quark?
344 0 : bool hasUp( const int & pid)
345 : {
346 0 : if( extraBits(pid) > 0 ) { return false; }
347 0 : if( fundamentalID(pid) > 0 ) { return false; }
348 0 : return findQ(pid,2);
349 : }
350 : // does this particle contain a down quark?
351 0 : bool hasDown( const int & pid)
352 : {
353 0 : if( extraBits(pid) > 0 ) { return false; }
354 0 : if( fundamentalID(pid) > 0 ) { return false; }
355 0 : return findQ(pid,1);
356 : }
357 : // does this particle contain a strange quark?
358 0 : bool hasStrange( const int & pid )
359 : {
360 0 : if( extraBits(pid) > 0 ) { return false; }
361 0 : if( fundamentalID(pid) > 0 ) { return false; }
362 0 : return findQ(pid,3);
363 : }
364 : // does this particle contain a charm quark?
365 0 : bool hasCharm( const int & pid )
366 : {
367 0 : if( extraBits(pid) > 0 ) { return false; }
368 0 : if( fundamentalID(pid) > 0 ) { return false; }
369 0 : return findQ(pid,4);
370 : }
371 : // does this particle contain a bottom quark?
372 0 : bool hasBottom( const int & pid )
373 : {
374 0 : if( extraBits(pid) > 0 ) { return false; }
375 0 : if( fundamentalID(pid) > 0 ) { return false; }
376 0 : return findQ(pid,5);
377 : }
378 : // does this particle contain a top quark?
379 0 : bool hasTop( const int & pid )
380 : {
381 0 : if( extraBits(pid) > 0 ) { return false; }
382 0 : if( fundamentalID(pid) > 0 ) { return false; }
383 0 : return findQ(pid,6);
384 : }
385 :
386 : // --- other information
387 : //
388 : // jSpin returns 2J+1, where J is the total spin
389 0 : int jSpin( const int & pid )
390 : {
391 0 : if( fundamentalID(pid) > 0 ) {
392 : // some of these are known
393 0 : int fund = fundamentalID(pid);
394 0 : if( fund > 0 && fund < 7 ) return 2;
395 0 : if( fund == 9 ) return 3;
396 0 : if( fund > 10 && fund < 17 ) return 2;
397 0 : if( fund > 20 && fund < 25 ) return 3;
398 : return 0;
399 0 : } else if( extraBits(pid) > 0 ) {
400 : return 0;
401 : }
402 0 : return abspid(pid)%10;
403 : }
404 : // sSpin returns 2S+1, where S is the spin
405 0 : int sSpin( const int & pid )
406 : {
407 0 : if( !isMeson(pid) ) { return 0; }
408 0 : int inl = digit(nl,pid);
409 : //int tent = digit(n,pid);
410 0 : int js = digit(nj,pid);
411 0 : if( digit(n,pid) == 9 ) { return 0; } // tentative ID
412 : //if( tent == 9 ) { return 0; } // tentative assignment
413 0 : if( inl == 0 && js >= 3 ) {
414 : return 1;
415 0 : } else if( inl == 0 && js == 1 ) {
416 : return 0;
417 0 : } else if( inl == 1 && js >= 3 ) {
418 : return 0;
419 0 : } else if( inl == 2 && js >= 3 ) {
420 : return 1;
421 0 : } else if( inl == 1 && js == 1 ) {
422 : return 1;
423 0 : } else if( inl == 3 && js >= 3 ) {
424 : return 1;
425 : }
426 : // default to zero
427 : return 0;
428 : }
429 : // lSpin returns 2L+1, where L is the orbital angular momentum
430 0 : int lSpin( const int & pid )
431 : {
432 0 : if( !isMeson(pid) ) { return 0; }
433 0 : int inl = digit(nl,pid);
434 : //int tent = digit(n,pid);
435 0 : int js = digit(nj,pid);
436 0 : if( digit(n,pid) == 9 ) { return 0; } // tentative ID
437 0 : if( inl == 0 && js == 3 ) {
438 : return 0;
439 0 : } else if( inl == 0 && js == 5 ) {
440 : return 1;
441 0 : } else if( inl == 0 && js == 7 ) {
442 : return 2;
443 0 : } else if( inl == 0 && js == 9 ) {
444 : return 3;
445 0 : } else if( inl == 0 && js == 1 ) {
446 : return 0;
447 0 : } else if( inl == 1 && js == 3 ) {
448 : return 1;
449 0 : } else if( inl == 1 && js == 5 ) {
450 : return 2;
451 0 : } else if( inl == 1 && js == 7 ) {
452 : return 3;
453 0 : } else if( inl == 1 && js == 9 ) {
454 : return 4;
455 0 : } else if( inl == 2 && js == 3 ) {
456 : return 1;
457 0 : } else if( inl == 2 && js == 5 ) {
458 : return 2;
459 0 : } else if( inl == 2 && js == 7 ) {
460 : return 3;
461 0 : } else if( inl == 2 && js == 9 ) {
462 : return 4;
463 0 : } else if( inl == 1 && js == 1 ) {
464 : return 1;
465 0 : } else if( inl == 3 && js == 3 ) {
466 : return 2;
467 0 : } else if( inl == 3 && js == 5 ) {
468 : return 3;
469 0 : } else if( inl == 3 && js == 7 ) {
470 : return 4;
471 0 : } else if( inl == 3 && js == 9 ) {
472 : return 5;
473 : }
474 : // default to zero
475 : return 0;
476 : }
477 :
478 : // 3 times the charge
479 20191855 : int threeCharge( const int & pid )
480 : {
481 : int charge=0;
482 : int ida, sid;
483 : unsigned short q1, q2, q3, ql;
484 : static int ch100[] = { -1, 2,-1, 2,-1, 2,-1, 2, 0, 0,
485 : -3, 0,-3, 0,-3, 0,-3, 0, 0, 0,
486 : 0, 0, 0, 3, 0, 0, 0, 0, 0, 0,
487 : 0, 0, 0, 3, 0, 0, 3, 0, 0, 0,
488 : 0, -1, 0, 0, 0, 0, 0, 0, 0, 0,
489 : 0, 6, 3, 6, 0, 0, 0, 0, 0, 0,
490 : 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
491 : 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
492 : 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
493 : 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
494 20191855 : q1 = digit(nq1,pid);
495 20191852 : q2 = digit(nq2,pid);
496 20191852 : q3 = digit(nq3,pid);
497 20191851 : ql = digit(nl,pid);
498 20191851 : ida = abspid(pid);
499 20191853 : sid = fundamentalID(pid);
500 20191854 : if( ida == 0 ) { // illegal
501 : return 0;
502 3329524 : } else if( isQBall(pid) ) { // QBall
503 0 : charge = 3*((abspid(pid)/10)%10000);
504 3329524 : } else if( extraBits(pid) > 0 ) { // ion
505 : return 0;
506 3329523 : } else if( isDyon(pid) ) { // Dyon
507 0 : charge = 3*( (abspid(pid)/10)%1000 );
508 : // this is half right
509 : // the charge sign will be changed below if pid < 0
510 0 : if( ql == 2 ) {
511 0 : charge = -charge;
512 : }
513 3329523 : } else if( sid > 0 && sid <= 100 ) { // use table
514 3329523 : charge = ch100[sid-1];
515 3329523 : if(ida==1000017 || ida==1000018) { charge = 0; }
516 3329523 : if(ida==1000034 || ida==1000052) { charge = 0; }
517 3329523 : if(ida==1000053 || ida==1000054) { charge = 0; }
518 3329523 : if(ida==5100061 || ida==5100062) { charge = 6; }
519 0 : } else if( digit(nj,pid) == 0 ) { // KL, Ks, or undefined
520 : return 0;
521 0 : } else if( isMeson(pid) ) { // mesons
522 0 : if( q2 == 3 || q2 == 5 ) {
523 0 : charge = ch100[q3-1] - ch100[q2-1];
524 : } else {
525 0 : charge = ch100[q2-1] - ch100[q3-1];
526 : }
527 0 : } else if( isRhadron(pid) ) { // Rhadron
528 0 : if (( q1 == 0 ) || ( q1 == 9 )) {
529 0 : if( q2 == 3 || q2 == 5 ) {
530 0 : charge = ch100[q3-1] - ch100[q2-1];
531 : } else {
532 0 : charge = ch100[q2-1] - ch100[q3-1];
533 : }
534 0 : } else if( ql == 0 ) {
535 0 : charge = ch100[q3-1] + ch100[q2-1] + ch100[q1-1];
536 0 : } else if ( digit(nr,pid) == 0 ) {
537 0 : charge = ch100[q3-1] + ch100[q2-1] + ch100[q1-1] + ch100[ql-1];
538 : }
539 0 : } else if( isDiQuark(pid) ) { // diquarks
540 0 : charge = ch100[q2-1] + ch100[q1-1];
541 0 : } else if( isBaryon(pid) ) { // baryons
542 0 : charge = ch100[q3-1] + ch100[q2-1] + ch100[q1-1];
543 : } else { // unknown
544 : return 0;
545 : }
546 3329523 : if( charge == 0 ) {
547 3315181 : return 0;
548 14342 : } else if( pid < 0 ) {
549 6739 : charge = -charge;
550 : }
551 : return charge;
552 : }
553 :
554 : // the actual charge
555 20191857 : double charge( const int & pid )
556 : {
557 20191857 : int tc = threeCharge(pid);
558 20191854 : if( isQBall(pid) ) {
559 0 : return double(tc)/30.;
560 : } else {
561 20191854 : return double(tc)/3.;
562 : }
563 : }
564 :
565 : } // HepPID
|
