| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // INCL++ intra-nuclear cascade model 26 // INCL++ intra-nuclear cascade model
27 // Alain Boudard, CEA-Saclay, France << 27 // Pekka Kaitaniemi, CEA and Helsinki Institute of Physics
28 // Joseph Cugnon, University of Liege, Belgium << 28 // Davide Mancusi, CEA
29 // Jean-Christophe David, CEA-Saclay, France << 29 // Alain Boudard, CEA
30 // Pekka Kaitaniemi, CEA-Saclay, France, and H << 30 // Sylvie Leray, CEA
31 // Sylvie Leray, CEA-Saclay, France << 31 // Joseph Cugnon, University of Liege
32 // Davide Mancusi, CEA-Saclay, France <<
33 // 32 //
34 #define INCLXX_IN_GEANT4_MODE 1 33 #define INCLXX_IN_GEANT4_MODE 1
35 34
36 #include "globals.hh" 35 #include "globals.hh"
37 36
38 #include "G4INCLParticleTable.hh" 37 #include "G4INCLParticleTable.hh"
39 #include "G4INCLNuclearMassTable.hh" <<
40 #include <algorithm> 38 #include <algorithm>
41 // #include <cassert> 39 // #include <cassert>
42 #include <cmath> 40 #include <cmath>
43 #include <cctype> 41 #include <cctype>
44 #include <sstream> 42 #include <sstream>
45 #ifdef INCLXX_IN_GEANT4_MODE <<
46 #include "G4SystemOfUnits.hh" <<
47 #endif <<
48 43
49 #ifdef INCLXX_IN_GEANT4_MODE 44 #ifdef INCLXX_IN_GEANT4_MODE
50 #include "G4PhysicalConstants.hh" 45 #include "G4PhysicalConstants.hh"
51 #include "G4SystemOfUnits.hh" 46 #include "G4SystemOfUnits.hh"
52 #endif 47 #endif
53 48
54 namespace G4INCL { 49 namespace G4INCL {
55 50
56 namespace ParticleTable { << 51 /// \brief Static instance of the NaturalIsotopicAbundances class
57 << 52 const NaturalIsotopicDistributions *ParticleTable::theNaturalIsotopicDistributions = NULL;
58 namespace { <<
59 <<
60 /// \brief Static instance of the Natura <<
61 const NaturalIsotopicDistributions *theN <<
62 <<
63 const G4double theINCLNucleonMass = 938. <<
64 const G4double theINCLPionMass = 138.0; <<
65 const G4double theINCLLambdaMass = 1115. <<
66 //const G4double theINCLSigmaMass = 1197 <<
67 //const G4double theINCLKaonMass = 497.6 <<
68 const G4double theINCLEtaMass = 547.862; <<
69 const G4double theINCLOmegaMass = 782.65 <<
70 const G4double theINCLEtaPrimeMass = 957 <<
71 const G4double theINCLPhotonMass = 0.0; <<
72 G4ThreadLocal G4double protonMass = 0.0; <<
73 G4ThreadLocal G4double neutronMass = 0.0 <<
74 G4ThreadLocal G4double piPlusMass = 0.0; <<
75 G4ThreadLocal G4double piMinusMass = 0.0 <<
76 G4ThreadLocal G4double piZeroMass = 0.0; <<
77 G4ThreadLocal G4double SigmaPlusMass = 0 <<
78 G4ThreadLocal G4double SigmaZeroMass = 0 <<
79 G4ThreadLocal G4double SigmaMinusMass = <<
80 G4ThreadLocal G4double LambdaMass = 0.0; <<
81 G4ThreadLocal G4double XiMinusMass = 0.0 <<
82 G4ThreadLocal G4double XiZeroMass = 0.0; <<
83 G4ThreadLocal G4double antiProtonMass = <<
84 G4ThreadLocal G4double antiNeutronMass = <<
85 G4ThreadLocal G4double antiSigmaPlusMass <<
86 G4ThreadLocal G4double antiSigmaZeroMass <<
87 G4ThreadLocal G4double antiSigmaMinusMas <<
88 G4ThreadLocal G4double antiLambdaMass = <<
89 G4ThreadLocal G4double antiXiMinusMass = <<
90 G4ThreadLocal G4double antiXiZeroMass = <<
91 G4ThreadLocal G4double KPlusMass = 0.0; <<
92 G4ThreadLocal G4double KZeroMass = 0.0; <<
93 G4ThreadLocal G4double KZeroBarMass = 0. <<
94 G4ThreadLocal G4double KShortMass = 0.0; <<
95 G4ThreadLocal G4double KLongMass = 0.0; <<
96 G4ThreadLocal G4double KMinusMass = 0.0; <<
97 G4ThreadLocal G4double etaMass = 0.0; <<
98 G4ThreadLocal G4double omegaMass = 0.0; <<
99 G4ThreadLocal G4double etaPrimeMass = 0. <<
100 G4ThreadLocal G4double photonMass = 0.0; <<
101 <<
102 // Hard-coded values of the real particl <<
103 G4ThreadLocal G4double theRealProtonMass <<
104 G4ThreadLocal G4double theRealNeutronMas <<
105 G4ThreadLocal G4double theRealChargedPiM <<
106 G4ThreadLocal G4double theRealPiZeroMass <<
107 G4ThreadLocal G4double theRealLambdaMass <<
108 G4ThreadLocal G4double theRealSigmaPlusM <<
109 G4ThreadLocal G4double theRealSigmaZeroM <<
110 G4ThreadLocal G4double theRealSigmaMinus <<
111 G4ThreadLocal G4double theRealAntiProton <<
112 G4ThreadLocal G4double theRealXiMinusMas <<
113 G4ThreadLocal G4double theRealXiZeroMass <<
114 G4ThreadLocal G4double theRealAntiNeutro <<
115 G4ThreadLocal G4double theRealAntiLambda <<
116 G4ThreadLocal G4double theRealAntiSigmaP <<
117 G4ThreadLocal G4double theRealAntiSigmaZ <<
118 G4ThreadLocal G4double theRealAntiSigmaM <<
119 G4ThreadLocal G4double theRealAntiXiMinu <<
120 G4ThreadLocal G4double theRealAntiXiZero <<
121 G4ThreadLocal G4double theRealChargedKao <<
122 G4ThreadLocal G4double theRealNeutralKao <<
123 G4ThreadLocal G4double theRealEtaMass = <<
124 G4ThreadLocal G4double theRealOmegaMass <<
125 G4ThreadLocal G4double theRealEtaPrimeMa <<
126 G4ThreadLocal G4double theRealPhotonMass <<
127 <<
128 // Width (second) <<
129 const G4double theChargedPiWidth = 2.603 <<
130 const G4double thePiZeroWidth = 8.52e-17 <<
131 const G4double theEtaWidth = 5.025e-19; <<
132 const G4double theOmegaWidth = 7.7528e-2 <<
133 const G4double theEtaPrimeWidth = 3.3243 <<
134 const G4double theChargedKaonWidth = 1.2 <<
135 const G4double theKShortWidth = 8.954e-1 <<
136 const G4double theKLongWidth = 5.116e-08 <<
137 const G4double theLambdaWidth = 2.632e-1 <<
138 const G4double theSigmaPlusWidth = 8.018 <<
139 const G4double theSigmaZeroWidth = 7.4e- <<
140 const G4double theSigmaMinusWidth = 1.47 <<
141 //const G4double theXiMinusWidth = 1.639 <<
142 //const G4double theXiZeroWidth = 2.90e- <<
143 //const G4double theAntiLambdaWidth = 2. <<
144 //const G4double theAntiSigmaPlusWidth = <<
145 //const G4double theAntiSigmaZeroWidth = <<
146 //const G4double theAntiSigmaMinusWidth <<
147 //const G4double theAntiXiMinusWidth = 1 <<
148 //const G4double theAntiXiZeroWidth = 2. <<
149 G4ThreadLocal G4double piPlusWidth = 0.0 <<
150 G4ThreadLocal G4double piMinusWidth = 0. <<
151 G4ThreadLocal G4double piZeroWidth = 0.0 <<
152 G4ThreadLocal G4double etaWidth = 0.0; <<
153 G4ThreadLocal G4double omegaWidth = 0.0; <<
154 G4ThreadLocal G4double etaPrimeWidth = 0 <<
155 G4ThreadLocal G4double LambdaWidth = 0.0 <<
156 G4ThreadLocal G4double SigmaPlusWidth = <<
157 G4ThreadLocal G4double SigmaZeroWidth = <<
158 G4ThreadLocal G4double SigmaMinusWidth = <<
159 G4ThreadLocal G4double KPlusWidth = 0.0; <<
160 G4ThreadLocal G4double KMinusWidth = 0.0 <<
161 G4ThreadLocal G4double KShortWidth = 0.0 <<
162 G4ThreadLocal G4double KLongWidth = 0.0; <<
163 G4ThreadLocal G4double XiMinusWidth = 0. <<
164 G4ThreadLocal G4double XiZeroWidth = 0.0 <<
165 G4ThreadLocal G4double antiLambdaWidth = <<
166 G4ThreadLocal G4double antiSigmaZeroWidt <<
167 G4ThreadLocal G4double antiSigmaMinusWid <<
168 G4ThreadLocal G4double antiSigmaPlusWidt <<
169 G4ThreadLocal G4double antiXiZeroWidth = <<
170 G4ThreadLocal G4double antiXiMinusWidth <<
171 <<
172 const G4int mediumNucleiTableSize = 30; <<
173 <<
174 const G4double mediumDiffuseness[mediumN <<
175 {0.0,0.0,0.0,0.0,0.0,1.78,1.77,1.77,1.69 <<
176 1.69,1.72,1.635,1.730,1.81,1.833,1.798 <<
177 1.93,0.567,0.571, 0.560,0.549,0.550,0. <<
178 0.580,0.575,0.569,0.537,0.0,0.0}; <<
179 const G4double mediumRadius[mediumNuclei <<
180 {0.0,0.0,0.0,0.0,0.0,0.334,0.327,0.479,0 <<
181 0.811,0.84,1.403,1.335,1.25,1.544,1.49 <<
182 2.58,2.77, 2.775,2.78,2.88,2.98,3.22,3 <<
183 3.14,0.0,0.0}; <<
184 <<
185 const G4double positionRMS[clusterTableZ <<
186 /* A= 0 1 2 3 4 <<
187 /* Z=0 */ {-1.0, -1.0, -1.0, -1.0, -1. <<
188 /* Z=1 */ {-1.0, -1.0, 2.10, 1.80, 1.7 <<
189 /* Z=2 */ {-1.0, -1.0, -1.0, 1.80, 1.6 <<
190 /* Z=3 */ {-1.0, -1.0, -1.0, -1.0, 1.7 <<
191 /* Z=4 */ {-1.0, -1.0, -1.0, -1.0, -1. <<
192 /* Z=5 */ {-1.0, -1.0, -1.0, -1.0, -1. <<
193 /* Z=6 */ {-1.0, -1.0, -1.0, -1.0, -1. <<
194 /* Z=7 */ {-1.0, -1.0, -1.0, -1.0, -1. <<
195 /* Z=8 */ {-1.0, -1.0, -1.0, -1.0, -1. <<
196 }; <<
197 <<
198 const G4double momentumRMS[clusterTableZ <<
199 /* A= 0 1 2 3 4 <<
200 /* Z=0 */ {-1.0, -1.0, -1.0, -1.0, -1. <<
201 /* Z=1 */ {-1.0, -1.0, 77.0, 110., 153 <<
202 /* Z=2 */ {-1.0, -1.0, -1.0, 110., 153 <<
203 /* Z=3 */ {-1.0, -1.0, -1.0, -1.0, 153 <<
204 /* Z=4 */ {-1.0, -1.0, -1.0, -1.0, -1. <<
205 /* Z=5 */ {-1.0, -1.0, -1.0, -1.0, -1. <<
206 /* Z=6 */ {-1.0, -1.0, -1.0, -1.0, -1. <<
207 /* Z=7 */ {-1.0, -1.0, -1.0, -1.0, -1. <<
208 /* Z=8 */ {-1.0, -1.0, -1.0, -1.0, -1. <<
209 }; <<
210 <<
211 const G4int elementTableSize = 113; // u <<
212 <<
213 /// \brief Table of chemical element nam <<
214 const std::string elementTable[elementTa <<
215 "", <<
216 "H", <<
217 "He", <<
218 "Li", <<
219 "Be", <<
220 "B", <<
221 "C", <<
222 "N", <<
223 "O", <<
224 "F", <<
225 "Ne", <<
226 "Na", <<
227 "Mg", <<
228 "Al", <<
229 "Si", <<
230 "P", <<
231 "S", <<
232 "Cl", <<
233 "Ar", <<
234 "K", <<
235 "Ca", <<
236 "Sc", <<
237 "Ti", <<
238 "V", <<
239 "Cr", <<
240 "Mn", <<
241 "Fe", <<
242 "Co", <<
243 "Ni", <<
244 "Cu", <<
245 "Zn", <<
246 "Ga", <<
247 "Ge", <<
248 "As", <<
249 "Se", <<
250 "Br", <<
251 "Kr", <<
252 "Rb", <<
253 "Sr", <<
254 "Y", <<
255 "Zr", <<
256 "Nb", <<
257 "Mo", <<
258 "Tc", <<
259 "Ru", <<
260 "Rh", <<
261 "Pd", <<
262 "Ag", <<
263 "Cd", <<
264 "In", <<
265 "Sn", <<
266 "Sb", <<
267 "Te", <<
268 "I", <<
269 "Xe", <<
270 "Cs", <<
271 "Ba", <<
272 "La", <<
273 "Ce", <<
274 "Pr", <<
275 "Nd", <<
276 "Pm", <<
277 "Sm", <<
278 "Eu", <<
279 "Gd", <<
280 "Tb", <<
281 "Dy", <<
282 "Ho", <<
283 "Er", <<
284 "Tm", <<
285 "Yb", <<
286 "Lu", <<
287 "Hf", <<
288 "Ta", <<
289 "W", <<
290 "Re", <<
291 "Os", <<
292 "Ir", <<
293 "Pt", <<
294 "Au", <<
295 "Hg", <<
296 "Tl", <<
297 "Pb", <<
298 "Bi", <<
299 "Po", <<
300 "At", <<
301 "Rn", <<
302 "Fr", <<
303 "Ra", <<
304 "Ac", <<
305 "Th", <<
306 "Pa", <<
307 "U", <<
308 "Np", <<
309 "Pu", <<
310 "Am", <<
311 "Cm", <<
312 "Bk", <<
313 "Cf", <<
314 "Es", <<
315 "Fm", <<
316 "Md", <<
317 "No", <<
318 "Lr", <<
319 "Rf", <<
320 "Db", <<
321 "Sg", <<
322 "Bh", <<
323 "Hs", <<
324 "Mt", <<
325 "Ds", <<
326 "Rg", <<
327 "Cn" <<
328 }; <<
329 <<
330 /// \brief Digit names to compose IUPAC <<
331 const std::string elementIUPACDigits = " <<
332 <<
333 #define INCL_DEFAULT_SEPARATION_ENERGY 6.83 <<
334 const G4double theINCLProtonSeparationEn <<
335 const G4double theINCLNeutronSeparationE <<
336 const G4double theINCLLambdaSeparationEn <<
337 //const G4double theINCLantiProtonSepara <<
338 const G4double theINCLantiProtonSeparati <<
339 G4ThreadLocal G4double protonSeparationE <<
340 G4ThreadLocal G4double neutronSeparation <<
341 G4ThreadLocal G4double lambdaSeparationE <<
342 //G4ThreadLocal G4double antiprotonSepar <<
343 //G4ThreadLocal G4double antiprotonSepar <<
344 #undef INCL_DEFAULT_SEPARATION_ENERGY <<
345 <<
346 G4ThreadLocal G4double rpCorrelationCoef <<
347 53
348 G4ThreadLocal G4double neutronSkin = 0.0 << 54 /// \brief Static pointer to the mass function for nuclei
349 G4ThreadLocal G4double neutronHalo = 0.0 << 55 ParticleTable::NuclearMassFn ParticleTable::getTableMass;
>> 56 /// \brief Static pointer to the mass function for particles
>> 57 ParticleTable::ParticleMassFn ParticleTable::getTableParticleMass;
>> 58 /// \brief Static pointer to the separation-energy function
>> 59 ParticleTable::SeparationEnergyFn ParticleTable::getSeparationEnergy;
>> 60
>> 61 const G4double ParticleTable::theINCLNucleonMass = 938.2796;
>> 62 const G4double ParticleTable::theINCLPionMass = 138.0;
>> 63 G4double ParticleTable::protonMass = 0.0;
>> 64 G4double ParticleTable::neutronMass = 0.0;
>> 65 G4double ParticleTable::piPlusMass = 0.0;
>> 66 G4double ParticleTable::piMinusMass = 0.0;
>> 67 G4double ParticleTable::piZeroMass = 0.0;
>> 68
>> 69 // e^2/(4 pi epsilon_0) [MeV fm]
>> 70 const G4double ParticleTable::eSquared = 1.439964;
>> 71
>> 72 // Hard-coded values of the real particle masses (MeV/c^2)
>> 73 G4double ParticleTable::theRealProtonMass = 938.27203;
>> 74 G4double ParticleTable::theRealNeutronMass = 939.56536;
>> 75 G4double ParticleTable::theRealChargedPiMass = 139.57018;
>> 76 G4double ParticleTable::theRealPiZeroMass = 134.9766;
>> 77
>> 78 const G4double ParticleTable::mediumDiffuseness[mediumNucleiTableSize] =
>> 79 {0.0,0.0,0.0,0.0,0.0,1.78,1.77,1.77,1.77,1.71,
>> 80 1.69,1.69,1.635,1.730,1.81,1.833,1.798,
>> 81 1.841,0.567,0.571, 0.560,0.549,0.550,0.551,
>> 82 0.580,0.575,0.569,0.537,0.0,0.0};
>> 83 const G4double ParticleTable::mediumRadius[mediumNucleiTableSize] =
>> 84 {0.0,0.0,0.0,0.0,0.0,0.334,0.327,0.479,0.631,0.838,
>> 85 0.811,1.07,1.403,1.335,1.25,1.544,1.498,1.513,
>> 86 2.58,2.77, 2.775,2.78,2.88,2.98,3.22,3.03,2.84,
>> 87 3.14,0.0,0.0};
>> 88 const G4double ParticleTable::positionRMS[clusterTableZSize][clusterTableASize] = {
>> 89 /* A= 0 1 2 3 4 5 6 7 8 9 10 11 12 */
>> 90 /* Z=0 */ {0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00},
>> 91 /* Z=1 */ {0.00, 0.00, 2.10, 1.80, 1.70, 1.83, 2.60, 2.50, 0.00, 0.00, 0.00, 0.00, 0.00},
>> 92 /* Z=2 */ {0.00, 0.00, 0.00, 1.80, 1.68, 1.70, 2.60, 2.50, 2.50, 2.50, 2.50, 0.00, 0.00},
>> 93 /* Z=3 */ {0.00, 0.00, 0.00, 0.00, 1.70, 1.83, 2.56, 2.40, 2.50, 2.50, 2.50, 2.50, 2.50},
>> 94 /* Z=4 */ {0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 2.60, 2.50, 2.50, 2.51, 2.50, 2.50, 2.50},
>> 95 /* Z=5 */ {0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 2.50, 2.50, 2.50, 2.50, 2.45, 2.40, 2.50},
>> 96 /* Z=6 */ {0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 2.50, 2.50, 2.50, 2.50, 2.47},
>> 97 /* Z=7 */ {0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 2.50, 2.50, 2.50},
>> 98 /* Z=8 */ {0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 2.50}
>> 99 };
>> 100
>> 101 const G4double ParticleTable::momentumRMS[clusterTableZSize][clusterTableASize] = {
>> 102 /* A= 0 1 2 3 4 5 6 7 8 9 10 11 12 */
>> 103 /* Z=0 */ {0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00},
>> 104 /* Z=1 */ {0.00, 0.00, 77.0, 110., 153., 100., 100., 100., 0.00, 0.00, 0.00, 0.00, 0.00},
>> 105 /* Z=2 */ {0.00, 0.00, 0.00, 110., 153., 100., 100., 100., 100., 100., 100., 0.00, 0.00},
>> 106 /* Z=3 */ {0.00, 0.00, 0.00, 0.00, 153., 100., 100., 100., 100., 100., 100., 100., 100.},
>> 107 /* Z=4 */ {0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 100., 100., 100., 100., 100., 100., 100.},
>> 108 /* Z=5 */ {0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 100., 100., 100., 100., 100., 100., 100.},
>> 109 /* Z=6 */ {0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 100., 100., 100., 100., 100.},
>> 110 /* Z=7 */ {0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 100., 100., 100.},
>> 111 /* Z=8 */ {0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 100.}
>> 112 };
>> 113
>> 114 /// \brief Table of chemical element names
>> 115 const std::string ParticleTable::elementTable[elementTableSize] = {
>> 116 "",
>> 117 "H",
>> 118 "He",
>> 119 "Li",
>> 120 "Be",
>> 121 "B",
>> 122 "C",
>> 123 "N",
>> 124 "O",
>> 125 "F",
>> 126 "Ne",
>> 127 "Na",
>> 128 "Mg",
>> 129 "Al",
>> 130 "Si",
>> 131 "P",
>> 132 "S",
>> 133 "Cl",
>> 134 "Ar",
>> 135 "K",
>> 136 "Ca",
>> 137 "Sc",
>> 138 "Ti",
>> 139 "V",
>> 140 "Cr",
>> 141 "Mn",
>> 142 "Fe",
>> 143 "Co",
>> 144 "Ni",
>> 145 "Cu",
>> 146 "Zn",
>> 147 "Ga",
>> 148 "Ge",
>> 149 "As",
>> 150 "Se",
>> 151 "Br",
>> 152 "Kr",
>> 153 "Rb",
>> 154 "Sr",
>> 155 "Y",
>> 156 "Zr",
>> 157 "Nb",
>> 158 "Mo",
>> 159 "Tc",
>> 160 "Ru",
>> 161 "Rh",
>> 162 "Pd",
>> 163 "Ag",
>> 164 "Cd",
>> 165 "In",
>> 166 "Sn",
>> 167 "Sb",
>> 168 "Te",
>> 169 "I",
>> 170 "Xe",
>> 171 "Cs",
>> 172 "Ba",
>> 173 "La",
>> 174 "Ce",
>> 175 "Pr",
>> 176 "Nd",
>> 177 "Pm",
>> 178 "Sm",
>> 179 "Eu",
>> 180 "Gd",
>> 181 "Tb",
>> 182 "Dy",
>> 183 "Ho",
>> 184 "Er",
>> 185 "Tm",
>> 186 "Yb",
>> 187 "Lu",
>> 188 "Hf",
>> 189 "Ta",
>> 190 "W",
>> 191 "Re",
>> 192 "Os",
>> 193 "Ir",
>> 194 "Pt",
>> 195 "Au",
>> 196 "Hg",
>> 197 "Tl",
>> 198 "Pb",
>> 199 "Bi",
>> 200 "Po",
>> 201 "At",
>> 202 "Rn",
>> 203 "Fr",
>> 204 "Ra",
>> 205 "Ac",
>> 206 "Th",
>> 207 "Pa",
>> 208 "U",
>> 209 "Np",
>> 210 "Pu",
>> 211 "Am",
>> 212 "Cm",
>> 213 "Bk",
>> 214 "Cf",
>> 215 "Es",
>> 216 "Fm",
>> 217 "Md",
>> 218 "No",
>> 219 "Lr",
>> 220 "Rf",
>> 221 "Db",
>> 222 "Sg",
>> 223 "Bh",
>> 224 "Hs",
>> 225 "Mt",
>> 226 "Ds",
>> 227 "Rg",
>> 228 "Cn"
>> 229 };
>> 230
>> 231 /// \brief Digit names to compose IUPAC element names
>> 232 const std::string ParticleTable::elementIUPACDigits = "nubtqphsoe";
>> 233
>> 234 /* Table for cluster decays
>> 235 * Definition of "Stable": halflife > 1 ms
>> 236 *
>> 237 * These table includes decay data for clusters that INCL presently does
>> 238 * not produce. It can't hurt.
>> 239 *
>> 240 * Unphysical nuclides (A<Z) are marked as stable, but should never be
>> 241 * produced by INCL. If you find them in the output, something is fishy.
>> 242 */
>> 243 const ParticleTable::ClusterDecayType ParticleTable::clusterDecayMode[clusterTableZSize][clusterTableASize] =
>> 244 {
>> 245 /* A = 0 1 2 3 4 5 6 7 8 9 10 11 12 */
>> 246 /* Z = 0 */ {StableCluster, StableCluster, NeutronDecay, NeutronUnbound, NeutronUnbound, NeutronUnbound, NeutronUnbound, NeutronUnbound, NeutronUnbound, NeutronUnbound, NeutronUnbound, NeutronUnbound, NeutronUnbound},
>> 247 /* Z = 1 */ {StableCluster, StableCluster, StableCluster, StableCluster, NeutronDecay, TwoNeutronDecay, NeutronDecay, TwoNeutronDecay, NeutronUnbound, NeutronUnbound, NeutronUnbound, NeutronUnbound, NeutronUnbound},
>> 248 /* Z = 2 */ {StableCluster, StableCluster, ProtonDecay, StableCluster, StableCluster, NeutronDecay, StableCluster, NeutronDecay, StableCluster, NeutronDecay, TwoNeutronDecay, NeutronUnbound, NeutronUnbound},
>> 249 /* Z = 3 */ {StableCluster, StableCluster, StableCluster, ProtonUnbound, ProtonDecay, ProtonDecay, StableCluster, StableCluster, StableCluster, StableCluster, NeutronDecay, StableCluster, NeutronDecay},
>> 250 /* Z = 4 */ {StableCluster, StableCluster, StableCluster, StableCluster, ProtonUnbound, ProtonDecay, TwoProtonDecay, StableCluster, AlphaDecay, StableCluster, StableCluster, StableCluster, StableCluster},
>> 251 /* Z = 5 */ {StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, ProtonUnbound, TwoProtonDecay, ProtonDecay, StableCluster, ProtonDecay, StableCluster, StableCluster, StableCluster},
>> 252 /* Z = 6 */ {StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, ProtonUnbound, ProtonUnbound, TwoProtonDecay, StableCluster, StableCluster, StableCluster, StableCluster},
>> 253 /* Z = 7 */ {StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, ProtonUnbound, ProtonUnbound, ProtonUnbound, ProtonDecay, ProtonDecay, StableCluster},
>> 254 /* Z = 8 */ {StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, ProtonUnbound, ProtonUnbound, ProtonUnbound, ProtonUnbound, ProtonDecay}
>> 255 };
>> 256
>> 257 /**
>> 258 * Precomputed factor 1.0/A
>> 259 */
>> 260 const G4double ParticleTable::clusterPosFact[maxClusterMass+1] = {0.0, 1.0, 0.5,
>> 261 0.33333, 0.25,
>> 262 0.2, 0.16667,
>> 263 0.14286, 0.125,
>> 264 0.11111, 0.1,
>> 265 0.09091, 0.083333};
>> 266
>> 267 /**
>> 268 * Precomputed factor (1.0/A)^2
>> 269 */
>> 270 const G4double ParticleTable::clusterPosFact2[maxClusterMass+1] = {0.0, 1.0, 0.25,
>> 271 0.11111, 0.0625,
>> 272 0.04, 0.0277778,
>> 273 0.020408, 0.015625,
>> 274 0.012346, 0.01,
>> 275 0.0082645, 0.0069444};
>> 276
>> 277 const G4int ParticleTable::clusterZMin[maxClusterMass+1] = {0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 3, 3};
>> 278 const G4int ParticleTable::clusterZMax[maxClusterMass+1] = {0, 0, 1, 2, 3, 3, 5, 5, 6, 6, 7, 7, 8};
>> 279 const G4double ParticleTable::clusterPhaseSpaceCut[maxClusterMass+1] = {0.0, 70000.0, 180000.0,
>> 280 90000.0, 90000.0,
>> 281 128941.0 ,145607.0,
>> 282 161365.0, 176389.0,
>> 283 190798.0, 204681.0,
>> 284 218109.0, 231135.0};
>> 285 const G4double ParticleTable::effectiveNucleonMass = 938.2796;
>> 286 const G4double ParticleTable::effectiveNucleonMass2 = 8.8036860777616e5;
>> 287 const G4double ParticleTable::effectiveDeltaMass = 1232.0;
>> 288 const G4double ParticleTable::effectivePionMass = 138.0;
>> 289 const G4double ParticleTable::effectiveDeltaDecayThreshold =
>> 290 ParticleTable::theRealNeutronMass + ParticleTable::theRealChargedPiMass
>> 291 + 0.5;
>> 292 const G4double ParticleTable::theINCLProtonSeparationEnergy = 6.83;
>> 293 const G4double ParticleTable::theINCLNeutronSeparationEnergy = ParticleTable::theINCLProtonSeparationEnergy;
>> 294 G4double ParticleTable::protonSeparationEnergy = theINCLProtonSeparationEnergy;
>> 295 G4double ParticleTable::neutronSeparationEnergy = theINCLNeutronSeparationEnergy;
350 296
351 #ifdef INCLXX_IN_GEANT4_MODE 297 #ifdef INCLXX_IN_GEANT4_MODE
352 G4ThreadLocal G4IonTable *theG4IonTable; << 298 G4IonTable *ParticleTable::theG4IonTable;
>> 299 #else
>> 300 std::vector< std::vector <G4bool> > ParticleTable::massTableMask;
>> 301 std::vector< std::vector <G4double> > ParticleTable::massTable;
353 #endif 302 #endif
354 303
355 /// \brief Default value for constant Fe << 304 void ParticleTable::initialize(Config const * const theConfig /*=0*/) {
356 G4ThreadLocal G4double constantFermiMome << 305 protonMass = theINCLNucleonMass;
357 << 306 neutronMass = theINCLNucleonMass;
358 /// \brief Transform a IUPAC char to an << 307 piPlusMass = theINCLPionMass;
359 char iupacToInt(char c) { << 308 piMinusMass = theINCLPionMass;
360 return (char)(((G4int)'0')+elementIUPA << 309 piZeroMass = theINCLPionMass;
361 } <<
362 <<
363 /// \brief Transform an integer digit (r <<
364 char intToIUPAC(char n) { return element <<
365 <<
366 /// \brief Get the singleton instance of <<
367 const NaturalIsotopicDistributions *getN <<
368 if(!theNaturalIsotopicDistributions) <<
369 theNaturalIsotopicDistributions = ne <<
370 return theNaturalIsotopicDistributions <<
371 } <<
372 <<
373 } // namespace <<
374 <<
375 void initialize(Config const * const theCo <<
376 protonMass = theINCLNucleonMass; <<
377 neutronMass = theINCLNucleonMass; <<
378 piPlusMass = theINCLPionMass; <<
379 piMinusMass = theINCLPionMass; <<
380 piZeroMass = theINCLPionMass; <<
381 <<
382 etaMass = theINCLEtaMass; <<
383 omegaMass = theINCLOmegaMass; <<
384 etaPrimeMass = theINCLEtaPrimeMass; <<
385 photonMass = theINCLPhotonMass; <<
386 <<
387 SigmaPlusMass = theRealSigmaPlusMass; <<
388 SigmaMinusMass = theRealSigmaMinusMass; <<
389 SigmaZeroMass = theRealSigmaZeroMass; <<
390 LambdaMass = theINCLLambdaMass; <<
391 KPlusMass = theRealChargedKaonMass; <<
392 KZeroMass = theRealNeutralKaonMass; <<
393 KZeroBarMass = theRealNeutralKaonMass; <<
394 KShortMass = theRealNeutralKaonMass; <<
395 KLongMass = theRealNeutralKaonMass; <<
396 KMinusMass = theRealChargedKaonMass; <<
397 <<
398 antiProtonMass = theRealAntiProtonMass; <<
399 XiZeroMass = theRealXiZeroMass; <<
400 XiMinusMass = theRealXiMinusMass; <<
401 antiNeutronMass = theRealAntiNeutronMass <<
402 antiSigmaPlusMass = theRealAntiSigmaPlus <<
403 antiSigmaMinusMass = theRealAntiSigmaMin <<
404 antiSigmaZeroMass = theRealAntiSigmaZero <<
405 antiLambdaMass = theRealAntiLambdaMass; <<
406 antiXiZeroMass = theRealAntiXiZeroMass; <<
407 antiXiMinusMass = theRealAntiXiMinusMass <<
408 <<
409 if(theConfig && theConfig->getUseRealMas <<
410 getTableMass = getRealMass; <<
411 getTableParticleMass = getRealMass; <<
412 } else { <<
413 getTableMass = getINCLMass; <<
414 getTableParticleMass = getINCLMass; <<
415 } <<
416 310
417 #ifndef INCLXX_IN_GEANT4_MODE 311 #ifndef INCLXX_IN_GEANT4_MODE
418 std::string dataFilePath; << 312 std::string dataFilePath;
419 if(theConfig) << 313 if(theConfig)
420 dataFilePath = theConfig->getINCLXXDat << 314 dataFilePath = theConfig->getINCLXXDataFilePath();
421 NuclearMassTable::initialize(dataFilePat << 315 readRealMasses(dataFilePath);
422 #endif <<
423 <<
424 #ifdef INCLXX_IN_GEANT4_MODE <<
425 G4ParticleTable *theG4ParticleTable = G4 <<
426 theG4IonTable = theG4ParticleTable->GetI <<
427 theRealProtonMass = theG4ParticleTable-> <<
428 theRealNeutronMass = theG4ParticleTable- <<
429 theRealChargedPiMass = theG4ParticleTabl <<
430 theRealPiZeroMass = theG4ParticleTable-> <<
431 <<
432 theRealEtaMass = theG4ParticleTable->Fin <<
433 theRealOmegaMass = theG4ParticleTable->F <<
434 theRealEtaPrimeMass = theG4ParticleTable <<
435 theRealPhotonMass = theG4ParticleTable-> <<
436 <<
437 theRealSigmaPlusMass = theG4ParticleTabl <<
438 theRealSigmaZeroMass = theG4ParticleTabl <<
439 theRealSigmaMinusMass = theG4ParticleTab <<
440 theRealLambdaMass = theG4ParticleTable-> <<
441 theRealChargedKaonMass = theG4ParticleTa <<
442 theRealNeutralKaonMass = theG4ParticleTa <<
443 <<
444 theRealAntiProtonMass = theG4ParticleTab <<
445 theRealAntiNeutronMass = theG4ParticleTa <<
446 theRealXiZeroMass = theG4ParticleTable-> <<
447 theRealXiMinusMass = theG4ParticleTable- <<
448 theRealAntiSigmaPlusMass = theG4Particle <<
449 theRealAntiSigmaZeroMass = theG4Particle <<
450 theRealAntiSigmaMinusMass = theG4Particl <<
451 theRealAntiLambdaMass = theG4ParticleTab <<
452 theRealAntiXiZeroMass = theG4ParticleTab <<
453 theRealAntiXiMinusMass = theG4ParticleTa <<
454 #endif 316 #endif
455 317
456 minDeltaMass = theRealNeutronMass + theR << 318 if(theConfig && theConfig->getUseRealMasses()) {
457 minDeltaMass2 = minDeltaMass*minDeltaMas << 319 getTableMass = ParticleTable::getRealMass;
458 minDeltaMassRndm = std::atan((minDeltaMa << 320 getTableParticleMass = ParticleTable::getRealMass;
459 << 321 } else {
460 piPlusWidth = theChargedPiWidth; << 322 getTableMass = ParticleTable::getINCLMass;
461 piMinusWidth = theChargedPiWidth; << 323 getTableParticleMass = ParticleTable::getINCLMass;
462 piZeroWidth = thePiZeroWidth; << 324 }
463 etaWidth = theEtaWidth; <<
464 omegaWidth = theOmegaWidth; <<
465 etaPrimeWidth = theEtaPrimeWidth; <<
466 <<
467 SigmaMinusWidth = theSigmaMinusWidth; <<
468 SigmaPlusWidth = theSigmaPlusWidth; <<
469 SigmaZeroWidth = theSigmaZeroWidth; <<
470 LambdaWidth = theLambdaWidth; <<
471 KPlusWidth = theChargedKaonWidth; <<
472 KMinusWidth = theChargedKaonWidth; <<
473 KShortWidth = theKShortWidth; <<
474 KLongWidth = theKLongWidth; <<
475 <<
476 // Initialise HFB tables <<
477 #ifdef INCLXX_IN_GEANT4_MODE 325 #ifdef INCLXX_IN_GEANT4_MODE
478 HFB::initialize(); << 326 G4ParticleTable *theG4ParticleTable = G4ParticleTable::GetParticleTable();
479 #else << 327 theG4IonTable = theG4ParticleTable->GetIonTable();
480 HFB::initialize(dataFilePath); << 328 theRealProtonMass = theG4ParticleTable->FindParticle("proton")->GetPDGMass() / MeV;
>> 329 theRealNeutronMass = theG4ParticleTable->FindParticle("neutron")->GetPDGMass() / MeV;
>> 330 theRealChargedPiMass = theG4ParticleTable->FindParticle("pi+")->GetPDGMass() / MeV;
>> 331 theRealPiZeroMass = theG4ParticleTable->FindParticle("pi0")->GetPDGMass() / MeV;
481 #endif 332 #endif
482 333
483 // Initialise the separation-energy func << 334 // Initialise the separation-energy function
484 if(!theConfig || theConfig->getSeparatio << 335 if(!theConfig || theConfig->getSeparationEnergyType()==INCLSeparationEnergy)
485 getSeparationEnergy = getSeparationEne << 336 getSeparationEnergy = ParticleTable::getSeparationEnergyINCL;
486 else if(theConfig->getSeparationEnergyTy << 337 else if(theConfig->getSeparationEnergyType()==RealSeparationEnergy)
487 getSeparationEnergy = getSeparationEne << 338 getSeparationEnergy = ParticleTable::getSeparationEnergyReal;
488 else if(theConfig->getSeparationEnergyTy << 339 else if(theConfig->getSeparationEnergyType()==RealForLightSeparationEnergy)
489 getSeparationEnergy = getSeparationEne << 340 getSeparationEnergy = ParticleTable::getSeparationEnergyRealForLight;
490 else { << 341 else {
491 INCL_FATAL("Unrecognized separation-en << 342 FATAL("Unrecognized separation-energy type in ParticleTable initialization: " << theConfig->getSeparationEnergyType() << std::endl);
492 return; << 343 }
493 } << 344
494 << 345 }
495 // Initialise the Fermi-momentum functio << 346
496 if(!theConfig || theConfig->getFermiMome << 347 G4int ParticleTable::getIsospin(const ParticleType t) {
497 getFermiMomentum = ParticleTable::getF << 348 // Actually this is the 3rd component of isospin (I_z) multiplied by 2!
498 if(theConfig) { << 349 if(t == Proton) {
499 const G4double aFermiMomentum = theC << 350 return 1;
500 if(aFermiMomentum>0.) << 351 } else if(t == Neutron) {
501 constantFermiMomentum = aFermiMome << 352 return -1;
502 else << 353 } else if(t == PiPlus) {
503 constantFermiMomentum = PhysicalCo << 354 return 2;
504 } else { << 355 } else if(t == PiMinus) {
505 constantFermiMomentum = PhysicalCons << 356 return -2;
506 } << 357 } else if(t == PiZero) {
507 } else if(theConfig->getFermiMomentumTyp << 358 return 0;
508 getFermiMomentum = ParticleTable::getF << 359 } else if(t == DeltaPlusPlus) {
509 else if(theConfig->getFermiMomentumType( << 360 return 3;
510 getFermiMomentum = ParticleTable::getF << 361 } else if(t == DeltaPlus) {
511 else { << 362 return 1;
512 INCL_FATAL("Unrecognized Fermi-momentu << 363 } else if(t == DeltaZero) {
513 return; << 364 return -1;
514 } << 365 } else if(t == DeltaMinus) {
515 << 366 return -3;
516 // Initialise the r-p correlation coeffi << 367 }
517 std::fill(rpCorrelationCoefficient, rpCo << 368
518 if(theConfig) { << 369 ERROR("Requested isospin of an unknown particle!");
519 rpCorrelationCoefficient[Proton] = the << 370 return -10; // Unknown
520 rpCorrelationCoefficient[Neutron] = th << 371 }
521 } << 372
522 << 373 std::string ParticleTable::getShortName(const ParticleSpecies s) {
523 // Initialise the neutron-skin parameter << 374 if(s.theType==Composite)
524 if(theConfig) { << 375 return getShortName(s.theA,s.theZ);
525 neutronSkin = theConfig->getNeutronSki << 376 else
526 neutronHalo = theConfig->getNeutronHal << 377 return getShortName(s.theType);
527 } << 378 }
528 << 379
529 } << 380 std::string ParticleTable::getName(const ParticleSpecies s) {
530 << 381 if(s.theType==Composite)
531 G4int getIsospin(const ParticleType t) { << 382 return getName(s.theA,s.theZ);
532 // Actually this is the 3rd component of << 383 else
533 if(t == Proton) { << 384 return getName(s.theType);
534 return 1; << 385 }
535 } else if(t == Neutron) { << 386
536 return -1; << 387 std::string ParticleTable::getName(const G4int A, const G4int Z) {
537 } else if(t == PiPlus) { << 388 std::stringstream stream;
538 return 2; << 389 stream << getElementName(Z) << "-" << A;
539 } else if(t == PiMinus) { << 390 return stream.str();
540 return -2; << 391 }
541 } else if(t == PiZero) { << 392
542 return 0; << 393 std::string ParticleTable::getShortName(const G4int A, const G4int Z) {
543 } else if(t == DeltaPlusPlus) { << 394 std::stringstream stream;
544 return 3; << 395 stream << getElementName(Z);
545 } else if(t == DeltaPlus) { << 396 if(A>0)
546 return 1; << 397 stream << A;
547 } else if(t == DeltaZero) { << 398 return stream.str();
548 return -1; << 399 }
549 } else if(t == DeltaMinus) { << 400
550 return -3; << 401 std::string ParticleTable::getName(const ParticleType p) {
551 } else if(t == Lambda) { << 402 if(p == G4INCL::Proton) {
552 return 0; << 403 return std::string("proton");
553 } else if(t == SigmaPlus) { << 404 } else if(p == G4INCL::Neutron) {
554 return 2; << 405 return std::string("neutron");
555 } else if(t == SigmaZero) { << 406 } else if(p == G4INCL::DeltaPlusPlus) {
556 return 0; << 407 return std::string("delta++");
557 } else if(t == SigmaMinus) { << 408 } else if(p == G4INCL::DeltaPlus) {
558 return -2; << 409 return std::string("delta+");
559 } else if(t == KPlus) { << 410 } else if(p == G4INCL::DeltaZero) {
560 return 1; << 411 return std::string("delta0");
561 } else if(t == KZero) { << 412 } else if(p == G4INCL::DeltaMinus) {
562 return -1; << 413 return std::string("delta-");
563 } else if(t == KZeroBar) { << 414 } else if(p == G4INCL::PiPlus) {
564 return 1; << 415 return std::string("pi+");
565 } else if(t == KShort) { << 416 } else if(p == G4INCL::PiZero) {
566 return 0; << 417 return std::string("pi0");
567 } else if(t == KLong) { << 418 } else if(p == G4INCL::PiMinus) {
568 return 0; << 419 return std::string("pi-");
569 } else if(t == KMinus) { << 420 } else if(p == G4INCL::Composite) {
570 return -1; << 421 return std::string("composite");
571 } else if(t == Eta) { << 422 }
572 return 0; << 423 return std::string("unknown");
573 } else if(t == Omega) { << 424 }
574 return 0; << 425
575 } else if(t == EtaPrime) { << 426 std::string ParticleTable::getShortName(const ParticleType p) {
576 return 0; << 427 if(p == G4INCL::Proton) {
577 } else if(t == Photon) { << 428 return std::string("p");
578 return 0; << 429 } else if(p == G4INCL::Neutron) {
579 } else if(t == antiProton) { << 430 return std::string("n");
580 return -1; << 431 } else if(p == G4INCL::DeltaPlusPlus) {
581 } else if(t == XiMinus) { << 432 return std::string("d++");
582 return -1; << 433 } else if(p == G4INCL::DeltaPlus) {
583 } else if(t == XiZero) { << 434 return std::string("d+");
584 return 1; << 435 } else if(p == G4INCL::DeltaZero) {
585 } else if(t == antiNeutron) { << 436 return std::string("d0");
586 return 1; << 437 } else if(p == G4INCL::DeltaMinus) {
587 } else if(t == antiLambda) { << 438 return std::string("d-");
588 return 0; << 439 } else if(p == G4INCL::PiPlus) {
589 } else if(t == antiSigmaPlus) { << 440 return std::string("pi+");
590 return -2; << 441 } else if(p == G4INCL::PiZero) {
591 } else if(t == antiSigmaZero) { << 442 return std::string("pi0");
592 return 0; << 443 } else if(p == G4INCL::PiMinus) {
593 } else if(t == antiSigmaMinus) { << 444 return std::string("pi-");
594 return 2; << 445 } else if(p == G4INCL::Composite) {
595 } else if(t == antiXiMinus) { << 446 return std::string("comp");
596 return 1; << 447 }
597 } else if(t == antiXiZero) { << 448 return std::string("unknown");
598 return -1; << 449 }
599 } << 450
600 INCL_ERROR("Requested isospin of an unkn << 451 G4double ParticleTable::getINCLMass(const ParticleType pt) {
601 return -10; // Unknown << 452 if(pt == Proton) {
602 } << 453 return protonMass;
603 << 454 } else if(pt == Neutron) {
604 std::string getShortName(const ParticleSpe << 455 return neutronMass;
605 if(sp.theType==Composite && sp.theS == 0 << 456 } else if(pt == PiPlus) {
606 return getShortName(sp.theA,sp.theZ); << 457 return piPlusMass;
607 else if(sp.theType==Composite) << 458 } else if(pt == PiMinus) {
608 return getName(sp.theA,sp.theZ,sp.theS << 459 return piMinusMass;
609 else << 460 } else if(pt == PiZero) {
610 return getShortName(sp.theType); << 461 return piZeroMass;
611 } << 462 } else {
612 << 463 ERROR("ParticleTable::getMass : Unknown particle type." << std::endl);
613 std::string getName(const ParticleSpecies << 464 return 0.0;
614 if(sp.theType==Composite && sp.theS == 0 << 465 }
615 return getName(sp.theA,sp.theZ); << 466 }
616 else if(sp.theType==Composite) << 467
617 return getName(sp.theA,sp.theZ,sp.theS << 468 G4double ParticleTable::getRealMass(const ParticleType t) {
618 else << 469 switch(t) {
619 return getName(sp.theType); << 470 case Proton:
620 } << 471 return theRealProtonMass;
621 << 472 break;
622 std::string getName(const G4int A, const G << 473 case Neutron:
623 std::stringstream stream; << 474 return theRealNeutronMass;
624 stream << getElementName(Z) << "-" << A; << 475 break;
625 return stream.str(); << 476 case PiPlus:
626 } << 477 case PiMinus:
627 << 478 return theRealChargedPiMass;
628 std::string getName(const G4int A, const G << 479 break;
629 std::stringstream stream; << 480 case PiZero:
630 if(S >= 0) // S < 0 for hypernuclei << 481 return theRealPiZeroMass;
631 return getName(A, Z); << 482 break;
632 else if(S == -1) << 483 default:
633 stream << getElementName(Z) << "-" << << 484 ERROR("Particle::getRealMass : Unknown particle type." << std::endl);
634 else <<
635 stream << getElementName(Z) << "-" << <<
636 return stream.str(); <<
637 } <<
638 <<
639 std::string getShortName(const G4int A, co <<
640 std::stringstream stream; <<
641 stream << getElementName(Z); <<
642 if(A>0) <<
643 stream << A; <<
644 return stream.str(); <<
645 } <<
646 <<
647 std::string getName(const ParticleType p) <<
648 if(p == G4INCL::Proton) { <<
649 return std::string("proton"); <<
650 } else if(p == G4INCL::Neutron) { <<
651 return std::string("neutron"); <<
652 } else if(p == G4INCL::DeltaPlusPlus) { <<
653 return std::string("delta++"); <<
654 } else if(p == G4INCL::DeltaPlus) { <<
655 return std::string("delta+"); <<
656 } else if(p == G4INCL::DeltaZero) { <<
657 return std::string("delta0"); <<
658 } else if(p == G4INCL::DeltaMinus) { <<
659 return std::string("delta-"); <<
660 } else if(p == G4INCL::PiPlus) { <<
661 return std::string("pi+"); <<
662 } else if(p == G4INCL::PiZero) { <<
663 return std::string("pi0"); <<
664 } else if(p == G4INCL::PiMinus) { <<
665 return std::string("pi-"); <<
666 } else if(p == G4INCL::Lambda) { <<
667 return std::string("lambda"); <<
668 } else if(p == G4INCL::SigmaPlus) { <<
669 return std::string("sigma+"); <<
670 } else if(p == G4INCL::SigmaZero) { <<
671 return std::string("sigma0"); <<
672 } else if(p == G4INCL::SigmaMinus) { <<
673 return std::string("sigma-"); <<
674 } else if(p == G4INCL::antiProton) { <<
675 return std::string("antiproton"); <<
676 } else if(p == G4INCL::XiMinus) { <<
677 return std::string("xi-"); <<
678 } else if(p == G4INCL::XiZero) { <<
679 return std::string("xi0"); <<
680 } else if(p == G4INCL::antiNeutron) { <<
681 return std::string("antineutron"); <<
682 } else if(p == G4INCL::antiSigmaPlus) { <<
683 return std::string("antisigma+"); <<
684 } else if(p == G4INCL::antiSigmaZero) { <<
685 return std::string("antisigma0"); <<
686 } else if(p == G4INCL::antiSigmaMinus) { <<
687 return std::string("antisigma-"); <<
688 } else if(p == G4INCL::antiLambda) { <<
689 return std::string("antilambda"); <<
690 } else if(p == G4INCL::antiXiMinus) { <<
691 return std::string("antixi-"); <<
692 } else if(p == G4INCL::antiXiZero) { <<
693 return std::string("antixi0"); <<
694 } else if(p == G4INCL::KPlus) { <<
695 return std::string("kaon+"); <<
696 } else if(p == G4INCL::KZero) { <<
697 return std::string("kaon0"); <<
698 } else if(p == G4INCL::KZeroBar) { <<
699 return std::string("kaon0bar"); <<
700 } else if(p == G4INCL::KMinus) { <<
701 return std::string("kaon-"); <<
702 } else if(p == G4INCL::KShort) { <<
703 return std::string("kaonshort"); <<
704 } else if(p == G4INCL::KLong) { <<
705 return std::string("kaonlong"); <<
706 } else if(p == G4INCL::Composite) { <<
707 return std::string("composite"); <<
708 } else if(p == G4INCL::Eta) { <<
709 return std::string("eta"); <<
710 } else if(p == G4INCL::Omega) { <<
711 return std::string("omega"); <<
712 } else if(p == G4INCL::EtaPrime) { <<
713 return std::string("etaprime"); <<
714 } else if(p == G4INCL::Photon) { <<
715 return std::string("photon"); <<
716 } <<
717 return std::string("unknown"); <<
718 } <<
719 <<
720 std::string getShortName(const ParticleTyp <<
721 if(p == G4INCL::Proton) { <<
722 return std::string("p"); <<
723 } else if(p == G4INCL::Neutron) { <<
724 return std::string("n"); <<
725 } else if(p == G4INCL::DeltaPlusPlus) { <<
726 return std::string("d++"); <<
727 } else if(p == G4INCL::DeltaPlus) { <<
728 return std::string("d+"); <<
729 } else if(p == G4INCL::DeltaZero) { <<
730 return std::string("d0"); <<
731 } else if(p == G4INCL::DeltaMinus) { <<
732 return std::string("d-"); <<
733 } else if(p == G4INCL::PiPlus) { <<
734 return std::string("pi+"); <<
735 } else if(p == G4INCL::PiZero) { <<
736 return std::string("pi0"); <<
737 } else if(p == G4INCL::PiMinus) { <<
738 return std::string("pi-"); <<
739 } else if(p == G4INCL::Lambda) { <<
740 return std::string("l"); <<
741 } else if(p == G4INCL::SigmaPlus) { <<
742 return std::string("s+"); <<
743 } else if(p == G4INCL::SigmaZero) { <<
744 return std::string("s0"); <<
745 } else if(p == G4INCL::SigmaMinus) { <<
746 return std::string("s-"); <<
747 } else if(p == G4INCL::antiProton) { <<
748 return std::string("pb"); <<
749 } else if(p == G4INCL::XiMinus) { <<
750 return std::string("x-"); <<
751 } else if(p == G4INCL::XiZero) { <<
752 return std::string("x0"); <<
753 } else if(p == G4INCL::antiNeutron) { <<
754 return std::string("nb"); <<
755 } else if(p == G4INCL::antiSigmaPlus) { <<
756 return std::string("s+b"); <<
757 } else if(p == G4INCL::antiSigmaZero) { <<
758 return std::string("s0b"); <<
759 } else if(p == G4INCL::antiSigmaMinus) { <<
760 return std::string("s-b"); <<
761 } else if(p == G4INCL::antiLambda) { <<
762 return std::string("lb"); <<
763 } else if(p == G4INCL::antiXiMinus) { <<
764 return std::string("x-b"); <<
765 } else if(p == G4INCL::antiXiZero) { <<
766 return std::string("x0b"); <<
767 } else if(p == G4INCL::KPlus) { <<
768 return std::string("k+"); <<
769 } else if(p == G4INCL::KZero) { <<
770 return std::string("k0"); <<
771 } else if(p == G4INCL::KZeroBar) { <<
772 return std::string("k0b"); <<
773 } else if(p == G4INCL::KMinus) { <<
774 return std::string("k-"); <<
775 } else if(p == G4INCL::KShort) { <<
776 return std::string("ks"); <<
777 } else if(p == G4INCL::KLong) { <<
778 return std::string("kl"); <<
779 } else if(p == G4INCL::Composite) { <<
780 return std::string("comp"); <<
781 } else if(p == G4INCL::Eta) { <<
782 return std::string("eta"); <<
783 } else if(p == G4INCL::Omega) { <<
784 return std::string("omega"); <<
785 } else if(p == G4INCL::EtaPrime) { <<
786 return std::string("etap"); <<
787 } else if(p == G4INCL::Photon) { <<
788 return std::string("photon"); <<
789 } <<
790 return std::string("unknown"); <<
791 } <<
792 <<
793 G4double getINCLMass(const ParticleType pt <<
794 if(pt == Proton) { <<
795 return protonMass; <<
796 } else if(pt == Neutron) { <<
797 return neutronMass; <<
798 } else if(pt == PiPlus) { <<
799 return piPlusMass; <<
800 } else if(pt == PiMinus) { <<
801 return piMinusMass; <<
802 } else if(pt == PiZero) { <<
803 return piZeroMass; <<
804 } else if(pt == SigmaPlus) { <<
805 return SigmaPlusMass; <<
806 } else if(pt == SigmaMinus) { <<
807 return SigmaMinusMass; <<
808 } else if(pt == SigmaZero) { <<
809 return SigmaZeroMass; <<
810 } else if(pt == Lambda) { <<
811 return LambdaMass; <<
812 } else if(pt == antiProton) { <<
813 return antiProtonMass; <<
814 } else if(pt == XiMinus) { <<
815 return XiMinusMass; <<
816 } else if(pt == XiZero) { <<
817 return XiZeroMass; <<
818 } else if(pt == antiNeutron) { <<
819 return antiNeutronMass; <<
820 } else if(pt == antiSigmaPlus) { <<
821 return antiSigmaPlusMass; <<
822 } else if(pt == antiSigmaMinus) { <<
823 return antiSigmaMinusMass; <<
824 } else if(pt == antiSigmaZero) { <<
825 return antiSigmaZeroMass; <<
826 } else if(pt == antiLambda) { <<
827 return antiLambdaMass; <<
828 } else if(pt == antiXiMinus) { <<
829 return antiXiMinusMass; <<
830 } else if(pt == antiXiZero) { <<
831 return antiXiZeroMass; <<
832 } else if(pt == KPlus) { <<
833 return KPlusMass; <<
834 } else if(pt == KZero) { <<
835 return KZeroMass; <<
836 } else if(pt == KZeroBar) { <<
837 return KZeroBarMass; <<
838 } else if(pt == KMinus) { <<
839 return KMinusMass; <<
840 } else if(pt == KShort) { <<
841 return KShortMass; <<
842 } else if(pt == KLong) { <<
843 return KLongMass; <<
844 } else if(pt == Eta) { <<
845 return etaMass; <<
846 } else if(pt == Omega) { <<
847 return omegaMass; <<
848 } else if(pt == EtaPrime) { <<
849 return etaPrimeMass; <<
850 } else if(pt == Photon) { <<
851 return photonMass; <<
852 } else { <<
853 INCL_ERROR("getMass : Unknown particle <<
854 return 0.0; 485 return 0.0;
855 } << 486 break;
856 } 487 }
857 << 488 }
858 G4double getRealMass(const ParticleType t) << 489
859 switch(t) { << 490 G4double ParticleTable::getRealMass(const G4int A, const G4int Z) {
860 case Proton: <<
861 return theRealProtonMass; <<
862 break; <<
863 case Neutron: <<
864 return theRealNeutronMass; <<
865 break; <<
866 case PiPlus: <<
867 case PiMinus: <<
868 return theRealChargedPiMass; <<
869 break; <<
870 case PiZero: <<
871 return theRealPiZeroMass; <<
872 break; <<
873 case Eta: <<
874 return theRealEtaMass; <<
875 break; <<
876 case Omega: <<
877 return theRealOmegaMass; <<
878 break; <<
879 case EtaPrime: <<
880 return theRealEtaPrimeMass; <<
881 break; <<
882 case Photon: <<
883 return theRealPhotonMass; <<
884 break; <<
885 case Lambda: <<
886 return theRealLambdaMass; <<
887 break; <<
888 case KPlus: <<
889 case KMinus: <<
890 return theRealChargedKaonMass; <<
891 break; <<
892 case KZero: <<
893 case KZeroBar: <<
894 case KShort: <<
895 case KLong: <<
896 return theRealNeutralKaonMass; <<
897 break; <<
898 case SigmaPlus: <<
899 return theRealSigmaPlusMass; <<
900 break; <<
901 case SigmaZero: <<
902 return theRealSigmaZeroMass; <<
903 break; <<
904 case SigmaMinus: <<
905 return theRealSigmaMinusMass; <<
906 break; <<
907 case antiProton: <<
908 return theRealAntiProtonMass; <<
909 break; <<
910 case XiMinus: <<
911 return theRealXiMinusMass; <<
912 break; <<
913 case XiZero: <<
914 return theRealXiZeroMass; <<
915 break; <<
916 case antiNeutron: <<
917 return theRealAntiNeutronMass; <<
918 break; <<
919 case antiSigmaPlus: <<
920 return theRealAntiSigmaPlusMass; <<
921 break; <<
922 case antiSigmaZero: <<
923 return theRealAntiSigmaZeroMass; <<
924 break; <<
925 case antiSigmaMinus: <<
926 return theRealAntiSigmaMinusMass; <<
927 break; <<
928 case antiXiMinus: <<
929 return theRealAntiXiMinusMass; <<
930 break; <<
931 case antiXiZero: <<
932 return theRealAntiXiZeroMass; <<
933 break; <<
934 case antiLambda: <<
935 return theRealAntiLambdaMass; <<
936 break; <<
937 default: <<
938 INCL_ERROR("Particle::getRealMass : <<
939 return 0.0; <<
940 break; <<
941 } <<
942 } <<
943 <<
944 G4double getRealMass(const G4int A, const <<
945 // assert(A>=0); 491 // assert(A>=0);
946 // For nuclei with Z<0 or Z>A, assume th << 492 // For nuclei with Z<0 or Z>A, assume that the exotic charge state is due to pions
947 if(Z<0 && S<0) << 493 if(Z<0)
948 return (A+S)*theRealNeutronMass - S*La << 494 return A*neutronMass - Z*getRealMass(PiMinus);
949 else if(Z>A && S<0) << 495 else if(Z>A)
950 return (A+S)*theRealProtonMass - S*Lam << 496 return A*protonMass + (A-Z)*getRealMass(PiPlus);
951 if(Z<0) << 497 else if(Z==0)
952 return (A)*theRealNeutronMass - Z*getR << 498 return A*getRealMass(Neutron);
953 else if(Z>A) << 499 else if(A==Z)
954 return (A)*theRealProtonMass + (A-Z)*g << 500 return A*getRealMass(Proton);
955 else if(Z==0 && S==0) << 501 else if(A>1) {
956 return A*theRealNeutronMass; <<
957 else if(A==Z) <<
958 return A*theRealProtonMass; <<
959 else if(Z==0 && S<0) <<
960 return (A+S)*theRealNeutronMass-S*Lamb <<
961 else if(A>1) { <<
962 #ifndef INCLXX_IN_GEANT4_MODE 502 #ifndef INCLXX_IN_GEANT4_MODE
963 return ::G4INCL::NuclearMassTable::get << 503 if(ParticleTable::hasMassTable(A,Z))
>> 504 return ParticleTable::massTable.at(Z).at(A);
>> 505 else {
>> 506 DEBUG("Real mass unavailable for isotope A=" << A << ", Z=" << Z
>> 507 << ", using Weizsaecker's formula"
>> 508 << std::endl);
>> 509 return getWeizsaeckerMass(A,Z);
>> 510 }
964 #else 511 #else
965 if(S<0) return theG4IonTable->GetNucle << 512 return theG4IonTable->GetNucleusMass(Z,A) / MeV;
966 else return theG4IonTable->GetNucle <<
967 #endif 513 #endif
968 } else << 514 } else
969 return 0.; << 515 return 0.;
970 } << 516 }
971 517
972 G4double getINCLMass(const G4int A, const << 518 G4double ParticleTable::getINCLMass(const G4int A, const G4int Z) {
973 // assert(A>=0); 519 // assert(A>=0);
974 // For nuclei with Z<0 or Z>A, assume th << 520 // For nuclei with Z<0 or Z>A, assume that the exotic charge state is due to pions
975 // Note that S<0 for lambda << 521 if(Z<0)
976 if(Z<0 && S<0) << 522 return A*neutronMass - Z*getINCLMass(PiMinus);
977 return (A+S)*neutronMass - S*LambdaMas << 523 else if(Z>A)
978 else if(Z>A && S<0) << 524 return A*protonMass + (A-Z)*getINCLMass(PiPlus);
979 return (A+S)*protonMass - S*LambdaMass << 525 else if(A>1)
980 else if(Z<0) << 526 return Z*(protonMass - protonSeparationEnergy) + (A-Z)*(neutronMass - neutronSeparationEnergy);
981 return (A)*neutronMass - Z*getINCLMass << 527 else if(A==1 && Z==0)
982 else if(Z>A) << 528 return getINCLMass(Neutron);
983 return (A)*protonMass + (A-Z)*getINCLM << 529 else if(A==1 && Z==1)
984 else if(A>1 && S<0) << 530 return getINCLMass(Proton);
985 return Z*(protonMass - protonSeparatio << 531 else
986 else if(A>1) << 532 return 0.;
987 return Z*(protonMass - protonSeparatio << 533 }
988 else if(A==1 && Z==0 && S==0) << 534
989 return getINCLMass(Neutron); << 535 G4double ParticleTable::getNuclearRadius(const G4int A, const G4int Z) {
990 else if(A==1 && Z==1 && S==0) << 536 // assert(A>0 && Z>=0 && Z<=A);
991 return getINCLMass(Proton); << 537 if(A >= 19 || (A < 6 && A >= 2)) {
992 else if(A==1 && Z==0 && S==-1) << 538 // For large (Woods-Saxon or Modified Harmonic Oscillator) or small
993 return getINCLMass(Lambda); << 539 // (Gaussian) nuclei, the radius parameter is just the nuclear radius
994 else << 540 return getRadiusParameter(A,Z);
995 return 0.; << 541 } else if(A < clusterTableASize && Z < clusterTableZSize && A >= 6) {
996 } << 542 const G4double thisRMS = positionRMS[Z][A];
997 << 543 if(thisRMS>0.0)
998 G4double getTableQValue(const G4int A1, co << 544 return thisRMS;
999 return getTableMass(A1,Z1,S1) + getTable << 545 else {
1000 } << 546 ERROR("ParticleTable::getRadiusParameter : Radius for nucleus A = " << A << " Z = " << Z << " is ==0.0" << std::endl);
1001 << 547 return 0.0;
1002 G4double getTableQValue(const G4int A1, c <<
1003 return getTableMass(A1,Z1,S1) + getTabl <<
1004 } <<
1005 <<
1006 G4double getTableSpeciesMass(const Partic <<
1007 if(p.theType == Composite) <<
1008 return (*getTableMass)(p.theA, p.theZ <<
1009 else <<
1010 return (*getTableParticleMass)(p.theT <<
1011 } <<
1012 <<
1013 G4int getMassNumber(const ParticleType t) <<
1014 <<
1015 switch(t) { <<
1016 case Proton: <<
1017 case Neutron: <<
1018 case DeltaPlusPlus: <<
1019 case DeltaPlus: <<
1020 case DeltaZero: <<
1021 case DeltaMinus: <<
1022 case SigmaPlus: <<
1023 case SigmaZero: <<
1024 case SigmaMinus: <<
1025 case Lambda: <<
1026 case XiZero: <<
1027 case XiMinus: <<
1028 return 1; <<
1029 break; <<
1030 case antiProton: <<
1031 case antiNeutron: <<
1032 case antiSigmaPlus: <<
1033 case antiSigmaZero: <<
1034 case antiSigmaMinus: <<
1035 case antiLambda: <<
1036 case antiXiZero: <<
1037 case antiXiMinus: <<
1038 return -1; <<
1039 break; <<
1040 case PiPlus: <<
1041 case PiMinus: <<
1042 case PiZero: <<
1043 case KPlus: <<
1044 case KZero: <<
1045 case KZeroBar: <<
1046 case KShort: <<
1047 case KLong: <<
1048 case KMinus: <<
1049 case Eta: <<
1050 case Omega: <<
1051 case EtaPrime: <<
1052 case Photon: <<
1053 return 0; <<
1054 break; <<
1055 default: <<
1056 return 0; <<
1057 break; <<
1058 } <<
1059 } <<
1060 <<
1061 G4int getChargeNumber(const ParticleType <<
1062 switch(t) { <<
1063 case DeltaPlusPlus: <<
1064 return 2; <<
1065 break; <<
1066 case Proton: <<
1067 case DeltaPlus: <<
1068 case PiPlus: <<
1069 case SigmaPlus: <<
1070 case KPlus: <<
1071 case antiSigmaMinus: <<
1072 case antiXiMinus: <<
1073 return 1; <<
1074 break; <<
1075 case Neutron: <<
1076 case DeltaZero: <<
1077 case PiZero: <<
1078 case SigmaZero: <<
1079 case Lambda: <<
1080 case KZero: <<
1081 case KZeroBar: <<
1082 case KShort: <<
1083 case KLong: <<
1084 case Eta: <<
1085 case Omega: <<
1086 case EtaPrime: <<
1087 case Photon: <<
1088 case XiZero: <<
1089 case antiNeutron: <<
1090 case antiLambda: <<
1091 case antiSigmaZero: <<
1092 case antiXiZero: <<
1093 return 0; <<
1094 break; <<
1095 case DeltaMinus: <<
1096 case PiMinus: <<
1097 case SigmaMinus: <<
1098 case KMinus: <<
1099 case antiProton: <<
1100 case XiMinus: <<
1101 case antiSigmaPlus: <<
1102 return -1; <<
1103 break; <<
1104 default: <<
1105 return 0; <<
1106 break; <<
1107 } <<
1108 } <<
1109 <<
1110 G4int getStrangenessNumber(const Particle <<
1111 switch(t) { <<
1112 case DeltaPlusPlus: <<
1113 case DeltaPlus: <<
1114 case DeltaZero: <<
1115 case DeltaMinus: <<
1116 case Proton: <<
1117 case Neutron: <<
1118 case PiPlus: <<
1119 case PiZero: <<
1120 case PiMinus: <<
1121 case Eta: <<
1122 case Omega: <<
1123 case EtaPrime: <<
1124 case Photon: <<
1125 case antiProton: <<
1126 case antiNeutron: <<
1127 return 0; <<
1128 break; <<
1129 case XiMinus: <<
1130 case XiZero: <<
1131 case antiXiMinus: <<
1132 case antiXiZero: <<
1133 return 2; <<
1134 break; <<
1135 case antiLambda: <<
1136 case antiSigmaPlus: <<
1137 case antiSigmaZero: <<
1138 case antiSigmaMinus: <<
1139 return 1; <<
1140 break; <<
1141 case Lambda: <<
1142 case SigmaPlus: <<
1143 case SigmaZero: <<
1144 case SigmaMinus: <<
1145 case KZeroBar: <<
1146 case KMinus: <<
1147 return -1; <<
1148 break; <<
1149 case KPlus: <<
1150 case KZero: <<
1151 return 1; <<
1152 break; <<
1153 case KShort: <<
1154 return 0; <<
1155 break; <<
1156 case KLong: <<
1157 return 0; <<
1158 break; <<
1159 default: <<
1160 return 0; <<
1161 break; <<
1162 } 548 }
1163 } << 549 } else if(A < 19) {
1164 << 550 const G4double theRadiusParameter = getRadiusParameter(A, Z);
1165 G4double getNuclearRadius(const ParticleT << 551 const G4double theDiffusenessParameter = getSurfaceDiffuseness(A, Z);
1166 // assert(A>=0); << 552 // The formula yields the nuclear RMS radius based on the parameters of
1167 if(A > 19 || (A < 6 && A >= 2)) { << 553 // the nuclear-density function
1168 // For large (Woods-Saxon or Modified << 554 return 1.581*theDiffusenessParameter*
1169 // (Gaussian) nuclei, the radius para << 555 (2.+5.*theRadiusParameter)/(2.+3.*theRadiusParameter);
1170 return getRadiusParameter(t,A,Z); << 556 } else {
1171 } else if(A < clusterTableASize && Z>=0 << 557 ERROR("ParticleTable::getNuclearRadius: No radius for nucleus A = " << A << " Z = " << Z << std::endl);
>> 558 return 0.0;
>> 559 }
>> 560 }
>> 561
>> 562 G4double ParticleTable::getRadiusParameter(const G4int A, const G4int Z) {
>> 563 // assert(A>0 && Z>=0 && Z<=A);
>> 564 if(A >= 28) {
>> 565 // phenomenological radius fit
>> 566 return (2.745e-4 * A + 1.063) * std::pow(A, 1.0/3.0);
>> 567 } else if(A < 6 && A >= 2) {
>> 568 if(Z<clusterTableZSize) {
1172 const G4double thisRMS = positionRMS[ 569 const G4double thisRMS = positionRMS[Z][A];
1173 if(thisRMS>0.0) 570 if(thisRMS>0.0)
1174 return thisRMS; 571 return thisRMS;
1175 else { 572 else {
1176 INCL_DEBUG("getNuclearRadius: Radiu << 573 ERROR("ParticleTable::getRadiusParameter : Radius for nucleus A = " << A << " Z = " << Z << " is ==0.0" << std::endl);
1177 << "returning radius for << 574 return 0.0;
1178 return positionRMS[6][12]; <<
1179 } <<
1180 } else if(A <= 19) { <<
1181 const G4double theRadiusParameter = g <<
1182 const G4double theDiffusenessParamete <<
1183 // The formula yields the nuclear RMS <<
1184 // the nuclear-density function <<
1185 return 1.225*theDiffusenessParameter* <<
1186 std::sqrt((2.+5.*theRadiusParameter <<
1187 } else { <<
1188 INCL_ERROR("getNuclearRadius: No radi <<
1189 return 0.0; <<
1190 } <<
1191 } <<
1192 <<
1193 G4double getLargestNuclearRadius(const G4 <<
1194 return Math::max(getNuclearRadius(Proto <<
1195 } <<
1196 <<
1197 G4double getRadiusParameter(const Particl <<
1198 // assert(A>0); <<
1199 if(A > 19) { <<
1200 // radius fit for lambdas <<
1201 if(t==Lambda){ <<
1202 G4double r0 = (1.128+0.439*std::pow( <<
1203 return r0; <<
1204 } <<
1205 // phenomenological radius fit <<
1206 G4double r0 = (2.745e-4 * A + 1.063) <<
1207 // HFB calculations <<
1208 if(getRPCorrelationCoefficient(t)<1.) <<
1209 G4double r0hfb = HFB::getRadiusParam <<
1210 if(r0hfb>0.)r0 = r0hfb; <<
1211 } <<
1212 // <<
1213 if(t==Neutron) <<
1214 r0 += neutronSkin; <<
1215 return r0; <<
1216 } else if(A < 6 && A >= 2) { <<
1217 if(Z<clusterTableZSize && Z>=0) { <<
1218 const G4double thisRMS = positionRM <<
1219 if(thisRMS>0.0) <<
1220 return thisRMS; <<
1221 else { <<
1222 INCL_DEBUG("getRadiusParameter: R <<
1223 << "returning radius f <<
1224 return positionRMS[6][12]; <<
1225 } <<
1226 } else { <<
1227 INCL_DEBUG("getRadiusParameter: Rad <<
1228 << "returning radius for <<
1229 return positionRMS[6][12]; <<
1230 } <<
1231 } else if(A <= 19 && A >= 6) { <<
1232 if(t==Lambda){ <<
1233 G4double r0 = (1.128+0.439*std::pow( <<
1234 return r0; <<
1235 } <<
1236 // HFB calculations <<
1237 if(getRPCorrelationCoefficient(t)<1.) <<
1238 G4double r0hfb = HFB::getSurfaceDiff <<
1239 if(r0hfb>0.)return r0hfb; <<
1240 } <<
1241 return mediumRadius[A-1]; <<
1242 // return 1.581*mediumDiffusenes <<
1243 } else { <<
1244 INCL_ERROR("getRadiusParameter: No ra <<
1245 return 0.0; <<
1246 } <<
1247 } <<
1248 <<
1249 G4double getMaximumNuclearRadius(const Pa <<
1250 const G4double XFOISA = 8.0; <<
1251 if(A > 19) { <<
1252 return getNuclearRadius(t,A,Z) + XFOI <<
1253 } else if(A <= 19 && A >= 6) { <<
1254 return 5.5 + 0.3 * (G4double(A) - 6.0 <<
1255 } else if(A >= 2) { <<
1256 return getNuclearRadius(t, A, Z) + 4. <<
1257 } else { <<
1258 INCL_ERROR("getMaximumNuclearRadius : <<
1259 return 0.0; <<
1260 } <<
1261 } <<
1262 <<
1263 G4double getSurfaceDiffuseness(const Part <<
1264 if(A > 19) { <<
1265 // phenomenological fit <<
1266 G4double a = 1.63e-4 * A + 0.510; <<
1267 // HFB calculations <<
1268 if(getRPCorrelationCoefficient(t)<1.) <<
1269 G4double ahfb = HFB::getSurfaceDiff <<
1270 if(ahfb>0.)a=ahfb; <<
1271 } <<
1272 // <<
1273 if(t==Lambda){ <<
1274 // Like for neutrons <<
1275 G4double ahfb = HFB::getSurfaceDiff <<
1276 if(ahfb>0.)a=ahfb; <<
1277 } <<
1278 if(t==Neutron) <<
1279 a += neutronHalo; <<
1280 return a; <<
1281 } else if(A <= 19 && A >= 6) { <<
1282 // HFB calculations <<
1283 if(getRPCorrelationCoefficient(t)<1.) <<
1284 G4double ahfb = HFB::getRadiusParam <<
1285 if(ahfb>0.)return ahfb; <<
1286 } 575 }
1287 return mediumDiffuseness[A-1]; <<
1288 } else if(A < 6 && A >= 2) { <<
1289 INCL_ERROR("getSurfaceDiffuseness: wa <<
1290 return 0.0; <<
1291 } else { 576 } else {
1292 INCL_ERROR("getSurfaceDiffuseness: No << 577 ERROR("ParticleTable::getRadiusParameter : No radius for nucleus A = " << A << " Z = " << Z << std::endl);
1293 return 0.0; 578 return 0.0;
1294 } 579 }
>> 580 } else if(A < 28 && A >= 6) {
>> 581 return mediumRadius[A-1];
>> 582 // return 1.581*mediumDiffuseness[A-1]*(2.+5.*mediumRadius[A-1])/(2.+3.*mediumRadius[A-1]);
>> 583 } else {
>> 584 ERROR("ParticleTable::getRadiusParameter: No radius for nucleus A = " << A << " Z = " << Z << std::endl);
>> 585 return 0.0;
>> 586 }
>> 587 }
>> 588
>> 589 G4double ParticleTable::getMaximumNuclearRadius(const G4int A, const G4int Z) {
>> 590 const G4double XFOISA = 8.0;
>> 591 if(A >= 19) {
>> 592 return getNuclearRadius(A,Z) + XFOISA * getSurfaceDiffuseness(A,Z);
>> 593 } else if(A < 19 && A >= 6) {
>> 594 return 5.5 + 0.3 * (G4double(A) - 6.0)/12.0;
>> 595 } else if(A >= 2) {
>> 596 return ParticleTable::getNuclearRadius(A, Z) + 4.5;
>> 597 } else {
>> 598 ERROR("ParticleTable::getMaximumNuclearRadius : No maximum radius for nucleus A = " << A << " Z = " << Z << std::endl);
>> 599 return 0.0;
1295 } 600 }
>> 601 }
1296 602
1297 G4double getMomentumRMS(const G4int A, co << 603 #ifdef INCLXX_IN_GEANT4_MODE
1298 // assert(Z>=0 && A>=0 && Z<=A); << 604 G4double ParticleTable::getSurfaceDiffuseness(const G4int A, const G4int /*Z*/ ) {
1299 return getFermiMomentum(A,Z) * Math::sq << 605 #else
1300 } << 606 G4double ParticleTable::getSurfaceDiffuseness(const G4int A, const G4int Z) {
>> 607 #endif // INCLXX_IN_GEANT4_MODE
1301 608
1302 G4double getSeparationEnergyINCL(const Pa << 609 if(A >= 28) {
1303 if(t==Proton) << 610 return 1.63e-4 * A + 0.510;
1304 return theINCLProtonSeparationEnergy; << 611 } else if(A < 28 && A >= 19) {
1305 else if(t==Neutron) << 612 return mediumDiffuseness[A-1];
1306 return theINCLNeutronSeparationEnergy << 613 } else if(A < 19 && A >= 6) {
1307 else if(t==Lambda) << 614 return mediumDiffuseness[A-1];
1308 return theINCLLambdaSeparationEnergy; << 615 } else if(A < 6 && A >= 2) {
1309 else if(t==antiProton) << 616 ERROR("ParticleTable::getSurfaceDiffuseness: was called for A = " << A << " Z = " << Z << std::endl);
1310 return theINCLantiProtonSeparationEne << 617 return 0.0;
1311 else { << 618 } else {
1312 INCL_ERROR("ParticleTable::getSeparat << 619 ERROR("ParticleTable::getSurfaceDiffuseness: No diffuseness for nucleus A = " << A << " Z = " << Z << std::endl);
1313 return 0.0; << 620 return 0.0;
1314 } << 621 }
1315 } << 622 }
>> 623
>> 624 std::string ParticleTable::getElementName(const G4int Z) {
>> 625 if(Z<1) {
>> 626 WARN("getElementName called with Z<1" << std::endl);
>> 627 return elementTable[0];
>> 628 } else if(Z<elementTableSize)
>> 629 return elementTable[Z];
>> 630 else
>> 631 return getIUPACElementName(Z);
>> 632 }
1316 633
1317 G4double getSeparationEnergyReal(const Pa << 634 #ifndef INCLXX_IN_GEANT4_MODE
1318 // Real separation energies for all nuc << 635 void ParticleTable::readRealMasses(std::string const &path) {
1319 if(t==Proton) << 636 // Clear the existing tables, if any
1320 return (*getTableParticleMass)(Proton << 637 massTableMask.clear();
1321 else if(t==Neutron) << 638 massTable.clear();
1322 return (*getTableParticleMass)(Neutro << 639
1323 else if(t==Lambda) << 640 // File name
1324 return (*getTableParticleMass)(Lambda << 641 std::string fileName(path + "/walletlifetime.dat");
1325 else { << 642 DEBUG("Reading real nuclear masses from file " << fileName << std::endl);
1326 INCL_ERROR("ParticleTable::getSeparat << 643
1327 return 0.0; << 644 // Open the file stream
>> 645 std::ifstream massTableIn(fileName.c_str());
>> 646 if(!massTableIn.good()) {
>> 647 FATAL("Cannot open " << fileName << " data file." << std::endl);
>> 648 return;
>> 649 }
>> 650
>> 651 // Read in the data
>> 652 unsigned int Z, A;
>> 653 const G4double amu = 931.494061; // atomic mass unit in MeV/c^2
>> 654 const G4double eMass = 0.5109988; // electron mass in MeV/c^2
>> 655 G4double excess;
>> 656 massTableIn >> A >> Z >> excess;
>> 657 do {
>> 658 // Dynamically determine the table size
>> 659 if(Z>=massTable.size()) {
>> 660 massTable.resize(Z+1);
>> 661 massTableMask.resize(Z+1);
>> 662 }
>> 663 if(A>=massTable[Z].size()) {
>> 664 massTable[Z].resize(A+1, 0.);
>> 665 massTableMask[Z].resize(A+1, false);
1328 } 666 }
1329 } <<
1330 <<
1331 G4double getSeparationEnergyRealForLight( <<
1332 // Real separation energies for light n <<
1333 if(Z<clusterTableZSize && A<clusterTabl <<
1334 return getSeparationEnergyReal(t, A, <<
1335 else <<
1336 return getSeparationEnergyINCL(t, A, <<
1337 } <<
1338 <<
1339 G4double getProtonSeparationEnergy() { re <<
1340 <<
1341 G4double getNeutronSeparationEnergy() { r <<
1342 <<
1343 G4double getLambdaSeparationEnergy() { re <<
1344 <<
1345 void setProtonSeparationEnergy(const G4do <<
1346 <<
1347 void setNeutronSeparationEnergy(const G4d <<
1348 <<
1349 void setLambdaSeparationEnergy(const G4do <<
1350 <<
1351 std::string getElementName(const G4int Z) <<
1352 if(Z<1) { <<
1353 INCL_WARN("getElementName called with <<
1354 return elementTable[0]; <<
1355 } else if(Z<elementTableSize) <<
1356 return elementTable[Z]; <<
1357 else <<
1358 return getIUPACElementName(Z); <<
1359 } <<
1360 <<
1361 std::string getIUPACElementName(const G4i <<
1362 std::stringstream elementStream; <<
1363 elementStream << Z; <<
1364 std::string elementName = elementStream <<
1365 std::transform(elementName.begin(), ele <<
1366 elementName[0] = (char)std::toupper(ele <<
1367 return elementName; <<
1368 } <<
1369 <<
1370 G4int parseElement(std::string pS) { <<
1371 // Normalize the element name <<
1372 std::transform(pS.begin(), pS.end(), pS <<
1373 pS[0] = (char)std::toupper(pS[0]); <<
1374 <<
1375 const std::string *iter = std::find(ele <<
1376 if(iter != elementTable+elementTableSiz <<
1377 return G4int(iter - elementTable); <<
1378 else <<
1379 return ParticleTable::parseIUPACEleme <<
1380 } <<
1381 <<
1382 G4int parseIUPACElement(std::string const <<
1383 // Normalise to lower case <<
1384 std::string elementName(sel); <<
1385 std::transform(elementName.begin(), ele <<
1386 // Return 0 if the element name contain <<
1387 if(elementName.find_first_not_of(elemen <<
1388 return 0; <<
1389 std::transform(elementName.begin(), ele <<
1390 std::stringstream elementStream(element <<
1391 G4int Z; <<
1392 elementStream >> Z; <<
1393 return Z; <<
1394 } <<
1395 <<
1396 IsotopicDistribution const &getNaturalIso <<
1397 return getNaturalIsotopicDistributions( <<
1398 } <<
1399 <<
1400 G4int drawRandomNaturalIsotope(const G4in <<
1401 return getNaturalIsotopicDistributions( <<
1402 } <<
1403 <<
1404 G4double getFermiMomentumConstant(const G <<
1405 return constantFermiMomentum; <<
1406 } <<
1407 <<
1408 G4double getFermiMomentumConstantLight(co <<
1409 // assert(Z>0 && A>0 && Z<=A); <<
1410 if(Z<clusterTableZSize && A<clusterTabl <<
1411 const G4double rms = momentumRMS[Z][A <<
1412 return ((rms>0.) ? rms : momentumRMS[ <<
1413 } else <<
1414 return getFermiMomentumConstant(A,Z); <<
1415 } <<
1416 <<
1417 G4double getFermiMomentumMassDependent(co <<
1418 // assert(A>0); <<
1419 static const G4double alphaParam = 259. <<
1420 static const G4double betaParam = 152. <<
1421 static const G4double gammaParam = 9.51 <<
1422 return alphaParam - betaParam*std::exp( <<
1423 } <<
1424 <<
1425 G4double getRPCorrelationCoefficient(cons <<
1426 // assert(t==Proton || t==Neutron || t==Lambd <<
1427 return rpCorrelationCoefficient[t]; <<
1428 } <<
1429 <<
1430 G4double getNeutronSkin() { return neutro <<
1431 667
1432 G4double getNeutronHalo() { return neutro << 668 massTable.at(Z).at(A) = A*amu + excess - Z*eMass;
>> 669 massTableMask.at(Z).at(A) = true;
1433 670
1434 G4ThreadLocal G4double minDeltaMass = 0.; << 671 massTableIn >> A >> Z >> excess;
1435 G4ThreadLocal G4double minDeltaMass2 = 0. << 672 } while(massTableIn.good());
1436 G4ThreadLocal G4double minDeltaMassRndm = <<
1437 G4ThreadLocal NuclearMassFn getTableMass <<
1438 G4ThreadLocal ParticleMassFn getTablePart <<
1439 G4ThreadLocal SeparationEnergyFn getSepar <<
1440 G4ThreadLocal FermiMomentumFn getFermiMom <<
1441 <<
1442 ParticleType getPionType(const G4int isos <<
1443 // assert(isosp == -2 || isosp == 0 || isosp <<
1444 if (isosp == -2) { <<
1445 return PiMinus; <<
1446 } <<
1447 else if (isosp == 0) { <<
1448 return PiZero; <<
1449 } <<
1450 else { <<
1451 return PiPlus; <<
1452 } <<
1453 } <<
1454 673
1455 ParticleType getNucleonType(const G4int i << 674 }
1456 // assert(isosp == -1 || isosp == 1); << 675 #endif
1457 if (isosp == -1) { <<
1458 return Neutron; <<
1459 } <<
1460 else { <<
1461 return Proton; <<
1462 } <<
1463 } <<
1464 <<
1465 ParticleType getDeltaType(const G4int iso <<
1466 // assert(isosp == -3 || isosp == -1 || isosp <<
1467 if (isosp == -3) { <<
1468 return DeltaMinus; <<
1469 } <<
1470 else if (isosp == -1) { <<
1471 return DeltaZero; <<
1472 } <<
1473 else if (isosp == 1) { <<
1474 return DeltaPlus; <<
1475 } <<
1476 else { <<
1477 return DeltaPlusPlus; <<
1478 } <<
1479 } <<
1480 <<
1481 ParticleType getSigmaType(const G4int iso <<
1482 // assert(isosp == -2 || isosp == 0 || isosp <<
1483 if (isosp == -2) { <<
1484 return SigmaMinus; <<
1485 } <<
1486 else if (isosp == 0) { <<
1487 return SigmaZero; <<
1488 } <<
1489 else { <<
1490 return SigmaPlus; <<
1491 } <<
1492 } <<
1493 <<
1494 ParticleType getXiType(const G4int isosp) <<
1495 // assert(isosp == -1 || isosp == 1); <<
1496 if (isosp == -1) { <<
1497 return XiMinus; <<
1498 } <<
1499 else { <<
1500 return XiZero; <<
1501 } <<
1502 } <<
1503 <<
1504 /*ParticleType getAntiNucleonType(const G4i <<
1505 // assert(isosp == -1); //|| isosp == 1 <<
1506 if (isosp == -1) { <<
1507 return antiProton; <<
1508 } <<
1509 else { <<
1510 return antiNeutron; <<
1511 } <<
1512 }*/ <<
1513 <<
1514 ParticleType getAntiSigmaType(const G4int <<
1515 // assert(isosp == -2 || isosp == 0 || isosp <<
1516 if (isosp == -2) { <<
1517 return antiSigmaPlus; <<
1518 } <<
1519 else if (isosp == 0) { <<
1520 return antiSigmaZero; <<
1521 } <<
1522 else { <<
1523 return antiSigmaMinus; <<
1524 } <<
1525 } <<
1526 <<
1527 ParticleType getAntiXiType(const G4int is <<
1528 // assert(isosp == -1 || isosp == 1); <<
1529 if (isosp == -1) { <<
1530 return antiXiZero; <<
1531 } <<
1532 else { <<
1533 return antiXiMinus; <<
1534 } <<
1535 } <<
1536 <<
1537 ParticleType getKaonType(const G4int isos <<
1538 // assert(isosp == -1 || isosp == 1); <<
1539 if (isosp == -1) { <<
1540 return KZero; <<
1541 } <<
1542 else { <<
1543 return KPlus; <<
1544 } <<
1545 } <<
1546 <<
1547 ParticleType getAntiKaonType(const G4int <<
1548 // assert(isosp == -1 || isosp == 1); <<
1549 if (isosp == -1) { <<
1550 return KMinus; <<
1551 } <<
1552 else { <<
1553 return KZeroBar; <<
1554 } <<
1555 } <<
1556 676
1557 G4double getWidth(const ParticleType pt) << 677 std::string ParticleTable::getIUPACElementName(const G4int Z) {
1558 // assert(pt == PiPlus || pt == PiMinus || pt << 678 std::stringstream elementStream;
1559 if(pt == PiPlus) { << 679 elementStream << Z;
1560 return piPlusWidth; << 680 std::string elementName = elementStream.str();
1561 } else if(pt == PiMinus) { << 681 std::transform(elementName.begin(), elementName.end(), elementName.begin(), ParticleTable::intToIUPAC);
1562 return piMinusWidth; << 682 elementName[0] = std::toupper(elementName.at(0));
1563 } else if(pt == PiZero) { << 683 return elementName;
1564 return piZeroWidth; << 684 }
1565 } else if(pt == Eta) { << 685
1566 return etaWidth; << 686 G4int ParticleTable::parseIUPACElement(std::string const &s) {
1567 } else if(pt == Omega) { << 687 // Normalise to lower case
1568 return omegaWidth; << 688 std::string elementName(s);
1569 } else if(pt == EtaPrime) { << 689 std::transform(elementName.begin(), elementName.end(), elementName.begin(), ::tolower);
1570 return etaPrimeWidth; << 690 // Return 0 if the element name contains anything but IUPAC digits
1571 } else if(pt == SigmaPlus) { << 691 if(elementName.find_first_not_of(elementIUPACDigits)!=std::string::npos)
1572 return SigmaPlusWidth; << 692 return 0;
1573 } else if(pt == SigmaZero) { << 693 std::transform(elementName.begin(), elementName.end(), elementName.begin(), ParticleTable::iupacToInt);
1574 return SigmaZeroWidth; << 694 std::stringstream elementStream(elementName);
1575 } else if(pt == SigmaMinus) { << 695 G4int Z;
1576 return SigmaMinusWidth; << 696 elementStream >> Z;
1577 } else if(pt == KPlus) { << 697 return Z;
1578 return KPlusWidth; << 698 }
1579 } else if(pt == KMinus) { <<
1580 return KMinusWidth; <<
1581 } else if(pt == KShort) { <<
1582 return KShortWidth; <<
1583 } else if(pt == KLong) { <<
1584 return KLongWidth; <<
1585 } else if(pt == Lambda) { <<
1586 return LambdaWidth; <<
1587 } else if(pt == XiMinus) { <<
1588 return XiMinusWidth; <<
1589 } else if(pt == XiZero) { <<
1590 return XiZeroWidth; <<
1591 } else if(pt == antiSigmaPlus) { <<
1592 return antiSigmaPlusWidth; <<
1593 } else if(pt == antiSigmaZero) { <<
1594 return antiSigmaZeroWidth; <<
1595 } else if(pt == antiSigmaMinus) { <<
1596 return antiSigmaMinusWidth; <<
1597 } else if(pt == antiLambda) { <<
1598 return antiLambdaWidth; <<
1599 } else if(pt == antiXiMinus) { <<
1600 return antiXiMinusWidth; <<
1601 } else if(pt == antiXiZero) { <<
1602 return antiXiZeroWidth; <<
1603 } else { <<
1604 INCL_ERROR("getWidth : Unknown <<
1605 return 0.0; <<
1606 } <<
1607 } <<
1608 <<
1609 } // namespace ParticleTable <<
1610 } // namespace G4INCL <<
1611 699
>> 700 }
1612 701