Context Navigation

G4QPDGCode.cc @ 1350

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

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

Context Navigation

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

Download in other formats: