Context Navigation

source: trunk/source/processes/hadronic/models/chiral_inv_phase_space/body/src/G4QPDGCode.cc @ 1071

Last change on this file since 1071 was 1055, checked in by garnier, 15 years ago
maj sur la beta de geant 4.9.3
File size: 88.1 KB

Line
1	//
2	// ********************************************************************
3	// * License and Disclaimer *
4	// * *
5	// * The Geant4 software is copyright of the Copyright Holders of *
6	// * the Geant4 Collaboration. It is provided under the terms and *
7	// * conditions of the Geant4 Software License, included in the file *
8	// * LICENSE and available at http://cern.ch/geant4/license . These *
9	// * include a list of copyright holders. *
10	// * *
11	// * Neither the authors of this software system, nor their employing *
12	// * institutes,nor the agencies providing financial support for this *
13	// * work make any representation or warranty, express or implied, *
14	// * regarding this software system or assume any liability for its *
15	// * use. Please see the license in the file LICENSE and URL above *
16	// * for the full disclaimer and the limitation of liability. *
17	// * *
18	// * This code implementation is the result of the scientific and *
19	// * technical work of the GEANT4 collaboration. *
20	// * By using, copying, modifying or distributing the software (or *
21	// * any work based on the software) you agree to acknowledge its *
22	// * use in resulting scientific publications, and indicate your *
23	// * acceptance of all terms of the Geant4 Software license. *
24	// ********************************************************************
25	//
26	//
27	// $Id: G4QPDGCode.cc,v 1.57 2009/02/23 09:49:24 mkossov Exp $
28	// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
29	//
30	// ---------------- G4QPDGCode ----------------
31	// by Mikhail Kossov, Sept 1999.
32	// class for Hadron definitions in CHIPS Model
33	// -------------------------------------------------------------------
34	// Short description: The PDG Code is made on the basis of the Quark
35	// Content (G4QuarkContent) of the hadronic state (including nuclear
36	// fragments). The PDG code of the ground state (e.g. pi, N, etc.) is
37	// calculated. It includes a complicated algortithm of the G.S. mass
38	// calculation for nuclear fragments (now it is synchronised with the
39	// G4 nuclear massess).
40	// -------------------------------------------------------------------
41
42	//#define debug
43	//#define pdebug
44	//#define qdebug
45	//#define idebug
46	//#define sdebug
47
48	#include "G4QPDGCodeVector.hh"
49	#include <cmath>
50	#include <cstdlib>
51	using namespace std;
52
53	G4QPDGCode::G4QPDGCode(G4int PDGCode): thePDGCode(PDGCode)
54	{
55	#ifdef sdebug
56	G4cout<<"G4QPDGCode:Constructer is called with PDGCode="<<PDGCode<<G4endl;
57	#endif
58	if(PDGCode) theQCode=MakeQCode(PDGCode);
59	else
60	{
61	#ifdef sdebug
62	G4cout<<"***G4QPDGCode: Constructed with PDGCode=0, QCode=-2"<<G4endl;
63	#endif
64	theQCode=-2;
65	}
66	#ifdef debug
67	G4cout<<"G4QPDGCode:Constructer(PDG) PDG="<<PDGCode<<", QCode="<<theQCode<<G4endl;
68	#endif
69	}
70
71	G4QPDGCode::G4QPDGCode(G4bool f, G4int QCode): theQCode(QCode)
72	{
73	if(f&&QCode<0)G4cerr<<"***G4QPDGCode::Constr. QCode="<<QCode<<G4endl;
74	thePDGCode = MakePDGCode(QCode);
75	#ifdef debug
76	G4cout<<"G4QPDGCode::Constr: PDGCode="<<thePDGCode<<G4endl;
77	#endif
78	}
79
80	G4QPDGCode::G4QPDGCode(G4QContent QCont)
81	{
82	InitByQCont(QCont);
83	}
84
85	G4QPDGCode::G4QPDGCode(const G4QPDGCode& rhs)
86	{
87	thePDGCode =rhs.thePDGCode;
88	theQCode =rhs.theQCode;
89	}
90
91	G4QPDGCode::G4QPDGCode(G4QPDGCode* rhs)
92	{
93	thePDGCode =rhs->thePDGCode;
94	theQCode =rhs->theQCode;
95	}
96
97	const G4QPDGCode& G4QPDGCode::operator=(const G4QPDGCode& rhs)
98	{
99	if(this != &rhs) // Beware of self assignment
100	{
101	thePDGCode =rhs.thePDGCode;
102	theQCode =rhs.theQCode;
103	}
104	return *this;
105	}
106
107	G4QPDGCode::~G4QPDGCode() {}
108
109	//G4int G4QPDGCode::nQHM=90;
110
111	// Standard output for QPDGCode
112	ostream& operator<<(ostream& lhs, G4QPDGCode& rhs)
113	// =========================================
114	{
115	lhs << "[ PDG=" << rhs.GetPDGCode() << ", Q=" << rhs.GetQCode() << "]";
116	return lhs;
117	}
118
119	// Standard output for const QPDGCode
120	ostream& operator<<(ostream& lhs, const G4QPDGCode& rhs)
121	// ===============================================
122	{
123	lhs << "[ PDG=" << rhs.GetPDGCode() << ", Q=" << rhs.GetQCode() << "]";
124	return lhs;
125	}
126
127	// Overloading of QPDGCode addition
128	G4int operator+(const G4QPDGCode& lhs, const G4QPDGCode& rhs)
129	// =======================================================
130	{
131	G4int s = lhs.GetPDGCode();
132	return s += rhs.GetPDGCode();
133	}
134	G4int operator+(const G4QPDGCode& lhs, const G4int& rhs)
135	// =======================================================
136	{
137	G4int s = lhs.GetPDGCode();
138	return s += rhs;
139	}
140	G4int operator+(const G4int& lhs, const G4QPDGCode& rhs)
141	// =======================================================
142	{
143	G4int s = lhs;
144	return s += rhs.GetPDGCode();
145	}
146
147	// Overloading of QPDGCode subtraction
148	G4int operator-(const G4QPDGCode& lhs, const G4QPDGCode& rhs)
149	// =======================================================
150	{
151	G4int s = lhs.GetPDGCode();
152	return s -= rhs.GetPDGCode();
153	}
154	G4int operator-(const G4QPDGCode& lhs, const G4int& rhs)
155	// =======================================================
156	{
157	G4int s = lhs.GetPDGCode();
158	return s -= rhs;
159	}
160	G4int operator-(const G4int& lhs, const G4QPDGCode& rhs)
161	// =======================================================
162	{
163	G4int s = lhs;
164	return s -= rhs.GetPDGCode();
165	}
166
167	// Overloading of QPDGCode multiplication
168	G4int operator*(const G4QPDGCode& lhs, const G4QPDGCode& rhs)
169	// =======================================================
170	{
171	G4int s = lhs.GetPDGCode();
172	return s *= rhs.GetPDGCode();
173	}
174	G4int operator*(const G4QPDGCode& lhs, const G4int& rhs)
175	// =======================================================
176	{
177	G4int s = lhs.GetPDGCode();
178	return s *= rhs;
179	}
180	G4int operator*(const G4int& lhs, const G4QPDGCode& rhs)
181	// =======================================================
182	{
183	G4int s = lhs;
184	return s *= rhs.GetPDGCode();
185	}
186
187	// Overloading of QPDGCode division
188	G4int operator/(const G4QPDGCode& lhs, const G4QPDGCode& rhs)
189	{// =======================================================
190	G4int s = lhs.GetPDGCode();
191	return s /= rhs.GetPDGCode();
192	}
193	G4int operator/(const G4QPDGCode& lhs, const G4int& rhs)
194	{// ==================================================
195	G4int s = lhs.GetPDGCode();
196	return s /= rhs;
197	}
198	G4int operator/(const G4int& lhs, const G4QPDGCode& rhs)
199	{// ==================================================
200	G4int s = lhs;
201	return s /= rhs.GetPDGCode();
202	}
203
204	// Overloading of QPDGCode residual
205	G4int operator%(const G4QPDGCode& lhs, const G4int& rhs)
206	{// ==================================================
207	G4int s = lhs.GetPDGCode();
208	return s %= rhs;
209	}
210
211	// TRUE if it is not RealNeutral (111,221,331 etc), FALSE if it is.
212	G4bool G4QPDGCode::TestRealNeutral(const G4int& PDGCode)
213	{// =================================================
214	if(PDGCode>0 && PDGCode<999) // RealNeutral are always positive && mesons
215	{
216	if(PDGCode==22) return false; // Photon
217	G4int p=PDGCode/10;
218	if(p/10==p%10) return false; // This is a RealNeutral
219	}
220	return true;
221	}
222
223	// Make a Q Code out of the PDG Code
224	G4int G4QPDGCode::MakePDGCode(const G4int& QCode)
225	{// ===========================================
226	//static const G4int modi = 81; // Q Codes for more than di-baryon nuclei
227	//static const G4int modi = 89; // Q Codes for more than di-baryon nuclei "IsoNuclei"
228	static const G4int modi = 122; // Q Codes for more than quarta-baryon nuclei "Lept/Hyper"
229	static G4int qC[modi] = { 11, 12, 13, 14, 15, 16, 22, 23, 24, 25, // 10
230	37, 110, 220, 330, 111, 211, 221, 311, 321, 331, // 20
231	2112, 2212, 3122, 3112, 3212, 3222, 3312, 3322, 113, 213, // 30
232	223, 313, 323, 333, 1114, 2114, 2214, 2224, 3124, 3114, // 40
233	3214, 3224, 3314, 3324, 3334, 115, 215, 225, 315, 325, // 50
234	335, 2116, 2216, 3126, 3116, 3216, 3226, 3316, 3326, 117, // 60
235	217, 227, 317, 327, 337, 1118, 2118, 2218, 2228, 3128, // 70
236	3118, 3218, 3228, 3318, 3328, 3338, 119, 219, 229, 319, // 80
237	329, 339, 90002999 , 89999003 , 90003998 ,
238	89998004 , 90003999 , 89999004 , 90004998 , 89998005 , // 90
239	90000001 , 90001000 , 91000000 , 90999001 , 91000999 ,
240	91999000 , 91999999 , 92999000 , 90000002 , 90001001 , //100
241	90002000 , 91000001 , 91001000 , 92000000 , 90999002 ,
242	91001999 , 90001002 , 90002001 , 91000002 , 91001001 , //110
243	91002000 , 92000001 , 92001000 , 90999003 , 90001003 ,
244	90002002 , 90003001 , 91001002 , 91002001 , 92000002 , //120
245	92001001 , 92002000};
246	static G4int aC[15] = {1,1000,999001,1000000,1000999,1999000,1999999, // sum 1
247	2,1001,2000,1000001,1001000,1999001,2000000,2000999}; // sum 2
248	if (QCode<0)
249	{
250	G4cerr<<"***G4QPDGCode::MakePDGCode: negative Q Code ="<<QCode<<G4endl;
251	return 0;
252	}
253	else if (QCode>=modi)
254	{
255	//G4int q=QCode-modi; // Starting BarNum=3
256	//G4int a=q/15+1; // BarNum/2
257	//G4int b=q%15;
258	G4int q=QCode-modi; // Starting BarNum=5
259	G4int a=q/15+2; // BarNum/2
260	G4int b=q%15;
261	return 90000000+a*1001+aC[b];
262	}
263	return qC[theQCode];
264	}
265
266	// Hadronic masses synhronized with the Geant4 hadronic masses
267	G4double G4QPDGCode:: QHaM(G4int nQ)
268	{// ===========================
269	static G4bool iniFlag=true;
270	static G4double m[nQHM]={.511, 0., 105.65837, 0., 1777., 0., 0., 91.188, 80.425, 140.00
271	,120.000, 800., 980., 1370., 134.98, 139.57, 547.75, 497.65, 493.68, 957.78
272	,939.5654,938.272, 1115.683, 1197.45, 1192.64, 1189.37, 1321.3, 1314.8, 775.8, 775.8
273	, 782.6, 896.1, 891.66, 1019.456, 1232., 1232., 1232., 1232., 1519.5, 1387.2
274	, 1383.7, 1382.8, 1535., 1531.8, 1672.45, 1318.3, 1318.3, 1275.4, 1432.4, 1425.6
275	, 1525., 1680., 1680., 1820., 1915., 1915., 1915., 2025., 2025., 1691.
276	, 1691., 1667., 1776., 1776., 1854., 1950., 1950., 1950., 1950., 2100.
277	, 2030., 2030., 2030., 2127., 2127., 2252., 2020., 2020., 2044., 2045.
278	, 2045., 2297., 2170.272, 2171.565, 2464., 2464., 3108.544, 3111.13,3402.272,3403.565};
279	if(iniFlag) // Initialization of the Geant4 hadronic masses
280	{
281	m[ 0]= G4Electron::Electron()->GetPDGMass();
282	m[ 1]= G4NeutrinoE::NeutrinoE()->GetPDGMass();
283	m[ 2]= G4MuonMinus::MuonMinus()->GetPDGMass();
284	m[ 3]= G4NeutrinoMu::NeutrinoMu()->GetPDGMass();
285	m[ 4]= G4TauMinus::TauMinus()->GetPDGMass();
286	m[ 5]=G4NeutrinoTau::NeutrinoTau()->GetPDGMass();
287	m[14]= G4PionZero::PionZero()->GetPDGMass();
288	m[15]= G4PionMinus::PionMinus()->GetPDGMass();
289	m[16]= G4Eta::Eta()->GetPDGMass();
290	m[17]= G4KaonZero::KaonZero()->GetPDGMass();
291	m[18]= G4KaonMinus::KaonMinus()->GetPDGMass();
292	m[19]= G4EtaPrime::EtaPrime()->GetPDGMass();
293	m[20]= G4Neutron::Neutron()->GetPDGMass();
294	m[21]= G4Proton::Proton()->GetPDGMass();
295	m[22]= G4Lambda::Lambda()->GetPDGMass();
296	m[23]= G4SigmaMinus::SigmaMinus()->GetPDGMass();
297	m[24]= G4SigmaZero::SigmaZero()->GetPDGMass();
298	m[25]= G4SigmaPlus::SigmaPlus()->GetPDGMass();
299	m[26]= G4XiMinus::XiMinus()->GetPDGMass();
300	m[27]= G4XiZero::XiZero()->GetPDGMass();
301	m[44]= G4OmegaMinus::OmegaMinus()->GetPDGMass();
302	iniFlag=false;
303	}
304	if(nQ<0 \|\| nQ>=nQHM)
305	{
306	G4cout<<"***G4QPDGCode::QHaM: negative Q-code or Q="<<nQ<<" >= nQmax = "<<nQHM<<G4endl;
307	return 0.;
308	}
309	return m[nQ];
310	}
311
312	// Make a Q Code out of the PDG Code
313	G4int G4QPDGCode::MakeQCode(const G4int& PDGCode)
314	{// ===========================================
315	static const G4int qr[10]={0,13,19,27,33,44,50,58,64,75};
316	G4int PDGC=abs(PDGCode); // Qcode is always not negative
317	G4int s=0;
318	G4int z=0;
319	G4int n=0;
320	if (PDGC>100000000) // Not supported
321	{
322	#ifdef debug
323	G4cout<<"***G4QPDGCode::MakeQCode: Unknown in Q-System code: "<<PDGCode<<G4endl;
324	#endif
325	return -2;
326	}
327	else if (PDGC>80000000&&PDGC<100000000) // Try to convert the NUCCoding to PDGCoding
328	{
329	if(PDGC==90000000) return 6;
330	ConvertPDGToZNS(PDGC, z, n, s);
331	G4int b=n+z+s; // Baryon number
332	#ifdef debug
333	G4cout<<"***G4QPDGCode::Z="<<z<<",N="<<n<<",S="<<s<<G4endl;
334	#endif
335	if(b<0) // ---> Baryons & Fragments
336	{
337	b=-b;
338	n=-n;
339	z=-z;
340	s=-s;
341	PDGC=90000000+s1000000+z1000+n; // New PDGC for anti-baryons
342	}
343	else if(!b) // --> Mesons
344	{
345	//G4bool anti=false; // For the PDG conversion
346	if(z<0) // --> Mesons conversion
347	{
348	n=-n;
349	z=-z;
350	s=-s;
351	//anti=true; // For the PDG conversion
352	}
353	if(!z)
354	{
355	if(s>0)
356	{
357	n=-n;
358	s=-s;
359	//anti=true; // For the PDG conversion
360	}
361	if (s==-1) return 17; // K0
362	else if(s==-2) return -1; // K0+K0 chipolino
363	else return -2; // Not supported by Q Code
364	}
365	else // --> z>0
366	{
367	if(z==1)
368	{
369	if (s==-1) return 18; // K+
370	else return 15; // pi+
371	}
372	else if(z==2) return -1; // Chipolino
373	else return -2; // Not supported by Q Code
374	}
375	} // End of meson case
376	if(b>0) // --> Baryon case
377	{
378	if(b==1)
379	{
380	if(!s) // --> Baryons
381	{
382	if(z==-1) return 34; // Delta-
383	else if(!z) return 91; // neutron
384	else if(z==1)return 91; // proton
385	else if(z==2)return 37; // Delta++
386	else if(z==3\|\|z==-2)return -1; // Delta+pi Chipolino
387	else return -2; // Not supported by Q Code
388	}
389	else if(s==1) // --> Hyperons
390	{
391	if(z==-1) return 93; // Sigma-
392	else if(!z) return 92; // Lambda (@@ 24->Sigma0)
393	else if(z==1)return 94; // Sigma+
394	else if(z==2\|\|z==-2) return -1; // Sigma+pi Chipolino
395	else return -2; // Not supported by Q Code
396	}
397	else if(s==2) // --> Xi Hyperons
398	{
399	if(z==-1) return 95; // Xi-
400	else if(!z) return 96; // Xi0
401	else if(z==1\|\|z==-2)return -1; // Xi+pi Chipolino
402	else return -2; // Not supported by Q Code
403	}
404	else if(s==3) // --> Xi Hyperons
405	{
406	if(z==-1) return 97; // Omega-
407	else if(!z\|\|z==-2) return -1; // Omega+pi Chipolino
408	else return -2; // Not supported by Q Code
409	}
410	}
411	else
412	{
413	if(b==2)
414	{
415	if (PDGC==90002999) return 82; // p DEL++
416	else if(PDGC==89999003) return 83; // n DEL-
417	else if(PDGC==90003998) return 84; // DEL++ DEL++
418	else if(PDGC==89998004) return 85; // DEL- DEL-
419	else if(PDGC==90999002) return 104; // n Sigma-
420	else if(PDGC==91001999) return 105; // p Sigma+
421	}
422	if(b==3)
423	{
424	if (PDGC==90003999) return 86; // p p DEL++
425	else if(PDGC==89999004) return 87; // n n DEL-
426	else if(PDGC==90004998) return 88; // p DEL++ DEL++
427	else if(PDGC==89998005) return 89; // n DEL- DEL-
428	else if(PDGC==90999003) return 113; // n n Sigma-
429	}
430	}
431	}
432	}
433	if (PDGC<80000000) // ----> Direct Baryons & Mesons
434	{
435	if (PDGC<100) // => Leptons and field bosons
436	{
437	if (PDGC==10) return -1; // Chipolino
438	else if(PDGC==11) return 0; // e-
439	else if(PDGC==12) return 1; // nu_e
440	else if(PDGC==13) return 2; // mu-
441	else if(PDGC==14) return 3; // nu_mu
442	else if(PDGC==15) return 4; // tau-
443	else if(PDGC==16) return 5; // nu_tau
444	else if(PDGC==22) return 6; // Photon
445	else if(PDGC==23) return 7; // Z0 boson
446	else if(PDGC==24) return 8; // W- boson
447	else if(PDGC==25) return 9; // H0 (neutral Higs boson)
448	else if(PDGC==37) return 10; // H- (charged Higs boson)
449	}
450	G4int r=PDGC%10; // 2s+1
451	G4int Q= 0;
452	if (!r)
453	{
454	// Internal CHIPS codes for the wide f_0 states must be 9000221, 9010221, 10221
455	if (PDGC==110) return 11; // Low R-P: Sigma (pi,pi S-wave)
456	else if(PDGC==220) return 12; // Midle Regeon-Pomeron
457	else if(PDGC==330) return 13; // High Regeon-Pomeron
458	#ifdef pdebug
459	G4cout<<"***G4QPDGCode::MakeQCode: (0) Unknown in Q-System code: "<<PDGCode<<G4endl;
460	#endif
461	return -2;
462	}
463	else Q=qr[r];
464	G4int p=PDGC/10; // Quark Content
465	if(r%2) // (2s+1 is odd) Mesons are all the same
466	{
467	if (p==11) return Q+=1;
468	else if(p==21) return Q+=2;
469	else if(p==22) return Q+=3;
470	else if(p==31) return Q+=4;
471	else if(p==32) return Q+=5;
472	else if(p==33) return Q+=6;
473	else
474	{
475	#ifdef pdebug
476	G4cout<<"***G4QPDGCode::MakeQCode:(1) Unknown in Q-System code: "<<PDGCode<<G4endl;
477	#endif
478	return -2;
479	}
480	}
481	else // (2s+1 is even) Baryons
482	{
483	G4int s=r/2;
484	if(s%2) // ((2s+1)/2 is odd) N Family
485	{
486	if (p==211) return Q+=1;
487	else if(p==221) return Q+=2;
488	else if(p==312) return Q+=3;
489	else if(p==311) return Q+=4;
490	else if(p==321) return Q+=5;
491	else if(p==322) return Q+=6;
492	else if(p==331) return Q+=7;
493	else if(p==332) return Q+=8;
494	else
495	{
496	#ifdef pdebug
497	G4cout<<"**G4QPDGCode::MakeQCode:(2) Unknown in Q-System code:"<<PDGCode<<G4endl;
498	#endif
499	return -2;
500	}
501	}
502	else // ((2s+1)/2 is odd) Delta Family
503	{
504	if (p==111) return Q+= 1;
505	else if(p==211) return Q+= 2;
506	else if(p==221) return Q+= 3;
507	else if(p==222) return Q+= 4;
508	else if(p==312) return Q+= 5;
509	else if(p==311) return Q+= 6;
510	else if(p==321) return Q+= 7;
511	else if(p==322) return Q+= 8;
512	else if(p==331) return Q+= 9;
513	else if(p==332) return Q+=10;
514	else if(p==333) return Q+=11;
515	else
516	{
517	#ifdef pdebug
518	G4cout<<"**G4QPDGCode::MakeQCode:(3) Unknown in Q-System code:"<<PDGCode<<G4endl;
519	#endif
520	return -2;
521	}
522	}
523	}
524	}
525	else // Nuclear Fragments
526	{
527	G4int d=n+n+z+s; // a#of d quarks
528	G4int u=n+z+z+s; // a#of u quarks
529	G4int t=d+u+s; // tot#of quarks
530	if(t%3 \|\| abs(t)<3) // b=0 are in mesons
531	{
532	#ifdef pdebug
533	G4cout<<"***G4QPDGCode::MakeQCode: Unknown PDGCode="<<PDGCode<<", t="<<t<<G4endl;
534	#endif
535	return -2;
536	}
537	else
538	{
539	G4int b=t/3; // baryon number
540	if(b==1) // baryons
541	{
542	if (s==0&&u==1&&d==2) return 90;
543	else if(s==0&&u==2&&d==1) return 91;
544	else if(s==1&&u==1&&d==1) return 92;
545	else
546	{
547	#ifdef pdebug
548	G4cout<<"**G4QPDGCode::MakeQCode:(5) Unknown in Q-System code:"<<PDGCode<<G4endl;
549	#endif
550	return -2;
551	}
552	}
553	else if(b==2) // di-baryons
554	{
555	if (s==0&&u==2&&d==4) return 98;
556	else if(s==0&&u==3&&d==3) return 99;
557	else if(s==0&&u==4&&d==2) return 100;
558	else if(s==1&&u==2&&d==3) return 101;
559	else if(s==1&&u==3&&d==2) return 102;
560	else if(s==2&&u==2&&d==2) return 103;
561	else
562	{
563	#ifdef pdebug
564	G4cout<<"**G4QPDGCode::MakeQCode:(6) Unknown in Q-System code:"<<PDGCode<<G4endl;
565	#endif
566	return -2;
567	}
568	}
569	else if(b==3) // tri-baryons
570	{
571	if (s==0&&u==4&&d==5) return 106; // pnn
572	else if(s==0&&u==5&&d==4) return 107; // npp
573	else if(s==1&&u==3&&d==5) return 108; // Lnn
574	else if(s==1&&u==4&&d==4) return 109; // Lnp
575	else if(s==1&&u==5&&d==3) return 110; // Lpp
576	else if(s==2&&u==3&&d==4) return 111; // LLn
577	else if(s==2&&u==4&&d==3) return 112; // LLp
578	else if(s==1&&u==2&&d==6) return 113; // SIG-nn
579	else
580	{
581	#ifdef pdebug
582	G4cout<<"**G4QPDGCode::MakeQCode:(7) Unknown in Q-System code:"<<PDGCode<<G4endl;
583	#endif
584	return -2;
585	}
586	}
587	G4int c=b/2;
588	G4int g=b%2;
589	G4int h=c*3;
590	//G4int Q=57+c*15;
591	//G4int Q=65+c*15; // "IsoNuclei"
592	G4int Q=83+c*15; // "Leptons/Hyperons"
593	u-=h;
594	d-=h;
595	if(g)
596	{
597	if (s==0&&u==1&&d==2) return Q+= 9;
598	else if(s==0&&u==2&&d==1) return Q+=10;
599	else if(s==1&&u==0&&d==2) return Q+=11;
600	else if(s==1&&u==1&&d==1) return Q+=12;
601	else if(s==1&&u==2&&d==0) return Q+=13;
602	else if(s==2&&u==0&&d==1) return Q+=14;
603	else if(s==2&&u==1&&d==0) return Q+=15;
604	else
605	{
606	#ifdef debug
607	G4cout<<"**G4QPDGCode::MakeQCode:(8) Unknown in Q-System code:"<<PDGCode<<G4endl;
608	#endif
609	return -2;
610	}
611	}
612	else
613	{
614	if (s==0&&u==-1&&d== 1) return Q+=1;
615	else if(s==0&&u== 0&&d== 0) return Q+=2;
616	else if(s==0&&u== 1&&d==-1) return Q+=3;
617	else if(s==1&&u==-1&&d== 0) return Q+=4;
618	else if(s==1&&u== 0&&d==-1) return Q+=5;
619	else if(s==2&&u==-2&&d== 0) return Q+=6;
620	else if(s==2&&u==-1&&d==-1) return Q+=7;
621	else if(s==2&&u== 0&&d==-2) return Q+=8;
622	else
623	{
624	#ifdef debug
625	G4cout<<"**G4QPDGCode::MakeQCode:(9) Unknown in Q-System code:"<<PDGCode<<G4endl;
626	#endif
627	return -2;
628	}
629	}
630	}
631	}
632	#ifdef pdebug
633	G4cout<<"***G4QPDGCode::MakeQCode: () Unknown in Q-System code: "<<PDGCode<<G4endl;
634	#endif
635	// return -2; not reachable statement
636	}
637
638	// Get the mean mass value for the PDG
639	G4double G4QPDGCode::GetMass()
640	{// =====================
641	G4int ab=theQCode;
642	#ifdef debug
643	G4cout<<"G4QPDGCode::GetMass: Mass for Q="<<ab<<",PDG="<<thePDGCode<<",N="<<nQHM<<G4endl;
644	#endif
645	if ( (ab < 0 && thePDGCode < 80000000) \|\| !thePDGCode) {
646	#ifdef debug
647	if(thePDGCode!=10)
648	G4cout<<"**G4QPDGCode::GetMass:m=100000.,QC="<<theQCode<<",PDG="<<thePDGCode<<G4endl;
649	#endif
650	return 100000.;
651	}
652	else if(ab>-1 && ab<nQHM)
653	{
654	#ifdef debug
655	G4cout<<"G4QPDGCode::GetMa:m="<<QHaM(ab)<<",Q="<<theQCode<<",PDG="<<thePDGCode<<G4endl;
656	#endif
657	return QHaM(ab); // Get mass from the table
658	}
659	//if(szn==0) return m[15]; // @@ mPi0 @@ MK ?
660	if(thePDGCode==90000000)
661	{
662	#ifdef debug
663	G4cout<<"G4QPDGCode::GetMass:***m=0, QC="<<theQCode<<",PDG="<<thePDGCode<<G4endl;
664	#endif
665	return 0.;
666	}
667	G4int s=0;
668	G4int z=0;
669	G4int n=0;
670	ConvertPDGToZNS(thePDGCode, z, n, s);
671	G4double m=GetNuclMass(z,n,s);
672	#ifdef debug
673	G4cout<<"G4QPDG::GetM:PDG="<<thePDGCode<<"=>Z="<<z<<",N="<<n<<",S="<<s<<",M="<<m<<G4endl;
674	#endif
675	return m;
676	}
677
678	// Get the width value for the PDG
679	G4double G4QPDGCode::GetWidth()
680	// =====================
681	{
682	//static const int nW = 72;
683	//static const int nW = 80; // "Isobars"
684	static const int nW = 90; // "Leptons/Hypernuclei"
685	static G4double width[nW] = {0.,0.,0.,0.,0.,0.,0.,2.495,2.118,10.
686	, 10., 800., 75., 350., 0., 0., .00118, 0., 0., .203
687	, 0., 0., 0., 0., 0., 0., 0., 0., 160., 160.
688	, 8.41, 50.5, 50.8, 4.43, 120., 120., 120., 120., 15.6, 39.
689	, 36., 35.8, 10., 9., 0., 107., 107., 185.5, 109., 98.5
690	, 76., 130., 130., 80., 120., 120., 120., 20., 20., 160.
691	, 160., 168., 159., 159., 87., 300., 300., 300., 300., 200.
692	, 180., 180., 180., 99., 99., 55., 387., 387., 208., 198.
693	, 198., 149., 120., 120., 170., 170., 120., 120., 170., 170.};
694	G4int ab=abs(theQCode);
695	if(ab<nW) return width[ab];
696	return 0.; // @@ May be real width should be implemented for nuclei e.g. pp
697	}
698
699	// Get a nuclear mass for Z (a#of protons), N (a#of neutrons), & S (a#of lambdas)
700	G4double G4QPDGCode::GetNuclMass(G4int z, G4int n, G4int s)
701	// ====================================================
702	{
703	static const G4double anb = .01; // Antibinding for Ksi-n,Sig-n,Sig+p,Sig-nn,
704	static const G4double mNeut= QHaM(20);
705	static const G4double mProt= QHaM(21);
706	static const G4double mLamb= QHaM(22);
707	static const G4double mPiC = QHaM(15);
708	static const G4double mKZ = QHaM(17);
709	static const G4double mKM = QHaM(18);
710	static const G4double mSiM = QHaM(23);
711	static const G4double mSiP = QHaM(25);
712	static const G4double mKsZ = QHaM(27);
713	static const G4double mKsM = QHaM(26);
714	static const G4double mOmM = QHaM(44);
715	static const G4double mKZa = mKZ +anb;
716	static const G4double mKMa = mKM +anb;
717	static const G4double mSigM= mSiM+anb;
718	static const G4double mSigP= mSiP+anb;
719	static const G4double mKsiZ= mKsZ+anb;
720	static const G4double mKsiM= mKsM+anb;
721	static const G4double mOmeg= mOmM+anb;
722	static const G4double mDiPi= mPiC+mPiC+anb;
723	static const G4double mDiKZ= mKZa+mKZ;
724	static const G4double mDiKM= mKMa+mKM;
725	static const G4double mDiPr= mProt+mProt;
726	static const G4double mDiNt= mNeut+mNeut;
727	static const G4double mSmPi= mSiM+mDiPi;
728	static const G4double mSpPi= mSiP+mDiPi;
729	static const G4double mOmN = mOmeg+mNeut;
730	static const G4double mSpP = mSigP+mProt;
731	static const G4double mSpPP= mSpP +mProt;
732	static const G4double mSmN = mSigM+mNeut;
733	static const G4double mSmNN= mSmN +mNeut;
734	// -------------- DAM Arrays ----------------------
735	static const G4int iNR=76; // Neutron maximum range for each Z
736	static const G4int nEl = 105; // Maximum Z of the associative memory is "nEl-1=104"
737	static const G4int iNF[nEl]={0,0,0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, // 14
738	1 , 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, // 29
739	16 , 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, // 44
740	31 , 32, 34, 36, 38, 40, 42, 44, 46, 48, 50, 52, 53, 54, 55, // 59
741	56 , 56, 57, 57, 58, 60, 61, 63, 64, 65, 66, 67, 68, 69, 70, // 74
742	71 , 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, // 89
743	86 , 87, 88, 89, 91, 94, 98,103,109,115,122,128,134,140,146};//104
744	#ifdef qdebug
745	static G4int iNmin[nEl]={0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, // 14
746	1 , 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, // 29
747	16 , 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, // 44
748	31 , 32, 34, 36, 38, 40, 42, 44, 46, 48, 50, 52, 53, 54, 55, // 59
749	56 , 56, 57, 57, 58, 60, 61, 63, 64, 65, 66, 67, 68, 69, 70, // 74
750	71 , 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, // 89
751	86 , 87, 88, 89, 91, 94, 98,103,109,115,122,128,134,140,146};//104
752	static G4int iNmax=iNR;
753	static G4int iNran[nEl]={19,20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, // 14
754	34 , 35, 36, 37, 38, 39, 40, 48, 48, 48, 48, 50, 50, 50, 52, // 29
755	53 , 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, // 44
756	68 , 69, 70, 70, 70, 71, 71, 71, 71, 71, 72, 72, 72, 72, 72, // 59
757	73 , 73, 73, 73, 74, 74, 74, 74, 74, 74, 74, 74, 74, 75, 76, // 74
758	76 , 76, 76, 76, 76, 75, 74, 73, 72, 71, 70, 70, 69, 69, 69, // 89
759	68 , 68, 68, 67, 63, 59, 55, 51, 47, 43, 39, 35, 31, 27, 23};//104
760	#endif
761	static const G4int iNL[nEl]={19,20,21,22,23,24, 25, 26, 27, 28, 29, 30, 31, 32, 33, // 14
762	34 , 35, 36, 37, 38, 39, 40, 48, 48, 48, 48, 50, 50, 50, 52, // 29
763	53 , 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, // 44
764	68 , 69, 70, 70, 70, 71, 71, 71, 71, 71, 72, 72, 72, 72, 72, // 59
765	73 , 73, 73, 73, 74, 74, 74, 74, 74, 74, 74, 74, 74, 75, 76, // 74
766	76 , 76, 76, 76, 76, 75, 74, 73, 72, 71, 70, 70, 69, 69, 69, // 89
767	68 , 68, 68, 67, 63, 59, 55, 51, 47, 43, 39, 35, 31, 27, 23};//104
768	// ******* S=-4 vectors ***********
769	static G4bool iNin6[nEl]={false,false,false,false,false,false,false,
770	false,false,false,false,false,false,false,false,false,false,false,false,false,false,
771	false,false,false,false,false,false,false,false,false,false,false,false,false,false,
772	false,false,false,false,false,false,false,false,false,false,false,false,false,false,
773	false,false,false,false,false,false,false,false,false,false,false,false,false,false,
774	false,false,false,false,false,false,false,false,false,false,false,false,false,false,
775	false,false,false,false,false,false,false,false,false,false,false,false,false,false,
776	false,false,false,false,false,false,false,false,false,false,false,false,false,false};
777	static G4double VZ6[nEl][iNR];
778	//******* S=-3 vectors ***********
779	static G4bool iNin7[nEl]={false,false,false,false,false,false,false,
780	false,false,false,false,false,false,false,false,false,false,false,false,false,false,
781	false,false,false,false,false,false,false,false,false,false,false,false,false,false,
782	false,false,false,false,false,false,false,false,false,false,false,false,false,false,
783	false,false,false,false,false,false,false,false,false,false,false,false,false,false,
784	false,false,false,false,false,false,false,false,false,false,false,false,false,false,
785	false,false,false,false,false,false,false,false,false,false,false,false,false,false,
786	false,false,false,false,false,false,false,false,false,false,false,false,false,false};
787	static G4double VZ7[nEl][iNR];
788	// ******* S=-2 vectors ***********
789	static G4bool iNin8[nEl]={false,false,false,false,false,false,false,
790	false,false,false,false,false,false,false,false,false,false,false,false,false,false,
791	false,false,false,false,false,false,false,false,false,false,false,false,false,false,
792	false,false,false,false,false,false,false,false,false,false,false,false,false,false,
793	false,false,false,false,false,false,false,false,false,false,false,false,false,false,
794	false,false,false,false,false,false,false,false,false,false,false,false,false,false,
795	false,false,false,false,false,false,false,false,false,false,false,false,false,false,
796	false,false,false,false,false,false,false,false,false,false,false,false,false,false};
797	static G4double VZ8[nEl][iNR];
798	// ******* S=-1 vectors ***********
799	static G4bool iNin9[nEl]={false,false,false,false,false,false,false,
800	false,false,false,false,false,false,false,false,false,false,false,false,false,false,
801	false,false,false,false,false,false,false,false,false,false,false,false,false,false,
802	false,false,false,false,false,false,false,false,false,false,false,false,false,false,
803	false,false,false,false,false,false,false,false,false,false,false,false,false,false,
804	false,false,false,false,false,false,false,false,false,false,false,false,false,false,
805	false,false,false,false,false,false,false,false,false,false,false,false,false,false,
806	false,false,false,false,false,false,false,false,false,false,false,false,false,false};
807	static G4double VZ9[nEl][iNR];
808	// ******* S=0 vectors ***********
809	static G4bool iNin0[nEl]={false,false,false,false,false,false,false,
810	false,false,false,false,false,false,false,false,false,false,false,false,false,false,
811	false,false,false,false,false,false,false,false,false,false,false,false,false,false,
812	false,false,false,false,false,false,false,false,false,false,false,false,false,false,
813	false,false,false,false,false,false,false,false,false,false,false,false,false,false,
814	false,false,false,false,false,false,false,false,false,false,false,false,false,false,
815	false,false,false,false,false,false,false,false,false,false,false,false,false,false,
816	false,false,false,false,false,false,false,false,false,false,false,false,false,false};
817	static G4double VZ0[nEl][iNR];
818	// ******* S=1 vectors ***********
819	static G4bool iNin1[nEl]={false,false,false,false,false,false,false,
820	false,false,false,false,false,false,false,false,false,false,false,false,false,false,
821	false,false,false,false,false,false,false,false,false,false,false,false,false,false,
822	false,false,false,false,false,false,false,false,false,false,false,false,false,false,
823	false,false,false,false,false,false,false,false,false,false,false,false,false,false,
824	false,false,false,false,false,false,false,false,false,false,false,false,false,false,
825	false,false,false,false,false,false,false,false,false,false,false,false,false,false,
826	false,false,false,false,false,false,false,false,false,false,false,false,false,false};
827	static G4double VZ1[nEl][iNR];
828	// ******* S=2 vectors ***********
829	static G4bool iNin2[nEl]={false,false,false,false,false,false,false,
830	false,false,false,false,false,false,false,false,false,false,false,false,false,false,
831	false,false,false,false,false,false,false,false,false,false,false,false,false,false,
832	false,false,false,false,false,false,false,false,false,false,false,false,false,false,
833	false,false,false,false,false,false,false,false,false,false,false,false,false,false,
834	false,false,false,false,false,false,false,false,false,false,false,false,false,false,
835	false,false,false,false,false,false,false,false,false,false,false,false,false,false,
836	false,false,false,false,false,false,false,false,false,false,false,false,false,false};
837	static G4double VZ2[nEl][iNR];
838	// ******* S=3 vectors ***********
839	static G4bool iNin3[nEl]={false,false,false,false,false,false,false,
840	false,false,false,false,false,false,false,false,false,false,false,false,false,false,
841	false,false,false,false,false,false,false,false,false,false,false,false,false,false,
842	false,false,false,false,false,false,false,false,false,false,false,false,false,false,
843	false,false,false,false,false,false,false,false,false,false,false,false,false,false,
844	false,false,false,false,false,false,false,false,false,false,false,false,false,false,
845	false,false,false,false,false,false,false,false,false,false,false,false,false,false,
846	false,false,false,false,false,false,false,false,false,false,false,false,false,false};
847	static G4double VZ3[nEl][iNR];
848	// ******* S=2 vectors ***********
849	static G4bool iNin4[nEl]={false,false,false,false,false,false,false,
850	false,false,false,false,false,false,false,false,false,false,false,false,false,false,
851	false,false,false,false,false,false,false,false,false,false,false,false,false,false,
852	false,false,false,false,false,false,false,false,false,false,false,false,false,false,
853	false,false,false,false,false,false,false,false,false,false,false,false,false,false,
854	false,false,false,false,false,false,false,false,false,false,false,false,false,false,
855	false,false,false,false,false,false,false,false,false,false,false,false,false,false,
856	false,false,false,false,false,false,false,false,false,false,false,false,false,false};
857	static G4double VZ4[nEl][iNR];
858	//
859	#ifdef qdebug
860	static G4int Smin=-1; // Dynamic Associated memory is implemented for S>-2 nuclei
861	static G4int Smax= 2; // Dynamic Associated memory is implemented for S< 3 nuclei
862	static G4int NZmin= 0; // Dynamic Associated memory is implemented for Z>-1 nuclei
863	static G4int NNmin= 0; // Dynamic Associated memory is implemented for N>-1 nuclei
864	static G4int NZS1max= 0; // Dynamic Associated memory is implemented for S<3, Z=-1, N<2
865	static G4int NNS1max= 0; // Dynamic Associated memory is implemented for S<3, Z=-1, N<2
866	#endif
867	// -------------------------------------------------------------------------------------
868	G4double rm=0.;
869	G4int nz=n+z;
870	G4int zns=nz+s;
871	if(nz+s<0)
872	{
873	z=-z;
874	n=-n;
875	s=-s;
876	nz=-nz;
877	zns=-zns;
878	}
879	if(z<0)
880	{
881	if(z==-1)
882	{
883	if(!s)
884	{
885	if(n==1) return mPiC; // pi-
886	else return mPiC+(n-1)mNeut; // Delta- + (N-1)n
887	}
888	else if(s==1) // Strange negative hadron
889	{
890	if(!n) return mKM; // K-
891	else if(n==1) return mSiM; // Sigma-
892	else if(n==2) return mSmN ; // Sigma- + n DiBaryon
893	else if(n==3) return mSmNN; // Sigma- +2n TriBaryon
894	else return mSigM+mNeut(n-1); // Sigma- + (N-1)n
895	}
896	else if(s==2) // --> Double-strange negative hadrons
897	{
898	if(!n) return mKsM; // Ksi-
899	else if(n==1) return mKsiM+mNeut; // Ksi- + n
900	else if(n==2) return mKsiM+mNeut+mNeut; // Ksi- + 2n
901	else return mKsiM+mNeutn; // Ksi- + Zn
902	}
903	else if(s==-2)
904	{
905	if (nz==2) return mDiKZ+mPiC; // 2K0 + Pi-
906	else return mDiKZ+mPiC+(nz-2)*mProt;
907	}
908	else if(s==3) // --> Triple-strange negative hadrons
909	{
910	if (n==-1) return mOmM; // Triple-strange Omega minus
911	else if(!n ) return mOmN; // Triple-strange (Omega-) + n DiBaryon
912	else if(n==-2) return mDiKZ+mKM; // Triple-strange K- + 2*K0
913	else return mOmeg+mNeut*(n+2);
914	}
915	else if(s==4)
916	{
917	if(n==-2) return mOmeg+mKM; // Omega- + K-
918	else if(n==-1) return mOmeg+mLamb;// Omega- + Lambda
919	else return mOmeg+mLamb+(n+1)mNeut; // Omega- + Lambda + (n+1)Neutrons
920	}
921	else if(!n) return mOmeg+(s-2)*mLamb; // Multy-Lambda + Omega minus
922	else
923	{
924	#ifdef qdebug
925	if(s>NZS1max)
926	{
927	NZS1max=s;
928	G4cout<<">>>>>>>>>>>>>G4QPDGCode::GetMass: Z=-1, S="<<s<<">2 with N="<<n<<G4endl;
929	}
930	#endif
931	return CalculateNuclMass(z,n,s);
932	}
933	}
934	else if(!s)
935	{
936	if(!nz)
937	{
938	if(n==2) return mDiPi;
939	else return mDiPi+(n-2)*mPiC;
940	}
941	else return mNeutnz-zmPiC+anb;
942	}
943	else if(!zns) // !!! s=0 is treated above !!
944	{
945	if(s>0) return anb+smKM+nmPiC; // sK- + nPi-
946	else return anb-smKZ-zmPiC; // (-s)aK0 + (-z)Pi-
947	}
948	else if(s==1)
949	{
950	if(!nz)
951	{
952	if(n==2) return mSmPi;
953	else return mSmPi+(n-2)*mPiC;
954	}
955	else return mSigM+nzmNeut-(z+1)mPiC;
956	}
957	else if(s==-1) return mKZa-zmPiC+(nz-1)mNeut; // aK0+(nz-1)n+(-z)*Pi-
958	else if(s==2)
959	{
960	if (nz==-1) return mKsiM+n*mPiC;
961	else if(!nz) return mKsiM+mNeut-(z+1)*mPiC;
962	else return mKsiM+(nz+1)mNeut-(z+1)mPiC;
963	}
964	else if(s==-2) return mDiKZ-zmPiC+(nz-2)mNeut;
965	else if(s==3)
966	{
967	if (nz==-2) return mOmeg+(n+1)mPiC; // Omega- + (n+1) Pi-
968	else if(nz==-1) return mOmeg+mNeut+nmPiC; // Omega- + n Pi-
969	else if(!nz) return mOmeg+mDiNt+(n-1)mPiC; // Omega- + 2N + (n-1)Pi-
970	else return mOmeg+(nz+2)mProt-(z+1)mPiC;
971	}
972	else if(s<-2) return anb-smKZ-zmPiC+(nz+s)*mNeut;
973	else if(s==4)
974	{
975	if (nz==-3) return mOmeg+mKM+(n+1)mPiC; // Om- + K- + (n+1)Pi-
976	else if(nz==-2) return mOmeg+mSigM+nmPiC; // Om- + Sig- + nPi-
977	else if(nz==-1) return mOmeg+mSigM+mNeut+(n-1)mPiC;//Om- + Sig- +N +(n-1)Pi-
978	else if(!nz) return mOmeg+mSigM+mDiNt+(n-2)mPiC;//Om- +Sig- +2N +(n-2)Pi-
979	else return mOmeg+mSigM+(nz+2)mDiNt-(z+2)mPiC;//+(nz+2)N-(z+2)Pi-
980	}
981	// s=5: 2K-, Ksi-; s=6: 3K-, Omega-; s>6 adds K- and Sigma- instead of Protons
982	else // !!All s<0 are done and s>4 can be easy extended if appear!!
983	{
984	#ifdef qdebug
985	if(z<NZmin)
986	{
987	NZmin=z;
988	G4cout<<">>>>>>>>>G4QPDGCode::GetMass: Z="<<z<<"<-1 with N="<<n<<", S="<<s<<G4endl;
989	}
990	#endif
991	return CalculateNuclMass(z,n,s);
992	}
993	}
994	else if(n<0)
995	{
996	if(n==-1)
997	{
998	if(!s)
999	{
1000	if(z==1) return mPiC; // pi+
1001	else return mPiC+(z-1)mProt; // Delta++ + (Z-1)p
1002	}
1003	else if(s==1) // --> Strange neutral hadrons
1004	{
1005	if(!z) return mKZ; // K0
1006	else if(z==1) return mSiP; // Sigma+
1007	else if(z==2) return mSpP ; // Sigma+ + p DiBaryon
1008	else if(z==3) return mSpPP; // Sigma+ +2p TriBaryon
1009	else return mSigP+mProt(z-1); // Sigma+ + (Z-1)p
1010	}
1011	else if(s==2) // --> Double-strange negative hadrons
1012	{
1013	if(!z) return mKsZ; // Ksi0
1014	else if(z==1) return mKsiZ+mProt; // Ksi- + p
1015	else if(z==2) return mKsiZ+mProt+mProt; // Ksi- + 2p
1016	else return mKsiZ+mProtz; // Ksi- + Zp
1017	}
1018	else if(s==-2)
1019	{
1020	if (nz==2) return mDiKM+mPiC; // 2K+ + Pi+
1021	else return mDiKM+mPiC+(nz-2)*mProt;
1022	}
1023	else if(s==3)
1024	{
1025	if(z==1) return mOmeg+mDiPr;
1026	else return mOmeg+(z+1)*mProt;
1027	}
1028	else if(s==4) return mOmeg+mLamb+(z+1)*mProt;
1029	else if(!z) return mKZa+(s-1)*mLamb; // Multy-Lambda + K0
1030	else
1031	{
1032	#ifdef qdebug
1033	if(s>NNS1max)
1034	{
1035	NNS1max=s;
1036	G4cout<<">>>>>>>>>>>>>G4QPDGCode::GetMass: N=-1, S="<<s<<">2 with Z="<<z<<G4endl;
1037	}
1038	#endif
1039	return CalculateNuclMass(z,n,s);
1040	}
1041	}
1042	else if(!s)
1043	{
1044	if(!nz)
1045	{
1046	if(z==2) return mDiPi;
1047	else return mDiPi+(z-2)*mPiC;
1048	}
1049	else return mProtnz-nmPiC+anb;
1050	}
1051	else if(!zns) // !!! s=0 is treated above !!
1052	{
1053	if(s>0) return anb+smKZ+zmPiC; // sK0 + nPi+
1054	else return anb-smKM-nmPiC; // (-s)aK+ + (-n)Pi+
1055	}
1056	else if(s==1)
1057	{
1058	if(!nz)
1059	{
1060	if(z==2) return mSpPi;
1061	else return mSpPi+(z-2)*mPiC;
1062	}
1063	else return mSigP+nzmProt-(n+1)mPiC;
1064	}
1065	else if(s==-1) return mKMa-nmPiC+(nz-1)mProt; // K+ + (nz-1)P + (-n)Pi+
1066	else if(s==2)
1067	{
1068	if (nz==-1) return mKsiZ+z*mPiC;
1069	else if(!nz) return mKsiZ+mProt-(n+1)*mPiC;
1070	else return mKsiZ+(nz+1)mProt-(n+1)mPiC;
1071	}
1072	else if(s==-2) return mDiKM-nmPiC+(nz-2)mProt;
1073	else if(s==3)
1074	{
1075	if (nz==-2) return mOmeg+(z+1)mPiC; // Omega + (z+1)Pi+
1076	else if(nz==-1) return mOmeg+mProt+zmPiC; // Omega- + P +zPi+
1077	else if(!nz) return mOmeg+mDiPr+(z-1)mPiC; // Omega- + 2P + (z-1) Pi+
1078	else return mOmeg+(nz+2)mProt-(n+1)mPiC;
1079	}
1080	else if(s<-2) return anb-smKM-nmPiC+(nz+s)*mProt;
1081	else if(s==4)
1082	{
1083	if (nz==-3) return mOmeg+mKZ+(z+1)mPiC; // Om- + K0 + (z+1)Pi+
1084	else if(nz==-2) return mOmeg+mSigP+zmPiC; // Om- + Sig+ + zPi+
1085	else if(nz==-1) return mOmeg+mSigP+mProt+(z-1)mPiC;// Om- +Sig+ +P +(z-1)Pi+
1086	else if(!nz) return mOmeg+mSigP+mDiPr+(z-2)mPiC;//Om- +Sig+ +2P +(z-2)Pi+
1087	else return mOmeg+mSigP+(nz+2)mProt-(n+2)mPiC;//+(nz+2)P-(n+2)Pi+
1088	}
1089	// s=5: 2KZ, Ksi0; s=6: 3KZ, Omega-; s>6 adds K0 and Sigma+ instead of Protons
1090	else // !!All s<0 are done and s>4 can be easy extended if appear!!
1091	{
1092	#ifdef qdebug
1093	if(n<NNmin)
1094	{
1095	NNmin=n;
1096	G4cout<<">>>>>>>>>G4QPDGCode::GetMass: N="<<n<<"<-1 with Z="<<z<<", S="<<s<<G4endl;
1097	}
1098	#endif
1099	return CalculateNuclMass(z,n,s);
1100	}
1101	}
1102	//return CalculateNuclMass(z,n,s); // @@ This is just to compare the calculation time @@
1103	if(!s) // **************> S=0 nucleus
1104	{
1105	if(nz==256) return 256000.;
1106	G4int Nmin=iNF[z]; // Minimun N(Z) for the Dynamic Associative Memory (DAM)
1107	if(!iNin0[z]) // ====> This element is already initialized
1108	{
1109	#ifdef idebug
1110	G4cout<<"**>G4QPDGCode::GetMass:Z="<<z<<", S=0 is initialized. F="<<iNin0[z]<<G4endl;
1111	#endif
1112	G4int iNfin=iNL[z];
1113	if(iNfin>iNR) iNfin=iNR;
1114	for (G4int in=0; in<iNfin; in++) VZ0[z][in] = CalculateNuclMass(z,in+Nmin,s);
1115	iNin0[z]=true;
1116	}
1117	G4int dNn=n-Nmin;
1118	if(dNn<0) // --> The minimum N must be reduced
1119	{
1120	#ifdef qdebug
1121	if(n<iNmin[z])
1122	{
1123	G4cout<<">>>>G4QPDGCode::GetMass:Z="<<z<<", S=0 with N="<<n<<"<"<<iNmin[z]<<G4endl;
1124	iNmin[z]=n;
1125	}
1126	#endif
1127	return CalculateNuclMass(z,n,s);
1128	}
1129	else if(dNn<iNL[z]) return VZ0[z][dNn]; // Found in DAM
1130	else // --> The maximum N must be increased
1131	{
1132	#ifdef qdebug
1133	if(dNn>iNmax)
1134	{
1135	G4cout<<"**>>G4QPDGCode::GetMass:Z="<<z<<", S=0 with dN="<<dNn<<">"<<iNmax<<G4endl;
1136	iNmax=dNn;
1137	}
1138	if(dNn>iNran[z])
1139	{
1140	G4cout<<">G4QPDGCode::GetMass:Z="<<z<<", S=0 with dN="<<dNn<<">"<<iNran[z]<<G4endl;
1141	iNran[z]=dNn;
1142	}
1143	#endif
1144	return CalculateNuclMass(z,n,s);
1145	}
1146	}
1147	else if(s==1) // ******************> S=1 nucleus
1148	{
1149
1150	G4int Nmin=iNF[z]; // Minimun N(Z) for the Dynamic Associative Memory (DAM)
1151	if(!iNin1[z]) // ====> This element is already initialized
1152	{
1153	#ifdef idebug
1154	G4cout<<"**>G4QPDGCode::GetMass:Z="<<z<<", S=1 is initialized. F="<<iNin1[z]<<G4endl;
1155	#endif
1156	G4int iNfin=iNL[z];
1157	if(iNfin>iNR) iNfin=iNR;
1158	for (G4int in=0; in<iNfin; in++) VZ1[z][in] = CalculateNuclMass(z,in+Nmin,s);
1159	iNin1[z]=true;
1160	}
1161	G4int dNn=n-Nmin;
1162	if(dNn<0) // --> The minimum N must be reduced
1163	{
1164	#ifdef qdebug
1165	if(n<iNmin[z])
1166	{
1167	G4cout<<">>>>G4QPDGCode::GetMass:Z="<<z<<", S=1 with N="<<n<<"<"<<iNmin[z]<<G4endl;
1168	iNmin[z]=n;
1169	}
1170	#endif
1171	return CalculateNuclMass(z,n,s);
1172	}
1173	else if(dNn<iNL[z]) return VZ1[z][dNn]; // Found in DAM
1174	else // --> The maximum N must be increased
1175	{
1176	#ifdef qdebug
1177	if(dNn>iNmax)
1178	{
1179	G4cout<<"**>>G4QPDGCode::GetMass:Z="<<z<<", S=1 with dN="<<dNn<<">"<<iNmax<<G4endl;
1180	iNmax=dNn;
1181	}
1182	if(dNn>iNran[z])
1183	{
1184	G4cout<<">G4QPDGCode::GetMass:Z="<<z<<", S=1 with dN="<<dNn<<">"<<iNran[z]<<G4endl;
1185	iNran[z]=dNn;
1186	}
1187	#endif
1188	return CalculateNuclMass(z,n,s);
1189	}
1190	}
1191	else if(s==-1) // ******************> S=-1 nucleus
1192	{
1193	G4int Nmin=iNF[z]; // Minimun N(Z) for the Dynamic Associative Memory (DAM)
1194	if(!iNin9[z]) // ====> This element is already initialized
1195	{
1196	#ifdef idebug
1197	G4cout<<"*>G4QPDGCode::GetMass:Z="<<z<<", S=-1 is initialized. F="<<iNin9[z]<<G4endl;
1198	#endif
1199	G4int iNfin=iNL[z];
1200	if(iNfin>iNR) iNfin=iNR;
1201	for (G4int in=0; in<iNfin; in++) VZ9[z][in] = CalculateNuclMass(z,in+Nmin,s);
1202	iNin9[z]=true;
1203	}
1204	G4int dNn=n-Nmin;
1205	if(dNn<0) // --> The minimum N must be reduced
1206	{
1207	#ifdef qdebug
1208	if(n<iNmin[z])
1209	{
1210	G4cout<<">>>G4QPDGCode::GetMass:Z="<<z<<" ,S=-1 with N="<<n<<"<"<<iNmin[z]<<G4endl;
1211	iNmin[z]=n;
1212	}
1213	#endif
1214	return CalculateNuclMass(z,n,s);
1215	}
1216	else if(dNn<iNL[z]) return VZ9[z][dNn]; // Found in DAM
1217	else // --> The maximum N must be increased
1218	{
1219	#ifdef qdebug
1220	if(dNn>iNmax)
1221	{
1222	G4cout<<"**>G4QPDGCode::GetMass:Z="<<z<<", S=-1 with dN="<<dNn<<">"<<iNmax<<G4endl;
1223	iNmax=dNn;
1224	}
1225	if(dNn>iNran[z])
1226	{
1227	G4cout<<"G4QPDGCode::GetMass:Z="<<z<<", S=-1 with dN="<<dNn<<">"<<iNran[z]<<G4endl;
1228	iNran[z]=dNn;
1229	}
1230	#endif
1231	return CalculateNuclMass(z,n,s);
1232	}
1233	}
1234	else if(s==2) // ******************> S=2 nucleus
1235	{
1236	G4int Nmin=iNF[z]; // Minimun N(Z) for the Dynamic Associative Memory (DAM)
1237	if(!iNin2[z]) // ====> This element is already initialized
1238	{
1239	#ifdef idebug
1240	G4cout<<"**>G4QPDGCode::GetMass:Z="<<z<<", S=2 is initialized. F="<<iNin2[z]<<G4endl;
1241	#endif
1242	G4int iNfin=iNL[z];
1243	if(iNfin>iNR) iNfin=iNR;
1244	for (G4int in=0; in<iNfin; in++) VZ2[z][in] = CalculateNuclMass(z,in+Nmin,s);
1245	iNin2[z]=true;
1246	}
1247	G4int dNn=n-Nmin;
1248	if(dNn<0) // --> The minimum N must be reduced
1249	{
1250	#ifdef qdebug
1251	if(n<iNmin[z])
1252	{
1253	G4cout<<">>>>G4QPDGCode::GetMass:Z="<<z<<", S=2 with N="<<n<<"<"<<iNmin[z]<<G4endl;
1254	iNmin[z]=n;
1255	}
1256	#endif
1257	return CalculateNuclMass(z,n,s);
1258	}
1259	else if(dNn<iNL[z]) return VZ2[z][dNn]; // Found in DAM
1260	else // --> The maximum N must be increased
1261	{
1262	#ifdef qdebug
1263	if(dNn>iNmax)
1264	{
1265	G4cout<<"**>>G4QPDGCode::GetMass:Z="<<z<<", S=2 with dN="<<dNn<<">"<<iNmax<<G4endl;
1266	iNmax=dNn;
1267	}
1268	if(dNn>iNran[z])
1269	{
1270	G4cout<<">G4QPDGCode::GetMass:Z="<<z<<", S=2 with dN="<<dNn<<">"<<iNran[z]<<G4endl;
1271	iNran[z]=dNn;
1272	}
1273	#endif
1274	return CalculateNuclMass(z,n,s);
1275	}
1276	}
1277	else if(s==-2) // ******************> S=-2 nucleus
1278	{
1279	G4int Nmin=iNF[z]; // Minimun N(Z) for the Dynamic Associative Memory (DAM)
1280	if(!iNin8[z]) // ====> This element is already initialized
1281	{
1282	#ifdef idebug
1283	G4cout<<"*>G4QPDGCode::GetMass:Z="<<z<<", S=-2 is initialized. F="<<iNin8[z]<<G4endl;
1284	#endif
1285	G4int iNfin=iNL[z];
1286	if(iNfin>iNR) iNfin=iNR;
1287	for (G4int in=0; in<iNfin; in++) VZ8[z][in] = CalculateNuclMass(z,in+Nmin,s);
1288	iNin8[z]=true;
1289	}
1290	G4int dNn=n-Nmin;
1291	if(dNn<0) // --> The minimum N must be reduced
1292	{
1293	#ifdef qdebug
1294	if(n<iNmin[z])
1295	{
1296	G4cout<<">>>G4QPDGCode::GetMass:Z="<<z<<", S=-2 with N="<<n<<"<"<<iNmin[z]<<G4endl;
1297	iNmin[z]=n;
1298	}
1299	#endif
1300	return CalculateNuclMass(z,n,s);
1301	}
1302	else if(dNn<iNL[z]) return VZ8[z][dNn]; // Found in DAM
1303	else // --> The maximum N must be increased
1304	{
1305	#ifdef qdebug
1306	if(dNn>iNmax)
1307	{
1308	G4cout<<"**>G4QPDGCode::GetMass:Z="<<z<<", S=-2 with dN="<<dNn<<">"<<iNmax<<G4endl;
1309	iNmax=dNn;
1310	}
1311	if(dNn>iNran[z])
1312	{
1313	G4cout<<"G4QPDGCode::GetMass:Z="<<z<<", S=-2 with dN="<<dNn<<">"<<iNran[z]<<G4endl;
1314	iNran[z]=dNn;
1315	}
1316	#endif
1317	return CalculateNuclMass(z,n,s);
1318	}
1319	}
1320	else if(s==-3) // ******************> S=-3 nucleus
1321	{
1322	G4int Nmin=iNF[z]; // Minimun N(Z) for the Dynamic Associative Memory (DAM)
1323	if(!iNin7[z]) // ====> This element is already initialized
1324	{
1325	#ifdef idebug
1326	G4cout<<"*>G4QPDGCode::GetMass:Z="<<z<<", S=-3 is initialized. F="<<iNin7[z]<<G4endl;
1327	#endif
1328	G4int iNfin=iNL[z];
1329	if(iNfin>iNR) iNfin=iNR;
1330	for (G4int in=0; in<iNfin; in++) VZ7[z][in] = CalculateNuclMass(z,in+Nmin,s);
1331	iNin7[z]=true;
1332	}
1333	G4int dNn=n-Nmin;
1334	if(dNn<0) // --> The minimum N must be reduced
1335	{
1336	#ifdef qdebug
1337	if(n<iNmin[z])
1338	{
1339	G4cout<<">>>G4QPDGCode::GetMass:Z="<<z<<", S=-3 with N="<<n<<"<"<<iNmin[z]<<G4endl;
1340	iNmin[z]=n;
1341	}
1342	#endif
1343	return CalculateNuclMass(z,n,s);
1344	}
1345	else if(dNn<iNL[z]) return VZ7[z][dNn]; // Found in DAM
1346	else // --> The maximum N must be increased
1347	{
1348	#ifdef qdebug
1349	if(dNn>iNmax)
1350	{
1351	G4cout<<"**>G4QPDGCode::GetMass:Z="<<z<<", S=-3 with dN="<<dNn<<">"<<iNmax<<G4endl;
1352	iNmax=dNn;
1353	}
1354	if(dNn>iNran[z])
1355	{
1356	G4cout<<"G4QPDGCode::GetMass:Z="<<z<<", S=-3 with dN="<<dNn<<">"<<iNran[z]<<G4endl;
1357	iNran[z]=dNn;
1358	}
1359	#endif
1360	return CalculateNuclMass(z,n,s);
1361	}
1362	}
1363	else if(s==3) // ******************> S=3 nucleus
1364	{
1365	G4int Nmin=iNF[z]; // Minimun N(Z) for the Dynamic Associative Memory (DAM)
1366	if(!iNin3[z]) // ====> This element is already initialized
1367	{
1368	#ifdef idebug
1369	G4cout<<"**>G4QPDGCode::GetMass:Z="<<z<<", S=3 is initialized. F="<<iNin3[z]<<G4endl;
1370	#endif
1371	G4int iNfin=iNL[z];
1372	if(iNfin>iNR) iNfin=iNR;
1373	for (G4int in=0; in<iNfin; in++) VZ3[z][in] = CalculateNuclMass(z,in+Nmin,s);
1374	iNin3[z]=true;
1375	}
1376	G4int dNn=n-Nmin;
1377	if(dNn<0) // --> The minimum N must be reduced
1378	{
1379	#ifdef qdebug
1380	if(n<iNmin[z])
1381	{
1382	G4cout<<">>>>G4QPDGCode::GetMass:Z="<<z<<", S=3 with N="<<n<<"<"<<iNmin[z]<<G4endl;
1383	iNmin[z]=n;
1384	}
1385	#endif
1386	return CalculateNuclMass(z,n,s);
1387	}
1388	else if(dNn<iNL[z]) return VZ3[z][dNn]; // Found in DAM
1389	else // --> The maximum N must be increased
1390	{
1391	#ifdef qdebug
1392	if(dNn>iNmax)
1393	{
1394	G4cout<<"**>>G4QPDGCode::GetMass:Z="<<z<<", S=3 with dN="<<dNn<<">"<<iNmax<<G4endl;
1395	iNmax=dNn;
1396	}
1397	if(dNn>iNran[z])
1398	{
1399	G4cout<<">G4QPDGCode::GetMass:Z="<<z<<", S=3 with dN="<<dNn<<">"<<iNran[z]<<G4endl;
1400	iNran[z]=dNn;
1401	}
1402	#endif
1403	return CalculateNuclMass(z,n,s);
1404	}
1405	}
1406	else if(s==-4) // ******************> S=-4 nucleus
1407	{
1408	G4int Nmin=iNF[z]; // Minimun N(Z) for the Dynamic Associative Memory (DAM)
1409	if(!iNin6[z]) // ====> This element is already initialized
1410	{
1411	#ifdef idebug
1412	G4cout<<"*>G4QPDGCode::GetMass:Z="<<z<<", S=-4 is initialized. F="<<iNin6[z]<<G4endl;
1413	#endif
1414	G4int iNfin=iNL[z];
1415	if(iNfin>iNR) iNfin=iNR;
1416	for (G4int in=0; in<iNfin; in++) VZ6[z][in] = CalculateNuclMass(z,in+Nmin,s);
1417	iNin6[z]=true;
1418	}
1419	G4int dNn=n-Nmin;
1420	if(dNn<0) // --> The minimum N must be reduced
1421	{
1422	#ifdef qdebug
1423	if(n<iNmin[z])
1424	{
1425	G4cout<<">>>G4QPDGCode::GetMass:Z="<<z<<", S=-4 with N="<<n<<"<"<<iNmin[z]<<G4endl;
1426	iNmin[z]=n;
1427	}
1428	#endif
1429	return CalculateNuclMass(z,n,s);
1430	}
1431	else if(dNn<iNL[z]) return VZ6[z][dNn]; // Found in DAM
1432	else // --> The maximum N must be increased
1433	{
1434	#ifdef qdebug
1435	if(dNn>iNmax)
1436	{
1437	G4cout<<"**>G4QPDGCode::GetMass:Z="<<z<<", S=-4 with dN="<<dNn<<">"<<iNmax<<G4endl;
1438	iNmax=dNn;
1439	}
1440	if(dNn>iNran[z])
1441	{
1442	G4cout<<"G4QPDGCode::GetMass:Z="<<z<<", S=-4 with dN="<<dNn<<">"<<iNran[z]<<G4endl;
1443	iNran[z]=dNn;
1444	}
1445	#endif
1446	return CalculateNuclMass(z,n,s);
1447	}
1448	}
1449	else if(s==4) // ******************> S=4 nucleus
1450	{
1451	G4int Nmin=iNF[z]; // Minimun N(Z) for the Dynamic Associative Memory (DAM)
1452	if(!iNin4[z]) // ====> This element is already initialized
1453	{
1454	#ifdef idebug
1455	G4cout<<"*>G4QPDGCode::GetMass:Z="<<z<<", S=4 is initialized. F="<<iNin4[z]<<G4endl;
1456	#endif
1457	G4int iNfin=iNL[z];
1458	if(iNfin>iNR) iNfin=iNR;
1459	for (G4int in=0; in<iNfin; in++) VZ4[z][in] = CalculateNuclMass(z,in+Nmin,s);
1460	iNin4[z]=true;
1461	}
1462	G4int dNn=n-Nmin;
1463	if(dNn<0) // --> The minimum N must be reduced
1464	{
1465	#ifdef qdebug
1466	if(n<iNmin[z])
1467	{
1468	G4cout<<">>>>G4QPDGCode::GetMass:Z="<<z<<", S=4 with N="<<n<<"<"<<iNmin[z]<<G4endl;
1469	iNmin[z]=n;
1470	}
1471	#endif
1472	return CalculateNuclMass(z,n,s);
1473	}
1474	else if(dNn<iNL[z]) return VZ4[z][dNn]; // Found in DAM
1475	else // --> The maximum N must be increased
1476	{
1477	#ifdef qdebug
1478	if(dNn>iNmax)
1479	{
1480	G4cout<<"**>>G4QPDGCode::GetMass:Z="<<z<<", S=4 with dN="<<dNn<<">"<<iNmax<<G4endl;
1481	iNmax=dNn;
1482	}
1483	if(dNn>iNran[z])
1484	{
1485	G4cout<<">G4QPDGCode::GetMass:Z="<<z<<", S=4 with dN="<<dNn<<">"<<iNran[z]<<G4endl;
1486	iNran[z]=dNn;
1487	}
1488	#endif
1489	return CalculateNuclMass(z,n,s);
1490	}
1491	}
1492	else
1493	{
1494	#ifdef qdebug
1495	if(s<Smin \|\| s>Smax)
1496	{
1497	if(s<Smin) Smin=s;
1498	if(s>Smax) Smax=s;
1499	G4cout<<">>G4QPDGCode::GetMass:Z="<<z<<" with S="<<s<<",N="<<n<<" (Improve)"<<G4endl;
1500	}
1501	#endif
1502	rm=CalculateNuclMass(z,n,s);
1503	}
1504	#ifdef debug
1505	G4cout<<"G4QPDGCode::GetMass:GetNuclMass="<<rm<<",Z="<<z<<",N="<<n<<",S="<<s<<G4endl;
1506	#endif
1507	return rm;
1508	}
1509
1510	// Calculate a nuclear mass for Z (a#of protons), N (a#of neutrons), & S (a#of lambdas)
1511	G4double G4QPDGCode::CalculateNuclMass(G4int z, G4int n, G4int s)
1512	// =========================================================
1513	{
1514	static const G4double mP = QHaM(21); // Proton
1515	static const G4double mN = QHaM(20); // Neutron
1516	static const G4double mL = QHaM(22); // Lambda
1517	static const G4double mD = G4Deuteron::Deuteron()->GetPDGMass(); // Deuteron (H-2)
1518	static const G4double mT = G4Triton::Triton()->GetPDGMass(); // Triton (H-3)
1519	static const G4double mHe3= G4He3::He3()->GetPDGMass(); // Hetrium (He-3)
1520	static const G4double mAl = G4Alpha::Alpha()->GetPDGMass(); // Alpha (He-4)
1521	static const G4double dmP = mP+mP; // DiProton
1522	static const G4double dmN = mN+mN; // DiNeutron
1523	static const G4double dmL = mL+mL; // DiLambda
1524	static const G4double dLN = mL+mN; // LambdaNeutron
1525	static const G4double dLP = mL+mP; // LambdaProton
1526	static const G4double mSm = QHaM(23); // Sigma-
1527	static const G4double mSp = QHaM(25); // Sigma+
1528	static const G4double dSP = mSp+mP; // ProtonSigma+
1529	static const G4double dSN = mSm+mN; // NeutronSigma-
1530	static const G4double dnS = dSN+mN; // 2NeutronsSigma-
1531	static const G4double dpS = dSP+mP; // 2ProtonsSigma+
1532	static const G4double mXm = QHaM(26); // Ksi-
1533	static const G4double mXz = QHaM(27); // Ksi0
1534	static const G4double mOm = QHaM(44); // Omega-
1535	static const G4double dXN = mXm+mN; // NeutronKsi-
1536	static const G4double dXP = mXz+mP; // ProtonKsi0
1537	static const G4double dOP = mOm+mP; // ProtonOmega-
1538	static const G4double dON = mOm+mN; // NeutronOmega-
1539	static const G4double mK = QHaM(18); // Charged Kaon
1540	static const G4double mK0 = QHaM(17); // Neutral Kaon
1541	static const G4double mPi = QHaM(15); // Charged Pion
1542	//////////static const G4double mPi0= QHaM(14);
1543	//static const G4int nSh = 164;
1544	//static G4double sh[nSh] = {0., // Fake element for C++ N=Z=0
1545	// -4.315548, 2.435504, -1.170501, 3.950887, 5.425238,
1546	// 13.342524, 15.547586, 22.583129, 23.983480, 30.561036,
1547	// 33.761971, 41.471027, 45.532156, 53.835880, 58.495514,
1548	// 65.693445, 69.903344, 76.899581, 81.329361, 88.979438,
1549	// 92.908703, 100.316636, 105.013393, 113.081686, 118.622601,
1550	// 126.979113, 132.714435, 141.413182, 146.433488, 153.746754,
1551	// 158.665225, 165.988967, 170.952395, 178.473011, 183.471531,
1552	// 191.231310, 196.504414, 204.617158, 210.251108, 218.373984,
1553	// 223.969281, 232.168660, 237.925619, 246.400505, 252.392471,
1554	// 260.938644, 267.191321, 276.107788, 282.722682, 291.881502,
1555	// 296.998590, 304.236025, 309.562296, 316.928655, 322.240263,
1556	// 329.927236, 335.480630, 343.233705, 348.923475, 356.911659,
1557	// 362.785757, 370.920926, 376.929998, 385.130316, 391.197741,
1558	// 399.451554, 405.679971, 414.101869, 420.346260, 428.832412,
1559	// 435.067445, 443.526983, 449.880034, 458.348602, 464.822352,
1560	// 473.313779, 479.744524, 488.320887, 495.025069, 503.841579,
1561	// 510.716379, 519.451976, 525.036156, 532.388151, 537.899017,
1562	// 545.252264, 550.802469, 558.402181, 564.101100, 571.963429,
1563	// 577.980340, 586.063802, 592.451334, 600.518525, 606.832027,
1564	// 614.831626, 621.205330, 629.237413, 635.489106, 643.434167,
1565	// 649.691284, 657.516479, 663.812101, 671.715021, 678.061128,
1566	// 686.002970, 692.343712, 700.360477, 706.624091, 714.617848,
1567	// 721.100390, 729.294717, 735.887170, 744.216084, 751.017094,
1568	// 759.551764, 766.377807, 775.080204, 781.965673, 790.552795,
1569	// 797.572494, 806.088030, 813.158751, 821.655631, 828.867137,
1570	// 836.860955, 842.183292, 849.195302, 854.731798, 861.898839,
1571	// 867.783606, 875.313342, 881.443441, 889.189065, 895.680189,
1572	// 903.679729, 910.368085, 918.579876, 925.543547, 933.790028,
1573	// 940.811396, 949.122548, 956.170201, 964.466810, 971.516490,
1574	// 979.766905, 986.844659, 995.113552,1002.212760,1010.418770,
1575	// 1018.302560,1025.781870,1033.263560,1040.747880,1048.234460,
1576	// 1055.723430,1063.214780,1070.708750,1078.204870,1085.703370,
1577	// 1093.204260,1100.707530,1108.213070};
1578	//static const G4double b1=8.09748564; // MeV
1579	//static const G4double b2=-0.76277387;
1580	//static const G4double b3=83.487332; // MeV
1581	//static const G4double b4=0.090578206;// 2*b4
1582	//static const G4double b5=0.676377211;// MeV
1583	//static const G4double b6=5.55231981; // MeV
1584	static const G4double b7=25.; // MeV (Lambda binding energy predexponent)
1585	// --- even-odd difference is 3.7(MeV)/X
1586	// --- S(X>151)=-57.56-5.542X*1.05
1587	static const G4double b8=10.5; // (Lambda binding energy exponent)
1588	//static const G4double b9=-1./3.;
1589	static const G4double a2=0.13; // LambdaBindingEnergy for deutron+LambdaSystem(MeV)
1590	static const G4double a3=2.2; // LambdaBindingEnergy for (t/He3)+LambdaSystem(MeV)
1591	static const G4double um=931.49432; // Umified atomic mass unit (MeV)
1592	//static const G4double me =0.511; // electron mass (MeV) :: n:8.071, p:7.289
1593	static const G4double eps =0.0001; // security addition for multybaryons
1594	//static G4double c[9][9]={// z=1 =2 =3 =4 =5 =6 =7 =8 =9
1595	// {13.136,14.931,25.320,38.000,45.000,55.000,65.000,75.000,85.000},//n=1
1596	// {14.950, 2.425,11.680,18.374,27.870,35.094,48.000,60.000,72.000}, //n=2
1597	// {25.930,11.390,14.086,15.770,22.921,28.914,39.700,49.000,60.000}, //n=3
1598	// {36.830,17.594,14.908, 4.942,12.416,15.699,24.960,32.048,45.000}, //n=4
1599	// {41.860,26.110,20.946,11.348,12.051,10.650,17.338,23.111,33.610}, //n=5
1600	// {45.000,31.598,24.954,12.607, 8.668, 0.000, 5.345, 8.006,16.780}, //n=6
1601	// {50.000,40.820,33.050,20.174,13.369, 3.125, 2.863, 2.855,10.680}, //n=7
1602	// {55.000,48.810,40.796,25.076,16.562, 3.020, 0.101,-4.737,1.9520}, //n=8
1603	// {60.000,55.000,50.100,33.660,23.664, 9.873, 5.683,-0.809,0.8730}}; //n=9
1604	if(z>107)
1605	{
1606	#ifdef debug
1607	G4cout<<"***G4QPDGCode::CalcNuclMass: Z="<<z<<">107, N="<<n<<", S="<<s<<G4endl;
1608	#endif
1609	return 256000.;
1610	}
1611	G4int Z=z;
1612	G4int N=n;
1613	G4int S=s;
1614	#ifdef debug
1615	G4cout<<"G4QPDGCode::CalcNuclMass called with Z="<<Z<<",N="<<N<<", S="<<S<<G4endl;
1616	#endif
1617	if(!N&&!Z&&!S)
1618	{
1619	#ifdef debug
1620	//G4cout<<"G4QPDGCode::GetNuclMass(0,0,0)="<<mPi0<<"#0"<<G4endl;
1621	#endif
1622	//return mPi0; // Pi0 mass (dangerose)
1623	return 0.; // Photon mass
1624	}
1625	else if(!N&&!Z&&S>1) return mL*S+eps;
1626	else if(!N&&Z>1&&!S) return mP*Z+eps;
1627	else if(N>1&&!Z&&!S) return mN*N+eps;
1628	G4int A=Z+N;
1629	G4int Bn=A+S;
1630	//if((Z<0\|\|N<0)&&!Bn)
1631	//{
1632	// if (N<0) return BnmL-ZmK - NmK0+eps S ;
1633	// else return BnmL+NmPi-AmK +eps(N+S);
1634	//}
1635	//if(A<0&&Bn>=0) // BnLAMBDA's + (-(N+Z))antiKaons
1636	//{
1637	// if (N<0&&Z<0) return BnmL-ZmK -NmK0+eps S ;
1638	// else if(N<0) return BnmL+ZmPi-AmK0+eps(Z+S);
1639	// else return BnmL+NmPi-AmK +eps(N+S);
1640	//}
1641	if(!Bn) // => "GS Mesons - octet" case (without eta&pi0)
1642	{
1643	if (!S&&Z<0) return mPi*N;
1644	else if(!S&&N<0) return mPi*Z;
1645	else if ( (N == 1 && S == -1) \|\| (N == -1 && S == 1) )
1646	return mK0; // Simple decision
1647	else if ( (S == 1 && Z == -1) \|\| (S == -1 && Z == 1) )
1648	return mK; // Simple decision
1649	else if(S>0) // General decision
1650	{
1651	if (-Z>S) return SmK-(S+Z)mPi+eps;
1652	else if(Z>=0) return SmK0+ZmPi+eps;
1653	else return (S+Z)mK0-ZmK+eps;
1654	}
1655	else if(S<0) // General decision
1656	{
1657	if (Z>-S) return -SmK+(S+Z)mPi+eps;
1658	else if(Z<=0) return -SmK0-ZmPi+eps;
1659	else return -(S+Z)mK0+ZmK+eps;
1660	}
1661	}
1662	else if(Bn==1) // => "GS Baryons - octet" case (withoutSigma0)
1663	{
1664	if (Z== 1 && N== 0 && S== 0) return mP;
1665	else if(Z== 0 && N== 1 && S== 0) return mN;
1666	else if(Z== 0 && N== 0 && S== 1) return mL;
1667	else if(Z== 1 && N==-1 && S== 1) return mSp; // Lower than Sigma+ (Simp.Decis)
1668	else if(Z==-1 && N== 1 && S== 1) return mSm; // Lower than Sigma- (Simp.Decis)
1669	else if(Z== 0 && N==-1 && S== 2) return mXz; // Lower than Xi0 (Simp.Decis)
1670	else if(Z==-1 && N== 0 && S== 2) return mXm; // Lower than Xi- (Simp.Decis)
1671	else if(Z==-1 && N==-1 && S== 3) return mOm; // Lower than Omega- (Simp.Decis)
1672	else if(!S&&Z<0) return mN-mPi*Z+eps; // Negative Isonuclei
1673	else if(!S&&N<0) return mP-mPi*N+eps; // Positive Isonuclei
1674	else if(S==1) // --> General decision
1675	{
1676	if (N>1) return mSm+(N-1)mPi+eps; // (Sigma-)+(nPI-)
1677	else if(Z>1) return mSp+(Z-1)mPi+eps; // (Sigma+)+(nPI+)
1678	}
1679	else if(S==2) // --> General decision
1680	{
1681	if (N>0) return mXm+NmPi+eps; // (Xi-)+(nPI-)
1682	else if(Z>0) return mXz+ZmPi+eps; // (Xi0)+(nPI+)
1683	}
1684	else if(S==3) // --> General decision
1685	{
1686	if (N>-1) return mOm+(N+1)mPi+eps; // (Omega-)+(nPI-)
1687	else if(Z>-1) return mOm+(Z+1)mPi+eps; // (Omega-)+(nPI+)
1688	}
1689	else if(S>3) // --> General Omega- decision
1690	{
1691	if (-Z>S-2) return mOm+(S-3)mK +(2-Z-S)mPi+eps;
1692	else if(Z>-1) return mOm+(S-3)*mK0+(Z+1)+mPi+eps;
1693	else return mOm+(S+Z-2)mK0-(Z+1)mK+eps;
1694	}
1695	}
1696	else if(Bn==2) // => "GS Baryons - decuplet" case (NP,LP, and LN are made below)
1697	{
1698	if (Z== 2 && N== 0 && S== 0) return dmP;
1699	//else if(Z== 1 && N== 1 && S== 0) return 1875.61282; // Exact deuteron PDG Mass
1700	else if(Z== 1 && N== 1 && S== 0) return mD; // Exact deuteron PDG Mass
1701	else if(Z== 0 && N== 2 && S== 0) return dmN;
1702	else if(Z== 2 && N==-1 && S== 1) return dSP;
1703	else if(Z== 1 && N== 0 && S== 1) return dLP;
1704	else if(Z== 0 && N== 1 && S== 1) return dLN;
1705	else if(Z==-1 && N== 2 && S== 1) return dSN;
1706	else if(Z== 1 && N==-1 && S== 2) return dXP;
1707	else if(Z== 0 && N== 0 && S== 2) return dmL;
1708	else if(Z==-1 && N== 1 && S== 2) return dXN;
1709	else if(Z== 0 && N==-1 && S== 3) return dOP;
1710	else if(Z==-1 && N== 0 && S== 3) return dON;
1711	else if(!S&&Z<0) return dmN-mPi*Z+eps; // Negative Isonuclei
1712	else if(!S&&N<0) return dmP-mPi*N+eps; // Positive Isonuclei
1713	else if(S==1) // --> General decision
1714	{
1715	if (N>2) return dSP+(N-2)mPi+eps; // (nSigma-)+(nPI-)
1716	else if(Z>2) return dSN+(Z-1)mPi+eps; // (pSigma+)+(nPI+)
1717	}
1718	else if(S==2) // --> General decision
1719	{
1720	if (N>1) return dXN+(N-1)mPi+eps; // (nXi-)+(nPI-)
1721	else if(Z>1) return dXP+(Z-1)mPi+eps; // (pXi0)+(nPI+)
1722	}
1723	else if(S==3) // --> General decision
1724	{
1725	if (N>0) return dON+NmPi+eps; // (nOmega-)+(nPI-)
1726	else if(Z>0) return dOP+ZmPi+eps; // (pOmega-)+(nPI+)
1727	}
1728	else if(S>3) // --> General Omega- decision
1729	{
1730	if (-Z>S-2) return dON+(S-3)mK +(2-Z-S)mPi+eps;
1731	else if(Z>0) return dOP+(S-3)*mK0+Z+mPi+eps;
1732	else return dOP+(S+Z-3)mK0-ZmK+eps;
1733	}
1734	//else if(S>0) // @@ Implement General Decision
1735	//{
1736	// //#ifdef debug
1737	// G4cout<<"***G4QPDGCode::GetNuclMass:B=2, Z="<<Z<<",N="<<N<<",S="<<S<<G4endl;
1738	// //#endif
1739	// return bigM; // Exotic dibaryons (?)
1740	//}
1741	}
1742	else if(Bn==3)
1743	{
1744	if(!S) // Bn=3
1745	{
1746	if (Z==1 && N== 2) return mT; // tritium
1747	else if(Z==2 && N== 1) return mHe3; // hetrium
1748	}
1749	if(S== 1 && Z==-1 && N== 3) // nnSig-
1750	{
1751	if (Z==-1 && N== 3) return dnS; // nnSig-
1752	else if(Z== 3 && N==-1) return dpS; // ppSig+
1753	}
1754	}
1755	else if(!S)
1756	{
1757	if(Z==2 && N==2) return mAl; // Alpha
1758	else if(Z<0) return AmN-ZmPi+eps; // Multybaryon Negative Isonuclei
1759	else if(Z>A) return AmP+(Z-A)mPi+eps; // Multybaryon Positive Isonuclei
1760	}
1761	// === Start mesonic extraction ===
1762	G4double km=0.; // Mass Sum of K mesons (G4E::DecayAntiStrang.)
1763	G4int Zm=Z;
1764	G4int Nm=N;
1765	G4int Sm=S;
1766	if(S<0&&Bn>0) // NEW: the free mass minimum
1767	{
1768	if(Zm>=-S) // Enough charge for K+'s
1769	{
1770	km=-S*mK; // Anti-Lambdas are compensated by protons
1771	Zm+=S;
1772	}
1773	else if(Zm>0)
1774	{
1775	km=ZmmK-(S+Zm)mK0; // Anti-Lambdas are partially compensated by neutrons
1776	Zm=0;
1777	Nm+=S+Zm;
1778	}
1779	}
1780	else Sm=0; // No alternative calculations
1781	// Old algorithm
1782	G4double k=0.; // Mass Sum of K mesons
1783	if(S<0&&Bn>0) // OLD @@ Can be improved by K0/K+ search of minimum
1784	{
1785	G4int sH=(-S)/2; // SmallHalfS \|\| The algorithm must be the same
1786	G4int bH=-S-sH; // BigHalhS \|\| as in G4QE::DecayAntiStrange
1787	if(Z>0 && Z>N)
1788	{
1789	if(Z>=bH) // => "Enough protons in nucleus" case
1790	{
1791	if(N>=sH)
1792	{
1793	k=bHmK+sHmK0;
1794	Z-=bH;
1795	N-=sH;
1796	}
1797	else
1798	{
1799	G4int dN=Z-N;
1800	if(dN>=-S)
1801	{
1802	k=-S*mK;
1803	Z+=S;
1804	}
1805	else
1806	{
1807	G4int sS=(-S-dN)/2;
1808	G4int bS=-S-dN-sS;
1809	sS+=dN;
1810	if(Z>=sS&&N>=bS)
1811	{
1812	k=sSmK+bSmK0;
1813	Z-=sS;
1814	N-=bS;
1815	}
1816	else if(Z<sS)
1817	{
1818	G4int dS=-S-Z;
1819	k=ZmK+dSmK0;
1820	N-=dS;
1821	Z=0;
1822	}
1823	else
1824	{
1825	G4int dS=-S-N;
1826	k=dSmK+NmK0;
1827	Z-=dS;
1828	N=0;
1829	}
1830	}
1831	}
1832	}
1833	else // Must not be here
1834	{
1835	#ifdef debug
1836	G4cout<<"***G4QPDGC::CalcNuclMass:Antimatter? Z="<<Z<<",N="<<N<<",S="<<S<<G4endl;
1837	#endif
1838	return 0.; // @@ Antiparticles aren't implemented @@
1839	}
1840	}
1841	else if(N>=bH)
1842	{
1843	if(Z>=sH)
1844	{
1845	k=sHmK+bHmK0;
1846	Z-=sH;
1847	N-=bH;
1848	}
1849	else
1850	{
1851	G4int dN=N-Z;
1852	if(dN>=-S)
1853	{
1854	k=-S*mK0;
1855	N+=S;
1856	}
1857	else
1858	{
1859	G4int sS=(-S-dN)/2;
1860	G4int bS=-S-dN-sS;
1861	bS+=dN;
1862	if(N>=bS&&Z>=sS)
1863	{
1864	k=sSmK+bSmK0;
1865	Z-=sS;
1866	N-=bS;
1867	}
1868	else if(N<bS)
1869	{
1870	G4int dS=-S-N;
1871	k=dSmK+NmK0;
1872	Z-=dS;
1873	N=0;
1874	}
1875	else
1876	{
1877	G4int dS=-S-Z;
1878	k=ZmK+dSmK0;
1879	N-=dS;
1880	Z=0;
1881	}
1882	}
1883	}
1884	}
1885	else // Must not be here
1886	{
1887	return 0.; // @@ Antiparticles aren't implemented @@
1888	#ifdef debug
1889	G4cout<<"***G4QPDGC::CalcNuclMass:Antimatter? N="<<N<<",Z="<<Z<<",S="<<S<<G4endl;
1890	#endif
1891	}
1892	S=0;
1893	}
1894	if(N<0)
1895	{
1896	k+=-N*mPi;
1897	Z+=N;
1898	N=0;
1899	}
1900	if(Z<0)
1901	{
1902	k+=-Z*mPi;
1903	N+=Z;
1904	Z=0;
1905	}
1906	A=Z+N;
1907	if (!A) return k+SmL+Seps; // @@ multy LAMBDA states are not implemented
1908	G4double m=k+A*um; // Expected mass in atomic units
1909	//G4double D=N-Z; // Isotopic shift of the nucleus
1910	if ( (A+S < 1 && k==0.) \|\| Z < 0 \|\| N < 0 )
1911	{ // @@ Can be generalized to anti-nuclei
1912	#ifdef debug
1913	G4cout<<"**G4QPDGCode::CalcNuclMass:A="<<A<<"<1 \|\| Z="<<Z<<"<0 \|\| N="<<N<<"<0"<<G4endl;
1914	//@@throw G4QException("***G4QPDGCode::GetNuclMass: Impossible nucleus");
1915	#endif
1916	return 0.; // @@ Temporary
1917	}
1918	if (!Z) return k+N(mN+.1)+S(mL+.1); // @@ n+LAMBDA states are not implemented
1919	else if(!N) return k+Z(mP+1.)+S(mL+.1); // @@ p+LAMBDA states are not implemented
1920	//else if(N<=9&&Z<=9) m+=1.433e-5pow(double(Z),2.39)-Zme+c[N-1][Z-1];// Geant4 Comp.now
1921	else
1922	{
1923	//G4double fA=A;
1924	// IsInTable is inside G4NucleiProperties::GetNuclearMass
1925	// G4ParticleTable::GetParticleTable()->FindIon(Zm,Am,0,Zm) creates new Ion!
1926	if(A==256 && Z==128) m=256000.;
1927	else m=k+G4NucleiProperties::GetNuclearMass(A,Z);
1928	//if(G4NucleiPropertiesTable::IsInTable(Z,A))
1929	// m=k+G4NucleiProperties::GetNuclearMass(A,Z);
1930	//else if(A==256 && Z==128) m=256000.;
1931	//else
1932	// m=k+G4ParticleTable::GetParticleTable()->FindIon(Z,A,0,Z)->GetPDGMass();
1933	//m+=-sh[Z]-sh[N]+b1DDpow(fA,b2)+b3(1.-2./(1.+exp(b4D)))+ZZ(b5pow(fA,b9)+b6/fA);
1934	}
1935	//@@//G4double maxM= k+ZmP+NmN+S*mL+eps; // @@ eps -- Wings of the Mass parabola
1936	//@@//if(m>maxM) m=maxM;
1937	G4double mm=m;
1938	if(Sm<0) // For the new algorithm of calculation
1939	{
1940	if(Nm<0)
1941	{
1942	km+=-Nm*mPi;
1943	Zm+=Nm;
1944	Nm=0;
1945	}
1946	if(Zm<0)
1947	{
1948	km+=-Zm*mPi;
1949	Nm+=Zm;
1950	Zm=0;
1951	}
1952	G4int Am=Zm+Nm;
1953	if(!Am) return km+eps;
1954	mm=km+Am*um; // Expected mass in atomic units
1955	//G4double Dm=Nm-Zm; // Isotopic shift of the nucleus
1956	if ( (Am < 1 && km==0.) \|\| Zm < 0 \|\| Nm < 0 )
1957	{ // @@ Can be generalized to anti-nuclei
1958	#ifdef debug
1959	G4cerr<<"*G4QPDGCode::CalcNucM:A="<<Am<<"<1 \|\| Z="<<Zm<<"<0 \|\| N="<<Nm<<"<0"<<G4endl;
1960	#endif
1961	}
1962	if (!Zm) return km+Nm*(mN+.1);
1963	else if(!Nm) return km+Zm*(mP+1.);
1964	//else if(Nm<=9&&Zm<=9) mm+=1.433e-5pow(double(Zm),2.39)-Zmme+c[Nm-1][Zm-1];// Geant4
1965	else
1966	{
1967	// IsInTable is inside G4NucleiProperties::GetNuclearMass
1968	// G4ParticleTable::GetParticleTable()->FindIon(Zm,Am,0,Zm) creates new Ion!
1969	mm=km+G4NucleiProperties::GetNuclearMass(Am,Zm);
1970	//G4double fA=Am;
1971	//if(G4NucleiPropertiesTable::IsInTable(Zm,Am))
1972	// mm=km+G4NucleiProperties::GetNuclearMass(Am,Zm);
1973	//else
1974	// mm=km+G4ParticleTable::GetParticleTable()->FindIon(Zm,Am,0,Zm)->GetPDGMass();
1975	// //mm+=-sh[Zm]-sh[Nm]+b1DmDmpow(fA,b2)+b3(1.-2./(1.+exp(b4*Dm)))
1976	// // +ZmZm(b5*pow(fA,b9)+b6/Am);
1977	}
1978	//@@//G4double mM= km+ZmmP+NmmN+eps;
1979	//@@//if(mm>mM) mm=mM;
1980	}
1981	if(m>mm) m=mm;
1982	if(S>0)
1983	{
1984	G4double bs=0.;
1985	if (A==2) bs=a2;
1986	else if(A==3) bs=a3;
1987	else if(A>3) bs=b7*exp(-b8/(A+1.));
1988	m+=S*(mL-bs);
1989	}
1990	#ifdef debug
1991	G4cout<<"G4QPDGCode::CalcNuclMass: >>>OUT<<< m="<<m<<G4endl;
1992	#endif
1993	//if(z==8&&n==9&&!s) G4cout<<"G4QPDGC::GetNucM:m="<<m<<",mm="<<mm<<G4endl;
1994	return m;
1995	}
1996
1997	// Get a Quark Content {d,u,s,ad,au,as} for the PDG
1998	G4QContent G4QPDGCode::GetQuarkContent() const
1999	// ======================================
2000	{
2001	G4int a=0; // particle
2002	if(thePDGCode<0) a=1; // anti-particle
2003	G4int d=0;
2004	G4int u=0;
2005	G4int s=0;
2006	G4int ad=0;
2007	G4int au=0;
2008	G4int as=0;
2009	G4int ab=abs(thePDGCode);
2010	if (!ab)
2011	{
2012	G4cerr<<"***G4QPDGCode::GetQuarkContent: PDG=0, return (0,0,0,0,0,0)"<<G4endl;
2013	return G4QContent(0,0,0,0,0,0);
2014	}
2015	else if(ab<4) // @@ for SU(6) : ab<7
2016	{
2017	if (thePDGCode== 1) return G4QContent(1,0,0,0,0,0);
2018	else if(thePDGCode== 2) return G4QContent(0,1,0,0,0,0);
2019	else if(thePDGCode== 3) return G4QContent(0,0,1,0,0,0);
2020	else if(thePDGCode==-1) return G4QContent(0,0,0,1,0,0);
2021	else if(thePDGCode==-2) return G4QContent(0,0,0,0,1,0);
2022	else if(thePDGCode==-3) return G4QContent(0,0,0,0,0,1);
2023	}
2024	else if(ab<99)
2025	{
2026	#ifdef debug
2027	if (ab==22) G4cout<<"-W-G4QPDGC::GetQuarkCont: For the Photon? - Return 0"<<G4endl;
2028	else if(ab==10) G4cout<<"-W-G4QPDGC::GetQuarkCont: For Chipolino? - Return 0"<<G4endl;
2029	else G4cout<<"-W-G4QPDGCode::GetQuarkCont: For PDG="<<thePDGCode<<" Return 0"<<G4endl;
2030	#endif
2031	return G4QContent(0,0,0,0,0,0); // Photon, bosons, leptons
2032	}
2033	else if(ab<80000000) // Baryons & Mesons
2034	{
2035	G4int c=ab/10; // temporary (quarks)
2036	G4int f=c%10; // (1) low quark
2037	G4int v=c/10; // (2,3) temporary(B) or (2) final(M) (high quarks, high quark)
2038	G4int t=0; // (3)prototype of highest quark (B)
2039	#ifdef sdebug
2040	G4cout<<"G4QPDGCode::GetQuarkContent: a="<<ab<<", c="<<c<<", f="<<f<<", v="<<v<<G4endl;
2041	#endif
2042	if(v>10) // Baryons
2043	{
2044	t=v/10; // (3) highest quark
2045	v%=10; // (2) high quark
2046	if (f==1)
2047	{
2048	if(a) ad++;
2049	else d++;
2050	}
2051	else if(f==2)
2052	{
2053	if(a) au++;
2054	else u++;
2055	}
2056	else if(f==3)
2057	{
2058	if(a) as++;
2059	else s++;
2060	}
2061	else G4cerr<<"*G4QPDGC::GetQCont:1 PDG="<<thePDGCode<<","<<f<<","<<v<<","<<t<<G4endl;
2062	if (v==1)
2063	{
2064	if(a) ad++;
2065	else d++;
2066	}
2067	else if(v==2)
2068	{
2069	if(a) au++;
2070	else u++;
2071	}
2072	else if(v==3)
2073	{
2074	if(a) as++;
2075	else s++;
2076	}
2077	else G4cerr<<"*G4QPDGC::GetQCont:2 PDG="<<thePDGCode<<","<<f<<","<<v<<","<<t<<G4endl;
2078	if (t==1)
2079	{
2080	if(a) ad++;
2081	else d++;
2082	}
2083	else if(t==2)
2084	{
2085	if(a) au++;
2086	else u++;
2087	}
2088	else if(t==3)
2089	{
2090	if(a) as++;
2091	else s++;
2092	}
2093	else G4cerr<<"*G4QPDGC::GetQCont:3 PDG="<<thePDGCode<<","<<f<<","<<v<<","<<t<<G4endl;
2094	return G4QContent(d,u,s,ad,au,as);
2095	}
2096	else // Mesons
2097	{
2098	if(f==v)
2099	{
2100	if (f==1) return G4QContent(1,0,0,1,0,0);
2101	else if(f==2) return G4QContent(0,1,0,0,1,0);
2102	else if(f==3) return G4QContent(0,0,1,0,0,1);
2103	else G4cerr<<"*G4QPDGC::GetQC:4 PDG="<<thePDGCode<<","<<f<<","<<v<<","<<t<<G4endl;
2104	}
2105	else
2106	{
2107	if (f==1 && v==2)
2108	{
2109	if(a)return G4QContent(1,0,0,0,1,0);
2110	else return G4QContent(0,1,0,1,0,0);
2111	}
2112	else if(f==1 && v==3)
2113	{
2114	if(a)return G4QContent(0,0,1,1,0,0);
2115	else return G4QContent(1,0,0,0,0,1);
2116	}
2117	else if(f==2 && v==3)
2118	{
2119	if(a)return G4QContent(0,0,1,0,1,0);
2120	else return G4QContent(0,1,0,0,0,1);
2121	}
2122	else G4cerr<<"*G4QPDGC::GetQC:5 PDG="<<thePDGCode<<","<<f<<","<<v<<","<<t<<G4endl;
2123	}
2124	}
2125	}
2126	else
2127	{
2128	G4int szn=ab-90000000;
2129	G4int ds=0;
2130	G4int dz=0;
2131	G4int dn=0;
2132	if(szn<-100000)
2133	{
2134	G4int ns=(-szn)/1000000+1;
2135	szn+=ns*1000000;
2136	ds+=ns;
2137	}
2138	else if(szn<-100)
2139	{
2140	G4int nz=(-szn)/1000+1;
2141	szn+=nz*1000;
2142	dz+=nz;
2143	}
2144	else if(szn<0)
2145	{
2146	G4int nn=-szn;
2147	szn=0;
2148	dn+=nn;
2149	}
2150	G4int sz =szn/1000;
2151	G4int n =szn%1000;
2152	if(n>700)
2153	{
2154	n-=1000;
2155	dz--;
2156	}
2157	G4int z =sz%1000-dz;
2158	if(z>700)
2159	{
2160	z-=1000;
2161	ds--;
2162	}
2163	s =sz/1000-ds;
2164	G4int b=z+n+s;
2165	d=n+b;
2166	u=z+b;
2167	if (d<0&&u<0&&s<0) return G4QContent(0,0,0,-d,-u,-s);
2168	else if (u<0&&s<0) return G4QContent(d,0,0,0,-u,-s);
2169	else if (d<0&&s<0) return G4QContent(0,u,0,-d,0,-s);
2170	else if (d<0&&u<0) return G4QContent(0,0,s,-d,-u,0);
2171	else if (u<0) return G4QContent(d,0,s,0,-u,0);
2172	else if (s<0) return G4QContent(d,u,0,0,0,-s);
2173	else if (d<0) return G4QContent(0,u,s,-d,0,0);
2174	else return G4QContent(d,u,s,0,0,0);
2175	}
2176	return G4QContent(0,0,0,0,0,0);
2177	}
2178
2179	// Quark Content of pseudo-meson for q_i->q_o Quark Exchange
2180	G4QContent G4QPDGCode::GetExQContent(G4int i, G4int o) const
2181	{// ==================================================
2182	G4QContent cQC(0,0,0,0,0,0);
2183	if (!i) cQC.IncAD();
2184	else if(i==1) cQC.IncAU();
2185	else if(i==2) cQC.IncAS();
2186	else G4cerr<<"***G4QPDGCode::GetExQContent: strange entering quark i="<<i<<G4endl;
2187	if (!o) cQC.IncD();
2188	else if(o==1) cQC.IncU();
2189	else if(o==2) cQC.IncS();
2190	else G4cerr<<"***G4QPDGCode::GetExQContent: strange exiting quark o="<<o<<G4endl;
2191	return cQC;
2192	}
2193
2194
2195	// Relative Cross Index for q_i->q_o (q=0(d),1(u),2(s), 7 means there is no such cluster)
2196	G4int G4QPDGCode::GetRelCrossIndex(G4int i, G4int o) const
2197	{// =====================================================
2198	// pD,nD,DD,DD,ppDnnDpDDnDD n, p, L,S-,S+,X-,X0,O-
2199	static const G4int b01[16]={ 7,15, 7, 7, 7, 7, 7, 7, 1, 7, 7,-1, 7, 7, 7, 7};
2200	static const G4int b02[16]={ 7, 7, 7, 7, 7, 7, 7, 7, 2, 7, 7, 7, 7, 7, 7, 7};
2201	static const G4int b10[16]={18, 7, 7, 7, 7, 7, 7, 7, 7,-1, 7, 7,-2, 7, 7, 7}; // Iso+Bary
2202	static const G4int b12[16]={ 7, 7, 7, 7, 7, 7, 7, 7, 7, 1, 7, 7, 7, 7, 7, 7};
2203	static const G4int b20[16]={ 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,-2, 7,-3, 7,-4, 7};
2204	static const G4int b21[16]={ 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,-1,-3, 7,-3, 7, 7};
2205	// nn,np,pp,Ln,Lp,LL,nS,pS,nnpnppLnnLpnLppLLnLLpnnS
2206	static const G4int d01[16]={ 1, 1, 7, 1, 7, 7,-3, 7, 1, 7, 1, 1, 7, 1, 7,-5};
2207	static const G4int d02[16]={ 3, 3, 7, 2, 7, 7, 7, 7, 3, 3, 3, 3, 7, 7, 7, 7};
2208	static const G4int d10[16]={ 7,-1,-1, 7,-1, 7, 7,-4, 7,-1, 7,-1,-1, 7,-1, 7}; //B=2,3
2209	static const G4int d12[16]={ 7, 2, 2, 7, 1, 7, 7, 7, 2, 2, 7, 2, 2, 7, 7, 7};
2210	static const G4int d20[16]={ 7, 7, 7,-3,-3,-2, 7,-5, 7, 7, 7,-3,-3,-3,-3, 7};
2211	static const G4int d21[16]={ 7, 7, 7,-2,-2,-1,-6, 7, 7, 7,-2,-2, 7,-2,-2, 7};
2212	// nn,np,pp,Ln,Lp,nn,np,pp! n, p,nn,np,pp,Ln,Lp
2213	static const G4int m01[15]={ 1, 1, 7, 1, 7, 1, 1, 7, 1, 7, 1, 1, 7, 1, 7};// Multibaryons
2214	static const G4int m02[15]={ 3, 3, 7, 3, 3, 7, 7, 7, 3, 3, 3, 3, 7, 7, 7};// @@ Regular !
2215	static const G4int m10[15]={ 7,-1,-1, 7,-1, 7,-1,-1, 7,-1, 7,-1,-1, 7,-1};// 01->1,10->-1
2216	static const G4int m12[15]={ 7, 2, 2, 2, 2, 7, 7, 7, 2, 2, 7, 2, 2, 7, 7};// 12->2,21->-2
2217	static const G4int m20[15]={ 7, 7, 7,-3,-3, 7,-3,-3, 7, 7, 7,-3,-3,-3,-3};// 02->3,20->-3
2218	static const G4int m21[15]={ 7, 7, 7,-2,-2,-2,-2, 7, 7, 7,-2,-2, 7,-2,-2};
2219	static const G4int fragmStart = 82; // "Isonuclei are added && Leptons are added"
2220
2221	if(theQCode<fragmStart) return 7;
2222	G4int sub=theQCode-fragmStart;
2223	if ( (sub > 1 && sub < 8) \|\| sub == 15) return 7;//Why they are in clusters?-Residuals(?)
2224	G4int rel=sub; // case of nuclear baryons and isonuclei
2225	if (sub>31)rel =(sub-32)%15; // case of heavy fragments (BaryNum>3)
2226	else if(sub>15)rel = sub-16; // case of nuclear di-baryon & tri-baryons
2227	#ifdef debug
2228	G4cout<<"G4QPDGCode::RelGetCrossIndex:i="<<i<<",o="<<o<<",su="<<sub<<",re="<<rel<<G4endl;
2229	#endif
2230	if (!i) // ==> input quark = 0(d) (d=-1/3)
2231	{
2232	if (!o) return 0; // -> output quark = 0(d) => 0 = the same cluster
2233	else if(o==1) // -> output quark = 1(u) (ch--)
2234	{
2235	if (sub<16) return b01[rel];
2236	else if(sub<32) return d01[rel];
2237	else return m01[rel];
2238	}
2239	else if(o==2)
2240	{
2241	if (sub<16) return b02[rel];
2242	else if(sub<32) return d02[rel];
2243	else return m02[rel];
2244	}
2245	}
2246	else if(i==1) // ==> input quark = 1(u) (u=2/3)
2247	{
2248	if (!o) // -> output quark = 0(d) (ch++)
2249	{
2250	if (sub<16) return b10[rel];
2251	else if(sub<32) return d10[rel];
2252	else return m10[rel];
2253	}
2254	else if(o==1) return 0; // -> output quark = 1(u) => 0 = the same cluster
2255	else if(o==2)
2256	{
2257	if (sub<16) return b12[rel];
2258	else if(sub<32) return d12[rel];
2259	else return m12[rel];
2260	}
2261	}
2262	else if(i==2)
2263	{
2264	if (!o)
2265	{
2266	if (sub<16) return b20[rel];
2267	else if(sub<32) return d20[rel];
2268	else return m20[rel];
2269	}
2270	else if(o==1)
2271	{
2272	if (sub<16) return b21[rel];
2273	else if(sub<32) return d21[rel];
2274	else return m21[rel];
2275	}
2276	else if(o==2) return 0;
2277	}
2278	return 7;
2279	}
2280
2281	// Get number of Combinations for q_i->q_o
2282	G4int G4QPDGCode::GetNumOfComb(G4int i, G4int o) const
2283	{// ================================================
2284	if(i>-1&&i<3)
2285	{
2286	G4int shiftQ=GetRelCrossIndex(i, o);
2287	G4int sQCode=theQCode; // QCode of the parent cluster
2288	if (shiftQ==7) return 0; // A parent cluster doesn't exist
2289	else if(!shiftQ) sQCode+=shiftQ; // Shift QCode of son to QCode of his parent
2290	G4QPDGCode parent; // Create a temporary (fake) parent cluster
2291	parent.InitByQCode(sQCode); // Initialize it by Q-Code
2292	G4QContent parentQC=parent.GetQuarkContent();// Quark Content of the parent cluster
2293	if (!o) return parentQC.GetD();
2294	else if(o==1) return parentQC.GetU();
2295	else if(o==2) return parentQC.GetS();
2296	else G4cerr<<"***G4QPDGCode:::GetNumOfComb: strange exiting quark o="<<o<<G4endl;
2297	}
2298	else G4cerr<<"***G4QPDGCode:::GetNumOfComb: strange entering quark i="<<i<<G4endl;
2299	return 0;
2300	}
2301
2302	// Get a total number of Combinations for q_i
2303	G4int G4QPDGCode::GetTotNumOfComb(G4int i) const
2304	{// ==========================================
2305	G4int tot=0;
2306	if(i>-1&&i<3) for(int j=0; j<3; j++) tot+=GetNumOfComb(i, j);
2307	else G4cerr<<"***G4QPDGCode:::GetTotNumOfComb: strange entering quark i="<<i<<G4endl;
2308	return tot;
2309	}
2310
2311	// Converts nuclear PDGCode to Z(#of protons), N(#of neutrons), S(#of lambdas) values
2312	void G4QPDGCode::ConvertPDGToZNS(G4int nucPDG, G4int& z, G4int& n, G4int& s)
2313	{// =======================================================================
2314	if(nucPDG>80000000&&nucPDG<100000000) // Condition of conversion
2315	{
2316	z=0;
2317	n=0;
2318	s=0;
2319	G4int r=nucPDG;
2320	if(r==90000000) return;
2321	G4int cn =r%1000; // candidate to #of neutrons
2322	if(cn)
2323	{
2324	if(cn>500) cn-=1000; // AntiNeutrons
2325	n=cn; // Increment neutrons
2326	r-=cn; // Subtract them from the residual
2327	if(r==90000000) return;
2328	}
2329	G4int cz =r%1000000; // candidate to #of neutrons
2330	if(cz)
2331	{
2332	if(cz>500000) cz-=1000000; // AntiProtons
2333	z=cz/1000; // Number of protons
2334	r-=cz; // Subtract them from the residual
2335	if(r==90000000) return;
2336	}
2337	G4int cs =r%10000000; // candidate to #of neutrons
2338	if(cs)
2339	{
2340	if(cs>5000000) cs-=10000000; // AntiLambda
2341	s=cs/1000000; // Number of Lambdas
2342	}
2343	}
2344	return;
2345	}

Note: See TracBrowser for help on using the repository browser.

Download in other formats: