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