Context Navigation

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

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

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

Download in other formats: