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