Geant4 11.1.1
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4INCLParticleTable.cc
Go to the documentation of this file.
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26// INCL++ intra-nuclear cascade model
27// Alain Boudard, CEA-Saclay, France
28// Joseph Cugnon, University of Liege, Belgium
29// Jean-Christophe David, CEA-Saclay, France
30// Pekka Kaitaniemi, CEA-Saclay, France, and Helsinki Institute of Physics, Finland
31// Sylvie Leray, CEA-Saclay, France
32// Davide Mancusi, CEA-Saclay, France
33//
34#define INCLXX_IN_GEANT4_MODE 1
35
36#include "globals.hh"
37
38#include "G4INCLParticleTable.hh"
39#include "G4INCLNuclearMassTable.hh"
40#include <algorithm>
41// #include <cassert>
42#include <cmath>
43#include <cctype>
44#include <sstream>
45#ifdef INCLXX_IN_GEANT4_MODE
46#include "G4SystemOfUnits.hh"
47#endif
48
49#ifdef INCLXX_IN_GEANT4_MODE
50#include "G4PhysicalConstants.hh"
51#include "G4SystemOfUnits.hh"
52#endif
53
54namespace G4INCL {
55
56 namespace ParticleTable {
57
58 namespace {
59
60 /// \brief Static instance of the NaturalIsotopicAbundances class
61 const NaturalIsotopicDistributions *theNaturalIsotopicDistributions = NULL;
62
63 const G4double theINCLNucleonMass = 938.2796;
64 const G4double theINCLPionMass = 138.0;
65 const G4double theINCLLambdaMass = 1115.683;
66// const G4double theINCLSigmaMass = 1197.45;
67// const G4double theINCLKaonMass = 497.614;
68 const G4double theINCLEtaMass = 547.862;
69 const G4double theINCLOmegaMass = 782.65;
70 const G4double theINCLEtaPrimeMass = 957.78;
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.0;
78 G4ThreadLocal G4double SigmaZeroMass = 0.0;
79 G4ThreadLocal G4double SigmaMinusMass = 0.0;
80 G4ThreadLocal G4double LambdaMass = 0.0;
81 G4ThreadLocal G4double KPlusMass = 0.0;
82 G4ThreadLocal G4double KZeroMass = 0.0;
83 G4ThreadLocal G4double KZeroBarMass = 0.0;
84 G4ThreadLocal G4double KShortMass = 0.0;
85 G4ThreadLocal G4double KLongMass = 0.0;
86 G4ThreadLocal G4double KMinusMass = 0.0;
87 G4ThreadLocal G4double etaMass = 0.0;
88 G4ThreadLocal G4double omegaMass = 0.0;
89 G4ThreadLocal G4double etaPrimeMass = 0.0;
90 G4ThreadLocal G4double photonMass = 0.0;
91
92 // Hard-coded values of the real particle masses (MeV/c^2)
93 G4ThreadLocal G4double theRealProtonMass = 938.27203;
94 G4ThreadLocal G4double theRealNeutronMass = 939.56536;
95 G4ThreadLocal G4double theRealChargedPiMass = 139.57018;
96 G4ThreadLocal G4double theRealPiZeroMass = 134.9766;
97 G4ThreadLocal G4double theRealLambdaMass = 1115.683;
98 G4ThreadLocal G4double theRealSigmaPlusMass = 1189.37;
99 G4ThreadLocal G4double theRealSigmaZeroMass = 1192.64;
100 G4ThreadLocal G4double theRealSigmaMinusMass = 1197.45;
101 G4ThreadLocal G4double theRealChargedKaonMass = 493.677;
102 G4ThreadLocal G4double theRealNeutralKaonMass = 497.614;
103 G4ThreadLocal G4double theRealEtaMass = 547.862;
104 G4ThreadLocal G4double theRealOmegaMass = 782.65;
105 G4ThreadLocal G4double theRealEtaPrimeMass = 957.78;
106 G4ThreadLocal G4double theRealPhotonMass = 0.0;
107
108 // Width (second)
109 const G4double theChargedPiWidth = 2.6033e-08;
110 const G4double thePiZeroWidth = 8.52e-17;
111 const G4double theEtaWidth = 5.025e-19; // 1.31 keV
112 const G4double theOmegaWidth = 7.7528e-23; // 8.49 MeV
113 const G4double theEtaPrimeWidth = 3.3243e-21; // 0.198 MeV
114 const G4double theChargedKaonWidth = 1.238e-08;
115 const G4double theKShortWidth = 8.954e-11;
116 const G4double theKLongWidth = 5.116e-08;
117 const G4double theLambdaWidth = 2.632e-10;
118 const G4double theSigmaPlusWidth = 8.018e-11;
119 const G4double theSigmaZeroWidth = 7.4e-20;
120 const G4double theSigmaMinusWidth = 1.479e-10;
121 G4ThreadLocal G4double piPlusWidth = 0.0;
122 G4ThreadLocal G4double piMinusWidth = 0.0;
123 G4ThreadLocal G4double piZeroWidth = 0.0;
124 G4ThreadLocal G4double etaWidth = 0.0;
125 G4ThreadLocal G4double omegaWidth = 0.0;
126 G4ThreadLocal G4double etaPrimeWidth = 0.0;
127 G4ThreadLocal G4double LambdaWidth = 0.0;
128 G4ThreadLocal G4double SigmaPlusWidth = 0.0;
129 G4ThreadLocal G4double SigmaZeroWidth = 0.0;
130 G4ThreadLocal G4double SigmaMinusWidth = 0.0;
131 G4ThreadLocal G4double KPlusWidth = 0.0;
132 G4ThreadLocal G4double KMinusWidth = 0.0;
133 G4ThreadLocal G4double KShortWidth = 0.0;
134 G4ThreadLocal G4double KLongWidth = 0.0;
135
136
137 const G4int mediumNucleiTableSize = 30;
138
139 const G4double mediumDiffuseness[mediumNucleiTableSize] =
140 {0.0,0.0,0.0,0.0,0.0,1.78,1.77,1.77,1.69,1.71,
141 1.69,1.72,1.635,1.730,1.81,1.833,1.798,
142 1.93,0.567,0.571, 0.560,0.549,0.550,0.551,
143 0.580,0.575,0.569,0.537,0.0,0.0};
144 const G4double mediumRadius[mediumNucleiTableSize] =
145 {0.0,0.0,0.0,0.0,0.0,0.334,0.327,0.479,0.631,0.838,
146 0.811,0.84,1.403,1.335,1.25,1.544,1.498,1.57,
147 2.58,2.77, 2.775,2.78,2.88,2.98,3.22,3.03,2.84,
148 3.14,0.0,0.0};
149
150 const G4double positionRMS[clusterTableZSize][clusterTableASize] = {
151 /* A= 0 1 2 3 4 5 6 7 8 9 10 11 12 */
152 /* Z=0 */ {-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0},
153 /* Z=1 */ {-1.0, -1.0, 2.10, 1.80, 1.70, 1.83, 2.60, 2.50, -1.0, -1.0, -1.0, -1.0, -1.0},
154 /* Z=2 */ {-1.0, -1.0, -1.0, 1.80, 1.68, 1.70, 2.60, 2.50, 2.50, 2.50, 2.50, -1.0, -1.0},
155 /* Z=3 */ {-1.0, -1.0, -1.0, -1.0, 1.70, 1.83, 2.56, 2.40, 2.50, 2.50, 2.50, 2.50, 2.50},
156 /* Z=4 */ {-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 2.60, 2.50, 2.50, 2.51, 2.50, 2.50, 2.50},
157 /* Z=5 */ {-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 2.50, 2.50, 2.50, 2.50, 2.45, 2.40, 2.50},
158 /* Z=6 */ {-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 2.50, 2.50, 2.50, 2.50, 2.47},
159 /* Z=7 */ {-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 2.50, 2.50, 2.50},
160 /* Z=8 */ {-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 2.50}
161 };
162
163 const G4double momentumRMS[clusterTableZSize][clusterTableASize] = {
164 /* A= 0 1 2 3 4 5 6 7 8 9 10 11 12 */
165 /* Z=0 */ {-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0},
166 /* Z=1 */ {-1.0, -1.0, 77.0, 110., 153., 100., 100., 100., -1.0, -1.0, -1.0, -1.0, -1.0},
167 /* Z=2 */ {-1.0, -1.0, -1.0, 110., 153., 100., 100., 100., 100., 100., 100., -1.0, -1.0},
168 /* Z=3 */ {-1.0, -1.0, -1.0, -1.0, 153., 100., 100., 100., 100., 100., 100., 100., 100.},
169 /* Z=4 */ {-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 100., 100., 100., 100., 100., 100., 100.},
170 /* Z=5 */ {-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 100., 100., 100., 100., 100., 100., 100.},
171 /* Z=6 */ {-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 100., 100., 100., 100., 100.},
172 /* Z=7 */ {-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 100., 100., 100.},
173 /* Z=8 */ {-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 100.}
174 };
175
176 const G4int elementTableSize = 113; // up to Cn
177
178 /// \brief Table of chemical element names
179 const std::string elementTable[elementTableSize] = {
180 "",
181 "H",
182 "He",
183 "Li",
184 "Be",
185 "B",
186 "C",
187 "N",
188 "O",
189 "F",
190 "Ne",
191 "Na",
192 "Mg",
193 "Al",
194 "Si",
195 "P",
196 "S",
197 "Cl",
198 "Ar",
199 "K",
200 "Ca",
201 "Sc",
202 "Ti",
203 "V",
204 "Cr",
205 "Mn",
206 "Fe",
207 "Co",
208 "Ni",
209 "Cu",
210 "Zn",
211 "Ga",
212 "Ge",
213 "As",
214 "Se",
215 "Br",
216 "Kr",
217 "Rb",
218 "Sr",
219 "Y",
220 "Zr",
221 "Nb",
222 "Mo",
223 "Tc",
224 "Ru",
225 "Rh",
226 "Pd",
227 "Ag",
228 "Cd",
229 "In",
230 "Sn",
231 "Sb",
232 "Te",
233 "I",
234 "Xe",
235 "Cs",
236 "Ba",
237 "La",
238 "Ce",
239 "Pr",
240 "Nd",
241 "Pm",
242 "Sm",
243 "Eu",
244 "Gd",
245 "Tb",
246 "Dy",
247 "Ho",
248 "Er",
249 "Tm",
250 "Yb",
251 "Lu",
252 "Hf",
253 "Ta",
254 "W",
255 "Re",
256 "Os",
257 "Ir",
258 "Pt",
259 "Au",
260 "Hg",
261 "Tl",
262 "Pb",
263 "Bi",
264 "Po",
265 "At",
266 "Rn",
267 "Fr",
268 "Ra",
269 "Ac",
270 "Th",
271 "Pa",
272 "U",
273 "Np",
274 "Pu",
275 "Am",
276 "Cm",
277 "Bk",
278 "Cf",
279 "Es",
280 "Fm",
281 "Md",
282 "No",
283 "Lr",
284 "Rf",
285 "Db",
286 "Sg",
287 "Bh",
288 "Hs",
289 "Mt",
290 "Ds",
291 "Rg",
292 "Cn"
293 };
294
295 /// \brief Digit names to compose IUPAC element names
296 const std::string elementIUPACDigits = "nubtqphsoe";
297
298#define INCL_DEFAULT_SEPARATION_ENERGY 6.83
299 const G4double theINCLProtonSeparationEnergy = INCL_DEFAULT_SEPARATION_ENERGY;
300 const G4double theINCLNeutronSeparationEnergy = INCL_DEFAULT_SEPARATION_ENERGY;
301 const G4double theINCLLambdaSeparationEnergy = INCL_DEFAULT_SEPARATION_ENERGY;
302 G4ThreadLocal G4double protonSeparationEnergy = INCL_DEFAULT_SEPARATION_ENERGY;
303 G4ThreadLocal G4double neutronSeparationEnergy = INCL_DEFAULT_SEPARATION_ENERGY;
304 G4ThreadLocal G4double lambdaSeparationEnergy = INCL_DEFAULT_SEPARATION_ENERGY;
305#undef INCL_DEFAULT_SEPARATION_ENERGY
306
307 G4ThreadLocal G4double rpCorrelationCoefficient[UnknownParticle];
308
309 G4ThreadLocal G4double neutronSkin = 0.0;
310 G4ThreadLocal G4double neutronHalo = 0.0;
311
312#ifdef INCLXX_IN_GEANT4_MODE
313 G4ThreadLocal G4IonTable *theG4IonTable;
314#endif
315
316 /// \brief Default value for constant Fermi momentum
317 G4ThreadLocal G4double constantFermiMomentum = 0.0;
318
319 /// \brief Transform a IUPAC char to an char representing an integer digit
320 char iupacToInt(char c) {
321 return (char)(((G4int)'0')+elementIUPACDigits.find(c));
322 }
323
324 /// \brief Transform an integer digit (represented by a char) to a IUPAC char
325 char intToIUPAC(char n) { return elementIUPACDigits.at(n); }
326
327 /// \brief Get the singleton instance of the natural isotopic distributions
328 const NaturalIsotopicDistributions *getNaturalIsotopicDistributions() {
329 if(!theNaturalIsotopicDistributions)
330 theNaturalIsotopicDistributions = new NaturalIsotopicDistributions;
331 return theNaturalIsotopicDistributions;
332 }
333
334 } // namespace
335
336 void initialize(Config const * const theConfig /*=0*/) {
337 protonMass = theINCLNucleonMass;
338 neutronMass = theINCLNucleonMass;
339 piPlusMass = theINCLPionMass;
340 piMinusMass = theINCLPionMass;
341 piZeroMass = theINCLPionMass;
342 /*
343 SigmaPlusMass = theINCLSigmaMass;
344 SigmaMinusMass = theINCLSigmaMass;
345 SigmaZeroMass = theINCLSigmaMass;
346 LambdaMass = theINCLLambdaMass;
347 KPlusMass = theINCLKaonMass;
348 KZeroMass = theINCLKaonMass;
349 KZeroBarMass = theINCLKaonMass;
350 KShortMass = theINCLKaonMass;
351 KLongMass = theINCLKaonMass;
352 KMinusMass = theINCLKaonMass;
353 */
354 SigmaPlusMass = theRealSigmaPlusMass;
355 SigmaMinusMass = theRealSigmaMinusMass;
356 SigmaZeroMass = theRealSigmaZeroMass;
357 LambdaMass = theINCLLambdaMass;
358 KPlusMass = theRealChargedKaonMass;
359 KZeroMass = theRealNeutralKaonMass;
360 KZeroBarMass = theRealNeutralKaonMass;
361 KShortMass = theRealNeutralKaonMass;
362 KLongMass = theRealNeutralKaonMass;
363 KMinusMass = theRealChargedKaonMass;
364
365 etaMass = theINCLEtaMass;
366 omegaMass = theINCLOmegaMass;
367 etaPrimeMass = theINCLEtaPrimeMass;
368 photonMass = theINCLPhotonMass;
369 if(theConfig && theConfig->getUseRealMasses()) {
370 getTableMass = getRealMass;
371 getTableParticleMass = getRealMass;
372 } else {
373 getTableMass = getINCLMass;
374 getTableParticleMass = getINCLMass;
375 }
376
377#ifndef INCLXX_IN_GEANT4_MODE
378 std::string dataFilePath;
379 if(theConfig)
380 dataFilePath = theConfig->getINCLXXDataFilePath();
381 NuclearMassTable::initialize(dataFilePath, getRealMass(Proton), getRealMass(Neutron));
382#endif
383
384#ifdef INCLXX_IN_GEANT4_MODE
385 G4ParticleTable *theG4ParticleTable = G4ParticleTable::GetParticleTable();
386 theG4IonTable = theG4ParticleTable->GetIonTable();
387 theRealProtonMass = theG4ParticleTable->FindParticle("proton")->GetPDGMass() / MeV;
388 theRealNeutronMass = theG4ParticleTable->FindParticle("neutron")->GetPDGMass() / MeV;
389 theRealChargedPiMass = theG4ParticleTable->FindParticle("pi+")->GetPDGMass() / MeV;
390 theRealPiZeroMass = theG4ParticleTable->FindParticle("pi0")->GetPDGMass() / MeV;
391 theRealEtaMass = theG4ParticleTable->FindParticle("eta")->GetPDGMass() / MeV;
392 theRealOmegaMass = theG4ParticleTable->FindParticle("omega")->GetPDGMass() / MeV;
393 theRealEtaPrimeMass = theG4ParticleTable->FindParticle("eta_prime")->GetPDGMass() / MeV;
394 theRealPhotonMass = theG4ParticleTable->FindParticle("gamma")->GetPDGMass() / MeV;
395 theRealSigmaPlusMass = theG4ParticleTable->FindParticle("sigma+")->GetPDGMass() / MeV;
396 theRealSigmaZeroMass = theG4ParticleTable->FindParticle("sigma0")->GetPDGMass() / MeV;
397 theRealSigmaMinusMass = theG4ParticleTable->FindParticle("sigma-")->GetPDGMass() / MeV;
398 theRealLambdaMass = theG4ParticleTable->FindParticle("lambda")->GetPDGMass() / MeV;
399 theRealChargedKaonMass = theG4ParticleTable->FindParticle("kaon+")->GetPDGMass() / MeV;
400 theRealNeutralKaonMass = theG4ParticleTable->FindParticle("kaon0")->GetPDGMass() / MeV;
401#endif
402
403 minDeltaMass = theRealNeutronMass + theRealChargedPiMass + 0.5;
404 minDeltaMass2 = minDeltaMass*minDeltaMass;
405 minDeltaMassRndm = std::atan((minDeltaMass-effectiveDeltaMass)*2./effectiveDeltaWidth);
406
407 piPlusWidth = theChargedPiWidth;
408 piMinusWidth = theChargedPiWidth;
409 piZeroWidth = thePiZeroWidth;
410 etaWidth = theEtaWidth;
411 omegaWidth = theOmegaWidth;
412 etaPrimeWidth = theEtaPrimeWidth;
413
414 SigmaMinusWidth = theSigmaMinusWidth;
415 SigmaPlusWidth = theSigmaPlusWidth;
416 SigmaZeroWidth = theSigmaZeroWidth;
417 LambdaWidth = theLambdaWidth;
418 KPlusWidth = theChargedKaonWidth;
419 KMinusWidth = theChargedKaonWidth;
420 KShortWidth = theKShortWidth;
421 KLongWidth = theKLongWidth;
422
423 // Initialise HFB tables
424#ifdef INCLXX_IN_GEANT4_MODE
425 HFB::initialize();
426#else
427 HFB::initialize(dataFilePath);
428#endif
429
430 // Initialise the separation-energy function
431 if(!theConfig || theConfig->getSeparationEnergyType()==INCLSeparationEnergy)
432 getSeparationEnergy = getSeparationEnergyINCL;
433 else if(theConfig->getSeparationEnergyType()==RealSeparationEnergy)
434 getSeparationEnergy = getSeparationEnergyReal;
435 else if(theConfig->getSeparationEnergyType()==RealForLightSeparationEnergy)
436 getSeparationEnergy = getSeparationEnergyRealForLight;
437 else {
438 INCL_FATAL("Unrecognized separation-energy type in ParticleTable initialization: " << theConfig->getSeparationEnergyType() << '\n');
439 return;
440 }
441
442 // Initialise the Fermi-momentum function
443 if(!theConfig || theConfig->getFermiMomentumType()==ConstantFermiMomentum) {
444 getFermiMomentum = ParticleTable::getFermiMomentumConstant;
445 if(theConfig) {
446 const G4double aFermiMomentum = theConfig->getFermiMomentum();
447 if(aFermiMomentum>0.)
448 constantFermiMomentum = aFermiMomentum;
449 else
450 constantFermiMomentum = PhysicalConstants::Pf;
451 } else {
452 constantFermiMomentum = PhysicalConstants::Pf;
453 }
454 } else if(theConfig->getFermiMomentumType()==ConstantLightFermiMomentum)
455 getFermiMomentum = ParticleTable::getFermiMomentumConstantLight;
456 else if(theConfig->getFermiMomentumType()==MassDependentFermiMomentum)
457 getFermiMomentum = ParticleTable::getFermiMomentumMassDependent;
458 else {
459 INCL_FATAL("Unrecognized Fermi-momentum type in ParticleTable initialization: " << theConfig->getFermiMomentumType() << '\n');
460 return;
461 }
462
463 // Initialise the r-p correlation coefficients
464 std::fill(rpCorrelationCoefficient, rpCorrelationCoefficient + UnknownParticle, 1.);
465 if(theConfig) {
466 rpCorrelationCoefficient[Proton] = theConfig->getRPCorrelationCoefficient(Proton);
467 rpCorrelationCoefficient[Neutron] = theConfig->getRPCorrelationCoefficient(Neutron);
468 }
469
470 // Initialise the neutron-skin parameters
471 if(theConfig) {
472 neutronSkin = theConfig->getNeutronSkin();
473 neutronHalo = theConfig->getNeutronHalo();
474 }
475
476 }
477
478 G4int getIsospin(const ParticleType t) {
479 // Actually this is the 3rd component of isospin (I_z) multiplied by 2!
480 if(t == Proton) {
481 return 1;
482 } else if(t == Neutron) {
483 return -1;
484 } else if(t == PiPlus) {
485 return 2;
486 } else if(t == PiMinus) {
487 return -2;
488 } else if(t == PiZero) {
489 return 0;
490 } else if(t == DeltaPlusPlus) {
491 return 3;
492 } else if(t == DeltaPlus) {
493 return 1;
494 } else if(t == DeltaZero) {
495 return -1;
496 } else if(t == DeltaMinus) {
497 return -3;
498 } else if(t == Lambda) {
499 return 0;
500 } else if(t == SigmaPlus) {
501 return 2;
502 } else if(t == SigmaZero) {
503 return 0;
504 } else if(t == SigmaMinus) {
505 return -2;
506 } else if(t == KPlus) {
507 return 1;
508 } else if(t == KZero) {
509 return -1;
510 } else if(t == KZeroBar) {
511 return 1;
512 } else if(t == KShort) {
513 return 0;
514 } else if(t == KLong) {
515 return 0;
516 } else if(t == KMinus) {
517 return -1;
518 } else if(t == Eta) {
519 return 0;
520 } else if(t == Omega) {
521 return 0;
522 } else if(t == EtaPrime) {
523 return 0;
524 } else if(t == Photon) {
525 return 0;
526 }
527 INCL_ERROR("Requested isospin of an unknown particle!");
528 return -10; // Unknown
529 }
530
531 std::string getShortName(const ParticleSpecies &sp) {
532 if(sp.theType==Composite && sp.theS == 0)
533 return getShortName(sp.theA,sp.theZ);
534 else if(sp.theType==Composite)
535 return getName(sp.theA,sp.theZ,sp.theS);
536 else
537 return getShortName(sp.theType);
538 }
539
540 std::string getName(const ParticleSpecies &sp) {
541 if(sp.theType==Composite && sp.theS == 0)
542 return getName(sp.theA,sp.theZ);
543 else if(sp.theType==Composite)
544 return getName(sp.theA,sp.theZ,sp.theS);
545 else
546 return getName(sp.theType);
547 }
548
549 std::string getName(const G4int A, const G4int Z) {
550 std::stringstream stream;
551 stream << getElementName(Z) << "-" << A;
552 return stream.str();
553 }
554
555 std::string getName(const G4int A, const G4int Z, const G4int S) {
556 std::stringstream stream;
557 if(S >= 0) // S < 0 for hypernuclei
558 return getName(A, Z);
559 else if(S == -1)
560 stream << getElementName(Z) << "-" << A << "_" << "Lambda";
561 else
562 stream << getElementName(Z) << "-" << A << "_" << S << "-Lambda";
563 return stream.str();
564 }
565
566 std::string getShortName(const G4int A, const G4int Z) {
567 std::stringstream stream;
568 stream << getElementName(Z);
569 if(A>0)
570 stream << A;
571 return stream.str();
572 }
573
574 std::string getName(const ParticleType p) {
575 if(p == G4INCL::Proton) {
576 return std::string("proton");
577 } else if(p == G4INCL::Neutron) {
578 return std::string("neutron");
579 } else if(p == G4INCL::DeltaPlusPlus) {
580 return std::string("delta++");
581 } else if(p == G4INCL::DeltaPlus) {
582 return std::string("delta+");
583 } else if(p == G4INCL::DeltaZero) {
584 return std::string("delta0");
585 } else if(p == G4INCL::DeltaMinus) {
586 return std::string("delta-");
587 } else if(p == G4INCL::PiPlus) {
588 return std::string("pi+");
589 } else if(p == G4INCL::PiZero) {
590 return std::string("pi0");
591 } else if(p == G4INCL::PiMinus) {
592 return std::string("pi-");
593 } else if(p == G4INCL::Lambda) {
594 return std::string("lambda");
595 } else if(p == G4INCL::SigmaPlus) {
596 return std::string("sigma+");
597 } else if(p == G4INCL::SigmaZero) {
598 return std::string("sigma0");
599 } else if(p == G4INCL::SigmaMinus) {
600 return std::string("sigma-");
601 } else if(p == G4INCL::KPlus) {
602 return std::string("kaon+");
603 } else if(p == G4INCL::KZero) {
604 return std::string("kaon0");
605 } else if(p == G4INCL::KZeroBar) {
606 return std::string("kaon0bar");
607 } else if(p == G4INCL::KMinus) {
608 return std::string("kaon-");
609 } else if(p == G4INCL::KShort) {
610 return std::string("kaonshort");
611 } else if(p == G4INCL::KLong) {
612 return std::string("kaonlong");
613 } else if(p == G4INCL::Composite) {
614 return std::string("composite");
615 } else if(p == G4INCL::Eta) {
616 return std::string("eta");
617 } else if(p == G4INCL::Omega) {
618 return std::string("omega");
619 } else if(p == G4INCL::EtaPrime) {
620 return std::string("etaprime");
621 } else if(p == G4INCL::Photon) {
622 return std::string("photon");
623 }
624 return std::string("unknown");
625 }
626
627 std::string getShortName(const ParticleType p) {
628 if(p == G4INCL::Proton) {
629 return std::string("p");
630 } else if(p == G4INCL::Neutron) {
631 return std::string("n");
632 } else if(p == G4INCL::DeltaPlusPlus) {
633 return std::string("d++");
634 } else if(p == G4INCL::DeltaPlus) {
635 return std::string("d+");
636 } else if(p == G4INCL::DeltaZero) {
637 return std::string("d0");
638 } else if(p == G4INCL::DeltaMinus) {
639 return std::string("d-");
640 } else if(p == G4INCL::PiPlus) {
641 return std::string("pi+");
642 } else if(p == G4INCL::PiZero) {
643 return std::string("pi0");
644 } else if(p == G4INCL::PiMinus) {
645 return std::string("pi-");
646 } else if(p == G4INCL::Lambda) {
647 return std::string("l");
648 } else if(p == G4INCL::SigmaPlus) {
649 return std::string("s+");
650 } else if(p == G4INCL::SigmaZero) {
651 return std::string("s0");
652 } else if(p == G4INCL::SigmaMinus) {
653 return std::string("s-");
654 } else if(p == G4INCL::KPlus) {
655 return std::string("k+");
656 } else if(p == G4INCL::KZero) {
657 return std::string("k0");
658 } else if(p == G4INCL::KZeroBar) {
659 return std::string("k0b");
660 } else if(p == G4INCL::KMinus) {
661 return std::string("k-");
662 } else if(p == G4INCL::KShort) {
663 return std::string("ks");
664 } else if(p == G4INCL::KLong) {
665 return std::string("kl");
666 } else if(p == G4INCL::Composite) {
667 return std::string("comp");
668 } else if(p == G4INCL::Eta) {
669 return std::string("eta");
670 } else if(p == G4INCL::Omega) {
671 return std::string("omega");
672 } else if(p == G4INCL::EtaPrime) {
673 return std::string("etap");
674 } else if(p == G4INCL::Photon) {
675 return std::string("photon");
676 }
677 return std::string("unknown");
678 }
679
680 G4double getINCLMass(const ParticleType pt) {
681 if(pt == Proton) {
682 return protonMass;
683 } else if(pt == Neutron) {
684 return neutronMass;
685 } else if(pt == PiPlus) {
686 return piPlusMass;
687 } else if(pt == PiMinus) {
688 return piMinusMass;
689 } else if(pt == PiZero) {
690 return piZeroMass;
691 } else if(pt == SigmaPlus) {
692 return SigmaPlusMass;
693 } else if(pt == SigmaMinus) {
694 return SigmaMinusMass;
695 } else if(pt == SigmaZero) {
696 return SigmaZeroMass;
697 } else if(pt == Lambda) {
698 return LambdaMass;
699 } else if(pt == KPlus) {
700 return KPlusMass;
701 } else if(pt == KZero) {
702 return KZeroMass;
703 } else if(pt == KZeroBar) {
704 return KZeroBarMass;
705 } else if(pt == KMinus) {
706 return KMinusMass;
707 } else if(pt == KShort) {
708 return KShortMass;
709 } else if(pt == KLong) {
710 return KLongMass;
711 } else if(pt == Eta) {
712 return etaMass;
713 } else if(pt == Omega) {
714 return omegaMass;
715 } else if(pt == EtaPrime) {
716 return etaPrimeMass;
717 } else if(pt == Photon) {
718 return photonMass;
719 } else {
720 INCL_ERROR("getMass : Unknown particle type." << '\n');
721 return 0.0;
722 }
723 }
724
725 G4double getRealMass(const ParticleType t) {
726 switch(t) {
727 case Proton:
728 return theRealProtonMass;
729 break;
730 case Neutron:
731 return theRealNeutronMass;
732 break;
733 case PiPlus:
734 case PiMinus:
735 return theRealChargedPiMass;
736 break;
737 case PiZero:
738 return theRealPiZeroMass;
739 break;
740 case SigmaPlus:
741 return theRealSigmaPlusMass;
742 break;
743 case SigmaZero:
744 return theRealSigmaZeroMass;
745 break;
746 case SigmaMinus:
747 return theRealSigmaMinusMass;
748 break;
749 case Lambda:
750 return theRealLambdaMass;
751 break;
752 case KPlus:
753 case KMinus:
754 return theRealChargedKaonMass;
755 break;
756 case KZero:
757 case KZeroBar:
758 case KShort:
759 case KLong:
760 return theRealNeutralKaonMass;
761 break;
762 case Eta:
763 return theRealEtaMass;
764 break;
765 case Omega:
766 return theRealOmegaMass;
767 break;
768 case EtaPrime:
769 return theRealEtaPrimeMass;
770 break;
771 case Photon:
772 return theRealPhotonMass;
773 break;
774 default:
775 INCL_ERROR("Particle::getRealMass : Unknown particle type." << '\n');
776 return 0.0;
777 break;
778 }
779 }
780
781 G4double getRealMass(const G4int A, const G4int Z, const G4int S) {
782// assert(A>=0);
783 // For nuclei with Z<0 or Z>A, assume that the exotic charge state is due to pions
784 if(Z<0 && S<0)
785 return (A+S)*theRealNeutronMass - S*LambdaMass - Z*getRealMass(PiMinus);
786 else if(Z>A && S<0)
787 return (A+S)*theRealProtonMass - S*LambdaMass + (A+S-Z)*getRealMass(PiPlus);
788 if(Z<0)
789 return (A)*theRealNeutronMass - Z*getRealMass(PiMinus);
790 else if(Z>A)
791 return (A)*theRealProtonMass + (A-Z)*getRealMass(PiPlus);
792 else if(Z==0 && S==0)
793 return A*theRealNeutronMass;
794 else if(A==Z)
795 return A*theRealProtonMass;
796 else if(Z==0 && S<0)
797 return (A+S)*theRealNeutronMass-S*LambdaMass;
798 else if(A>1) {
799#ifndef INCLXX_IN_GEANT4_MODE
800 return ::G4INCL::NuclearMassTable::getMass(A,Z,S);
801#else
802 if(S<0) return theG4IonTable->GetNucleusMass(Z,A,std::abs(S)) / MeV;
803 else return theG4IonTable->GetNucleusMass(Z,A) / MeV;
804#endif
805 } else
806 return 0.;
807 }
808
809 G4double getINCLMass(const G4int A, const G4int Z, const G4int S) {
810// assert(A>=0);
811 // For nuclei with Z<0 or Z>A, assume that the exotic charge state is due to pions
812 // Note that S<0 for lambda
813 if(Z<0 && S<0)
814 return (A+S)*neutronMass - S*LambdaMass - Z*getINCLMass(PiMinus);
815 else if(Z>A && S<0)
816 return (A+S)*protonMass - S*LambdaMass + (A+S-Z)*getINCLMass(PiPlus);
817 else if(Z<0)
818 return (A)*neutronMass - Z*getINCLMass(PiMinus);
819 else if(Z>A)
820 return (A)*protonMass + (A-Z)*getINCLMass(PiPlus);
821 else if(A>1 && S<0)
822 return Z*(protonMass - protonSeparationEnergy) + (A+S-Z)*(neutronMass - neutronSeparationEnergy) + std::abs(S)*(LambdaMass - lambdaSeparationEnergy);
823 else if(A>1)
824 return Z*(protonMass - protonSeparationEnergy) + (A-Z)*(neutronMass - neutronSeparationEnergy);
825 else if(A==1 && Z==0 && S==0)
826 return getINCLMass(Neutron);
827 else if(A==1 && Z==1 && S==0)
828 return getINCLMass(Proton);
829 else if(A==1 && Z==0 && S==-1)
830 return getINCLMass(Lambda);
831 else
832 return 0.;
833 }
834
835 G4double getTableQValue(const G4int A1, const G4int Z1, const G4int S1, const G4int A2, const G4int Z2, const G4int S2) {
836 return getTableMass(A1,Z1,S1) + getTableMass(A2,Z2,S2) - getTableMass(A1+A2,Z1+Z2,S1+S2);
837 }
838
839 G4double getTableQValue(const G4int A1, const G4int Z1, const G4int S1, const G4int A2, const G4int Z2, const G4int S2, const G4int A3, const G4int Z3, const G4int S3) {
840 return getTableMass(A1,Z1,S1) + getTableMass(A2,Z2,S2) - getTableMass(A3,Z3,S3) - getTableMass(A1+A2-A3,Z1+Z2-Z3,S1+S2-S3);
841 }
842
843 G4double getTableSpeciesMass(const ParticleSpecies &p) {
844 if(p.theType == Composite)
845 return (*getTableMass)(p.theA, p.theZ, p.theS);
846 else
847 return (*getTableParticleMass)(p.theType);
848 }
849
850 G4int getMassNumber(const ParticleType t) {
851 switch(t) {
852 case Proton:
853 case Neutron:
854 case DeltaPlusPlus:
855 case DeltaPlus:
856 case DeltaZero:
857 case DeltaMinus:
858 case SigmaPlus:
859 case SigmaZero:
860 case SigmaMinus:
861 case Lambda:
862 return 1;
863 break;
864 case PiPlus:
865 case PiMinus:
866 case PiZero:
867 case KPlus:
868 case KZero:
869 case KZeroBar:
870 case KShort:
871 case KLong:
872 case KMinus:
873 case Eta:
874 case Omega:
875 case EtaPrime:
876 case Photon:
877 return 0;
878 break;
879 default:
880 return 0;
881 break;
882 }
883 }
884
885 G4int getChargeNumber(const ParticleType t) {
886 switch(t) {
887 case DeltaPlusPlus:
888 return 2;
889 break;
890 case Proton:
891 case DeltaPlus:
892 case PiPlus:
893 case SigmaPlus:
894 case KPlus:
895 return 1;
896 break;
897 case Neutron:
898 case DeltaZero:
899 case PiZero:
900 case SigmaZero:
901 case Lambda:
902 case KZero:
903 case KZeroBar:
904 case KShort:
905 case KLong:
906 case Eta:
907 case Omega:
908 case EtaPrime:
909 case Photon:
910 return 0;
911 break;
912 case DeltaMinus:
913 case PiMinus:
914 case SigmaMinus:
915 case KMinus:
916 return -1;
917 break;
918 default:
919 return 0;
920 break;
921 }
922 }
923
924 G4int getStrangenessNumber(const ParticleType t) {
925 switch(t) {
926 case DeltaPlusPlus:
927 case DeltaPlus:
928 case DeltaZero:
929 case DeltaMinus:
930 case Proton:
931 case Neutron:
932 case PiPlus:
933 case PiZero:
934 case PiMinus:
935 case Eta:
936 case Omega:
937 case EtaPrime:
938 case Photon:
939 return 0;
940 break;
941 case Lambda:
942 case SigmaPlus:
943 case SigmaZero:
944 case SigmaMinus:
945 case KZeroBar:
946 case KMinus:
947 return -1;
948 break;
949 case KPlus:
950 case KZero:
951 return 1;
952 break;
953 case KShort:
954 return 0;
955 break;
956 case KLong:
957 return 0;
958 break;
959 default:
960 return 0;
961 break;
962 }
963 }
964
965 G4double getNuclearRadius(const ParticleType t, const G4int A, const G4int Z) {
966// assert(A>=0);
967 if(A > 19 || (A < 6 && A >= 2)) {
968 // For large (Woods-Saxon or Modified Harmonic Oscillator) or small
969 // (Gaussian) nuclei, the radius parameter is just the nuclear radius
970 return getRadiusParameter(t,A,Z);
971 } else if(A < clusterTableASize && Z>=0 && Z < clusterTableZSize && A >= 6) {
972 const G4double thisRMS = positionRMS[Z][A];
973 if(thisRMS>0.0)
974 return thisRMS;
975 else {
976 INCL_DEBUG("getNuclearRadius: Radius for nucleus A = " << A << " Z = " << Z << " is not available" << '\n'
977 << "returning radius for C12");
978 return positionRMS[6][12];
979 }
980 } else if(A <= 19) {
981 const G4double theRadiusParameter = getRadiusParameter(t, A, Z);
982 const G4double theDiffusenessParameter = getSurfaceDiffuseness(t, A, Z);
983 // The formula yields the nuclear RMS radius based on the parameters of
984 // the nuclear-density function
985 return 1.225*theDiffusenessParameter*
986 std::sqrt((2.+5.*theRadiusParameter)/(2.+3.*theRadiusParameter));
987 } else {
988 INCL_ERROR("getNuclearRadius: No radius for nucleus A = " << A << " Z = " << Z << '\n');
989 return 0.0;
990 }
991 }
992
993 G4double getLargestNuclearRadius(const G4int A, const G4int Z) {
994 return Math::max(getNuclearRadius(Proton, A, Z), getNuclearRadius(Neutron, A, Z));
995 }
996
997 G4double getRadiusParameter(const ParticleType t, const G4int A, const G4int Z) {
998// assert(A>0);
999 if(A > 19) {
1000 // radius fit for lambdas
1001 if(t==Lambda){
1002 G4double r0 = (1.128+0.439*std::pow(A,-2./3.)) * std::pow(A, 1.0/3.0);
1003 return r0;
1004 }
1005 // phenomenological radius fit
1006 G4double r0 = (2.745e-4 * A + 1.063) * std::pow(A, 1.0/3.0);
1007 // HFB calculations
1008 if(getRPCorrelationCoefficient(t)<1.){
1009 G4double r0hfb = HFB::getRadiusParameterHFB(t,A,Z);
1010 if(r0hfb>0.)r0 = r0hfb;
1011 }
1012 //
1013 if(t==Neutron)
1014 r0 += neutronSkin;
1015 return r0;
1016 } else if(A < 6 && A >= 2) {
1017 if(Z<clusterTableZSize && Z>=0) {
1018 const G4double thisRMS = positionRMS[Z][A];
1019 if(thisRMS>0.0)
1020 return thisRMS;
1021 else {
1022 INCL_DEBUG("getRadiusParameter: Radius for nucleus A = " << A << " Z = " << Z << " is not available" << '\n'
1023 << "returning radius for C12");
1024 return positionRMS[6][12];
1025 }
1026 } else {
1027 INCL_DEBUG("getRadiusParameter: Radius for nucleus A = " << A << " Z = " << Z << " is not available" << '\n'
1028 << "returning radius for C12");
1029 return positionRMS[6][12];
1030 }
1031 } else if(A <= 19 && A >= 6) {
1032 if(t==Lambda){
1033 G4double r0 = (1.128+0.439*std::pow(A,-2./3.)) * std::pow(A, 1.0/3.0);
1034 return r0;
1035 }
1036 // HFB calculations
1037 if(getRPCorrelationCoefficient(t)<1.){
1038 G4double r0hfb = HFB::getSurfaceDiffusenessHFB(t,A,Z);
1039 if(r0hfb>0.)return r0hfb;
1040 }
1041 return mediumRadius[A-1];
1042 // return 1.581*mediumDiffuseness[A-1]*(2.+5.*mediumRadius[A-1])/(2.+3.*mediumRadius[A-1]);
1043 } else {
1044 INCL_ERROR("getRadiusParameter: No radius for nucleus A = " << A << " Z = " << Z << '\n');
1045 return 0.0;
1046 }
1047 }
1048
1049 G4double getMaximumNuclearRadius(const ParticleType t, const G4int A, const G4int Z) {
1050 const G4double XFOISA = 8.0;
1051 if(A > 19) {
1052 return getNuclearRadius(t,A,Z) + XFOISA * getSurfaceDiffuseness(t,A,Z);
1053 } else if(A <= 19 && A >= 6) {
1054 return 5.5 + 0.3 * (G4double(A) - 6.0)/12.0;
1055 } else if(A >= 2) {
1056 return getNuclearRadius(t, A, Z) + 4.5;
1057 } else {
1058 INCL_ERROR("getMaximumNuclearRadius : No maximum radius for nucleus A = " << A << " Z = " << Z << '\n');
1059 return 0.0;
1060 }
1061 }
1062
1063 G4double getSurfaceDiffuseness(const ParticleType t, const G4int A, const G4int Z) {
1064 if(A > 19) {
1065 // phenomenological fit
1066 G4double a = 1.63e-4 * A + 0.510;
1067 // HFB calculations
1068 if(getRPCorrelationCoefficient(t)<1.){
1069 G4double ahfb = HFB::getSurfaceDiffusenessHFB(t,A,Z);
1070 if(ahfb>0.)a=ahfb;
1071 }
1072 //
1073 if(t==Lambda){
1074 // Like for neutrons
1075 G4double ahfb = HFB::getSurfaceDiffusenessHFB(Neutron,A,Z);
1076 if(ahfb>0.)a=ahfb;
1077 }
1078 if(t==Neutron)
1079 a += neutronHalo;
1080 return a;
1081 } else if(A <= 19 && A >= 6) {
1082 // HFB calculations
1083 if(getRPCorrelationCoefficient(t)<1.){
1084 G4double ahfb = HFB::getRadiusParameterHFB(t,A,Z);
1085 if(ahfb>0.)return ahfb;
1086 }
1087 return mediumDiffuseness[A-1];
1088 } else if(A < 6 && A >= 2) {
1089 INCL_ERROR("getSurfaceDiffuseness: was called for A = " << A << " Z = " << Z << '\n');
1090 return 0.0;
1091 } else {
1092 INCL_ERROR("getSurfaceDiffuseness: No diffuseness for nucleus A = " << A << " Z = " << Z << '\n');
1093 return 0.0;
1094 }
1095 }
1096
1097 G4double getMomentumRMS(const G4int A, const G4int Z) {
1098// assert(Z>=0 && A>=0 && Z<=A);
1099 return getFermiMomentum(A,Z) * Math::sqrtThreeFifths;
1100 }
1101
1102 G4double getSeparationEnergyINCL(const ParticleType t, const G4int /*A*/, const G4int /*Z*/) {
1103 if(t==Proton)
1104 return theINCLProtonSeparationEnergy;
1105 else if(t==Neutron)
1106 return theINCLNeutronSeparationEnergy;
1107 else if(t==Lambda)
1108 return theINCLLambdaSeparationEnergy;
1109 else {
1110 INCL_ERROR("ParticleTable::getSeparationEnergyINCL : Unknown particle type." << '\n');
1111 return 0.0;
1112 }
1113 }
1114
1115 G4double getSeparationEnergyReal(const ParticleType t, const G4int A, const G4int Z) {
1116 // Real separation energies for all nuclei
1117 if(t==Proton)
1118 return (*getTableParticleMass)(Proton) + (*getTableMass)(A-1,Z-1,0) - (*getTableMass)(A,Z,0);
1119 else if(t==Neutron)
1120 return (*getTableParticleMass)(Neutron) + (*getTableMass)(A-1,Z,0) - (*getTableMass)(A,Z,0);
1121 else if(t==Lambda)
1122 return (*getTableParticleMass)(Lambda) + (*getTableMass)(A-1,Z,0) - (*getTableMass)(A,Z,-1);
1123 else {
1124 INCL_ERROR("ParticleTable::getSeparationEnergyReal : Unknown particle type." << '\n');
1125 return 0.0;
1126 }
1127 }
1128
1129 G4double getSeparationEnergyRealForLight(const ParticleType t, const G4int A, const G4int Z) {
1130 // Real separation energies for light nuclei, fixed values for heavy nuclei
1131 if(Z<clusterTableZSize && A<clusterTableASize)
1132 return getSeparationEnergyReal(t, A, Z);
1133 else
1134 return getSeparationEnergyINCL(t, A, Z);
1135 }
1136
1137 G4double getProtonSeparationEnergy() { return protonSeparationEnergy; }
1138
1139 G4double getNeutronSeparationEnergy() { return neutronSeparationEnergy; }
1140
1141 G4double getLambdaSeparationEnergy() { return lambdaSeparationEnergy; }
1142
1143 void setProtonSeparationEnergy(const G4double sen) { protonSeparationEnergy = sen; }
1144
1145 void setNeutronSeparationEnergy(const G4double sen) { neutronSeparationEnergy = sen; }
1146
1147 void setLambdaSeparationEnergy(const G4double sen) { lambdaSeparationEnergy = sen; }
1148
1149 std::string getElementName(const G4int Z) {
1150 if(Z<1) {
1151 INCL_WARN("getElementName called with Z<1" << '\n');
1152 return elementTable[0];
1153 } else if(Z<elementTableSize)
1154 return elementTable[Z];
1155 else
1156 return getIUPACElementName(Z);
1157 }
1158
1159 std::string getIUPACElementName(const G4int Z) {
1160 std::stringstream elementStream;
1161 elementStream << Z;
1162 std::string elementName = elementStream.str();
1163 std::transform(elementName.begin(), elementName.end(), elementName.begin(), intToIUPAC);
1164 elementName[0] = (char)std::toupper(elementName.at(0));
1165 return elementName;
1166 }
1167
1168 G4int parseElement(std::string pS) {
1169 // Normalize the element name
1170 std::transform(pS.begin(), pS.end(), pS.begin(), ::tolower);
1171 pS[0] = (char)std::toupper(pS[0]);
1172
1173 const std::string *iter = std::find(elementTable, elementTable+elementTableSize, pS);
1174 if(iter != elementTable+elementTableSize)
1175 return G4int(iter - elementTable);
1176 else
1177 return ParticleTable::parseIUPACElement(pS);
1178 }
1179
1180 G4int parseIUPACElement(std::string const &sel) {
1181 // Normalise to lower case
1182 std::string elementName(sel);
1183 std::transform(elementName.begin(), elementName.end(), elementName.begin(), ::tolower);
1184 // Return 0 if the element name contains anything but IUPAC digits
1185 if(elementName.find_first_not_of(elementIUPACDigits)!=std::string::npos)
1186 return 0;
1187 std::transform(elementName.begin(), elementName.end(), elementName.begin(), iupacToInt);
1188 std::stringstream elementStream(elementName);
1189 G4int Z;
1190 elementStream >> Z;
1191 return Z;
1192 }
1193
1194 IsotopicDistribution const &getNaturalIsotopicDistribution(const G4int Z) {
1195 return getNaturalIsotopicDistributions()->getIsotopicDistribution(Z);
1196 }
1197
1198 G4int drawRandomNaturalIsotope(const G4int Z) {
1199 return getNaturalIsotopicDistributions()->drawRandomIsotope(Z);
1200 }
1201
1202 G4double getFermiMomentumConstant(const G4int /*A*/, const G4int /*Z*/) {
1203 return constantFermiMomentum;
1204 }
1205
1206 G4double getFermiMomentumConstantLight(const G4int A, const G4int Z) {
1207// assert(Z>0 && A>0 && Z<=A);
1208 if(Z<clusterTableZSize && A<clusterTableASize) {
1209 const G4double rms = momentumRMS[Z][A];
1210 return ((rms>0.) ? rms : momentumRMS[6][12]) * Math::sqrtFiveThirds;
1211 } else
1212 return getFermiMomentumConstant(A,Z);
1213 }
1214
1215 G4double getFermiMomentumMassDependent(const G4int A, const G4int /*Z*/) {
1216// assert(A>0);
1217 static const G4double alphaParam = 259.416; // MeV/c
1218 static const G4double betaParam = 152.824; // MeV/c
1219 static const G4double gammaParam = 9.5157E-2;
1220 return alphaParam - betaParam*std::exp(-gammaParam*((G4double)A));
1221 }
1222
1223 G4double getRPCorrelationCoefficient(const ParticleType t) {
1224// assert(t==Proton || t==Neutron || t==Lambda);
1225 return rpCorrelationCoefficient[t];
1226 }
1227
1228 G4double getNeutronSkin() { return neutronSkin; }
1229
1230 G4double getNeutronHalo() { return neutronHalo; }
1231
1232 G4ThreadLocal G4double minDeltaMass = 0.;
1233 G4ThreadLocal G4double minDeltaMass2 = 0.;
1234 G4ThreadLocal G4double minDeltaMassRndm = 0.;
1235 G4ThreadLocal NuclearMassFn getTableMass = NULL;
1236 G4ThreadLocal ParticleMassFn getTableParticleMass = NULL;
1237 G4ThreadLocal SeparationEnergyFn getSeparationEnergy = NULL;
1238 G4ThreadLocal FermiMomentumFn getFermiMomentum = NULL;
1239
1240 ParticleType getPionType(const G4int isosp) {
1241// assert(isosp == -2 || isosp == 0 || isosp == 2);
1242 if (isosp == -2) {
1243 return PiMinus;
1244 }
1245 else if (isosp == 0) {
1246 return PiZero;
1247 }
1248 else {
1249 return PiPlus;
1250 }
1251 }
1252
1253 ParticleType getNucleonType(const G4int isosp) {
1254// assert(isosp == -1 || isosp == 1);
1255 if (isosp == -1) {
1256 return Neutron;
1257 }
1258 else {
1259 return Proton;
1260 }
1261 }
1262
1263 ParticleType getDeltaType(const G4int isosp) {
1264// assert(isosp == -3 || isosp == -1 || isosp == 1 || isosp == 3);
1265 if (isosp == -3) {
1266 return DeltaMinus;
1267 }
1268 else if (isosp == -1) {
1269 return DeltaZero;
1270 }
1271 else if (isosp == 1) {
1272 return DeltaPlus;
1273 }
1274 else {
1275 return DeltaPlusPlus;
1276 }
1277 }
1278
1279
1280 ParticleType getSigmaType(const G4int isosp) {
1281// assert(isosp == -2 || isosp == 0 || isosp == 2);
1282 if (isosp == -2) {
1283 return SigmaMinus;
1284 }
1285 else if (isosp == 0) {
1286 return SigmaZero;
1287 }
1288 else {
1289 return SigmaPlus;
1290 }
1291 }
1292
1293 ParticleType getKaonType(const G4int isosp) {
1294// assert(isosp == -1 || isosp == 1);
1295 if (isosp == -1) {
1296 return KZero;
1297 }
1298 else {
1299 return KPlus;
1300 }
1301 }
1302
1303 ParticleType getAntiKaonType(const G4int isosp) {
1304// assert(isosp == -1 || isosp == 1);
1305 if (isosp == -1) {
1306 return KMinus;
1307 }
1308 else {
1309 return KZeroBar;
1310 }
1311 }
1312
1313 G4double getWidth(const ParticleType pt) {
1314// assert(pt == PiPlus || pt == PiMinus || pt == PiZero || pt == Eta || pt == Omega || pt == EtaPrime || pt == KShort || pt == KLong || pt== KPlus || pt == KMinus || pt == Lambda || pt == SigmaPlus || pt == SigmaZero || pt == SigmaMinus);
1315 if(pt == PiPlus) {
1316 return piPlusWidth;
1317 } else if(pt == PiMinus) {
1318 return piMinusWidth;
1319 } else if(pt == PiZero) {
1320 return piZeroWidth;
1321 } else if(pt == Eta) {
1322 return etaWidth;
1323 } else if(pt == Omega) {
1324 return omegaWidth;
1325 } else if(pt == EtaPrime) {
1326 return etaPrimeWidth;
1327 } else if(pt == SigmaPlus) {
1328 return SigmaPlusWidth;
1329 } else if(pt == SigmaZero) {
1330 return SigmaZeroWidth;
1331 } else if(pt == SigmaMinus) {
1332 return SigmaMinusWidth;
1333 } else if(pt == KPlus) {
1334 return KPlusWidth;
1335 } else if(pt == KMinus) {
1336 return KMinusWidth;
1337 } else if(pt == KShort) {
1338 return KShortWidth;
1339 } else if(pt == KLong) {
1340 return KLongWidth;
1341 } else {
1342 INCL_ERROR("getWidth : Unknown particle type." << '\n');
1343 return 0.0;
1344 }
1345 }
1346
1347 } // namespace ParticleTable
1348} // namespace G4INCL
1349
S
G4double S(G4double temp)
Definition: G4DNAElectronHoleRecombination.cc:76
INCL_ERROR
#define INCL_ERROR(x)
Definition: G4INCLLogger.hh:225
INCL_WARN
#define INCL_WARN(x)
Definition: G4INCLLogger.hh:232
INCL_FATAL
#define INCL_FATAL(x)
Definition: G4INCLLogger.hh:216
INCL_DEBUG
#define INCL_DEBUG(x)
Definition: G4INCLLogger.hh:240
G4INCLNuclearMassTable.hh
Functions that encapsulate a mass table.
INCL_DEFAULT_SEPARATION_ENERGY
#define INCL_DEFAULT_SEPARATION_ENERGY
Definition: G4INCLParticleTable.cc:298
G4double
double G4double
Definition: G4Types.hh:83
G4int
int G4int
Definition: G4Types.hh:85
Z
const G4int Z[17]
Definition: G4WaterStopping.cc:51
A
const G4double A[17]
Definition: G4WaterStopping.cc:53
G4INCL::Config::getNeutronHalo
G4double getNeutronHalo() const
Get the neutron-halo size.
Definition: G4INCLConfig.hh:309
G4INCL::Config::getFermiMomentumType
FermiMomentumType getFermiMomentumType() const
Get the Fermi-momentum type.
Definition: G4INCLConfig.hh:252
G4INCL::Config::getSeparationEnergyType
SeparationEnergyType getSeparationEnergyType() const
Get the separation-energy type.
Definition: G4INCLConfig.hh:249
G4INCL::Config::getRPCorrelationCoefficient
G4double getRPCorrelationCoefficient(const ParticleType t) const
Get the r-p correlation coefficient.
Definition: G4INCLConfig.hh:282
G4INCL::Config::getINCLXXDataFilePath
std::string const & getINCLXXDataFilePath() const
Definition: G4INCLConfig.hh:219
G4INCL::Config::getNeutronSkin
G4double getNeutronSkin() const
Get the neutron-skin thickness.
Definition: G4INCLConfig.hh:303
G4INCL::Config::getFermiMomentum
G4double getFermiMomentum() const
Get the Fermi momentum.
Definition: G4INCLConfig.hh:258
G4INCL::Config::getUseRealMasses
G4bool getUseRealMasses() const
Whether to use real masses.
Definition: G4INCLConfig.hh:211
G4INCL::IsotopicDistribution
Class that stores isotopic abundances for a given element.
Definition: G4INCLNaturalIsotopicDistributions.hh:64
G4ParticleTable::GetIonTable
G4IonTable * GetIonTable() const
Definition: G4ParticleTable.cc:698
G4ParticleTable::FindParticle
G4ParticleDefinition * FindParticle(G4int PDGEncoding)
Definition: G4ParticleTable.cc:586
G4ParticleTable::GetParticleTable
static G4ParticleTable * GetParticleTable()
Definition: G4ParticleTable.cc:87
G4INCL::HFB::initialize
void initialize()
Definition: G4INCLHFB.cc:81
G4INCL::HFB::getSurfaceDiffusenessHFB
G4double getSurfaceDiffusenessHFB(const ParticleType t, const G4int A, const G4int Z)
Definition: G4INCLHFB.cc:140
G4INCL::HFB::getRadiusParameterHFB
G4double getRadiusParameterHFB(const ParticleType t, const G4int A, const G4int Z)
Get the radius and diffuseness parameters from HFB calculations.
Definition: G4INCLHFB.cc:130
G4INCL::Math::max
T max(const T t1, const T t2)
brief Return the largest of the two arguments
Definition: G4INCLGlobals.hh:112
G4INCL::Math::sqrtFiveThirds
const G4double sqrtFiveThirds
Definition: G4INCLGlobals.hh:76
G4INCL::Math::sqrtThreeFifths
const G4double sqrtThreeFifths
Definition: G4INCLGlobals.hh:77
G4INCL::ParticleTable::getMassNumber
G4int getMassNumber(const ParticleType t)
Get mass number from particle type.
Definition: G4INCLParticleTable.cc:850
G4INCL::ParticleTable::getFermiMomentum
G4ThreadLocal FermiMomentumFn getFermiMomentum
Definition: G4INCLParticleTable.cc:1238
G4INCL::ParticleTable::effectiveDeltaWidth
const G4double effectiveDeltaWidth
Definition: G4INCLParticleTable.hh:72
G4INCL::ParticleTable::parseElement
G4int parseElement(std::string pS)
Get the name of the element from the atomic number.
Definition: G4INCLParticleTable.cc:1168
G4INCL::ParticleTable::minDeltaMass2
G4ThreadLocal G4double minDeltaMass2
Definition: G4INCLParticleTable.cc:1233
G4INCL::ParticleTable::FermiMomentumFn
G4double(* FermiMomentumFn)(const G4int, const G4int)
Definition: G4INCLParticleTable.hh:232
G4INCL::ParticleTable::getTableMass
G4ThreadLocal NuclearMassFn getTableMass
Static pointer to the mass function for nuclei.
Definition: G4INCLParticleTable.cc:1235
G4INCL::ParticleTable::getSeparationEnergy
G4ThreadLocal SeparationEnergyFn getSeparationEnergy
Static pointer to the separation-energy function.
Definition: G4INCLParticleTable.cc:1237
G4INCL::ParticleTable::getTableQValue
G4double getTableQValue(const G4int A1, const G4int Z1, const G4int S1, const G4int A2, const G4int Z2, const G4int S2)
Get Q-value (in MeV/c^2)
Definition: G4INCLParticleTable.cc:835
G4INCL::ParticleTable::getTableParticleMass
G4ThreadLocal ParticleMassFn getTableParticleMass
Static pointer to the mass function for particles.
Definition: G4INCLParticleTable.cc:1236
G4INCL::ParticleTable::initialize
void initialize(Config const *const theConfig=0)
Initialize the particle table.
Definition: G4INCLParticleTable.cc:336
G4INCL::ParticleTable::effectiveDeltaMass
const G4double effectiveDeltaMass
Definition: G4INCLParticleTable.hh:71
G4INCL::ParticleTable::getFermiMomentumMassDependent
G4double getFermiMomentumMassDependent(const G4int A, const G4int)
Return the value Fermi momentum from a fit.
Definition: G4INCLParticleTable.cc:1215
G4INCL::ParticleTable::getTableSpeciesMass
G4double getTableSpeciesMass(const ParticleSpecies &p)
Definition: G4INCLParticleTable.cc:843
G4INCL::ParticleTable::drawRandomNaturalIsotope
G4int drawRandomNaturalIsotope(const G4int Z)
Definition: G4INCLParticleTable.cc:1198
G4INCL::ParticleTable::getSeparationEnergyReal
G4double getSeparationEnergyReal(const ParticleType t, const G4int A, const G4int Z)
Return the real separation energy.
Definition: G4INCLParticleTable.cc:1115
G4INCL::ParticleTable::getNeutronSeparationEnergy
G4double getNeutronSeparationEnergy()
Getter for neutronSeparationEnergy.
Definition: G4INCLParticleTable.cc:1139
G4INCL::ParticleTable::minDeltaMass
G4ThreadLocal G4double minDeltaMass
Definition: G4INCLParticleTable.cc:1232
G4INCL::ParticleTable::getRadiusParameter
G4double getRadiusParameter(const ParticleType t, const G4int A, const G4int Z)
Definition: G4INCLParticleTable.cc:997
G4INCL::ParticleTable::getLargestNuclearRadius
G4double getLargestNuclearRadius(const G4int A, const G4int Z)
Definition: G4INCLParticleTable.cc:993
G4INCL::ParticleTable::getKaonType
ParticleType getKaonType(const G4int isosp)
Get the type of kaon.
Definition: G4INCLParticleTable.cc:1293
G4INCL::ParticleTable::getNeutronHalo
G4double getNeutronHalo()
Get the size of the neutron halo.
Definition: G4INCLParticleTable.cc:1230
G4INCL::ParticleTable::getRealMass
G4double getRealMass(const G4INCL::ParticleType t)
Get particle mass (in MeV/c^2)
Definition: G4INCLParticleTable.cc:725
G4INCL::ParticleTable::getSigmaType
ParticleType getSigmaType(const G4int isosp)
Get the type of sigma.
Definition: G4INCLParticleTable.cc:1280
G4INCL::ParticleTable::getINCLMass
G4double getINCLMass(const G4int A, const G4int Z, const G4int S)
Get INCL nuclear mass (in MeV/c^2)
Definition: G4INCLParticleTable.cc:809
G4INCL::ParticleTable::ParticleMassFn
G4double(* ParticleMassFn)(const ParticleType)
Definition: G4INCLParticleTable.hh:220
G4INCL::ParticleTable::getStrangenessNumber
G4int getStrangenessNumber(const ParticleType t)
Get strangeness number from particle type.
Definition: G4INCLParticleTable.cc:924
G4INCL::ParticleTable::getMaximumNuclearRadius
G4double getMaximumNuclearRadius(const ParticleType t, const G4int A, const G4int Z)
Definition: G4INCLParticleTable.cc:1049
G4INCL::ParticleTable::getRPCorrelationCoefficient
G4double getRPCorrelationCoefficient(const ParticleType t)
Get the value of the r-p correlation coefficient.
Definition: G4INCLParticleTable.cc:1223
G4INCL::ParticleTable::parseIUPACElement
G4int parseIUPACElement(std::string const &pS)
Parse a IUPAC element name.
Definition: G4INCLParticleTable.cc:1180
G4INCL::ParticleTable::getSeparationEnergyINCL
G4double getSeparationEnergyINCL(const ParticleType t, const G4int, const G4int)
Return INCL's default separation energy.
Definition: G4INCLParticleTable.cc:1102
G4INCL::ParticleTable::setNeutronSeparationEnergy
void setNeutronSeparationEnergy(const G4double s)
Setter for protonSeparationEnergy.
Definition: G4INCLParticleTable.cc:1145
G4INCL::ParticleTable::getFermiMomentumConstant
G4double getFermiMomentumConstant(const G4int, const G4int)
Return the constant value of the Fermi momentum.
Definition: G4INCLParticleTable.cc:1202
G4INCL::ParticleTable::getName
std::string getName(const ParticleType t)
Get the native INCL name of the particle.
Definition: G4INCLParticleTable.cc:574
G4INCL::ParticleTable::minDeltaMassRndm
G4ThreadLocal G4double minDeltaMassRndm
Definition: G4INCLParticleTable.cc:1234
G4INCL::ParticleTable::SeparationEnergyFn
G4double(* SeparationEnergyFn)(const ParticleType, const G4int, const G4int)
Definition: G4INCLParticleTable.hh:227
G4INCL::ParticleTable::getNeutronSkin
G4double getNeutronSkin()
Get the thickness of the neutron skin.
Definition: G4INCLParticleTable.cc:1228
G4INCL::ParticleTable::getIUPACElementName
std::string getIUPACElementName(const G4int Z)
Get the name of an unnamed element from the IUPAC convention.
Definition: G4INCLParticleTable.cc:1159
G4INCL::ParticleTable::getIsospin
G4int getIsospin(const ParticleType t)
Get the isospin of a particle.
Definition: G4INCLParticleTable.cc:478
G4INCL::ParticleTable::getNucleonType
ParticleType getNucleonType(const G4int isosp)
Get the type of nucleon.
Definition: G4INCLParticleTable.cc:1253
G4INCL::ParticleTable::getSurfaceDiffuseness
G4double getSurfaceDiffuseness(const ParticleType t, const G4int A, const G4int Z)
Definition: G4INCLParticleTable.cc:1063
G4INCL::ParticleTable::getFermiMomentumConstantLight
G4double getFermiMomentumConstantLight(const G4int A, const G4int Z)
Return the constant value of the Fermi momentum - special for light.
Definition: G4INCLParticleTable.cc:1206
G4INCL::ParticleTable::setProtonSeparationEnergy
void setProtonSeparationEnergy(const G4double s)
Setter for protonSeparationEnergy.
Definition: G4INCLParticleTable.cc:1143
G4INCL::ParticleTable::getPionType
ParticleType getPionType(const G4int isosp)
Get the type of pion.
Definition: G4INCLParticleTable.cc:1240
G4INCL::ParticleTable::getDeltaType
ParticleType getDeltaType(const G4int isosp)
Get the type of delta.
Definition: G4INCLParticleTable.cc:1263
G4INCL::ParticleTable::setLambdaSeparationEnergy
void setLambdaSeparationEnergy(const G4double sen)
Definition: G4INCLParticleTable.cc:1147
G4INCL::ParticleTable::NuclearMassFn
G4double(* NuclearMassFn)(const G4int, const G4int, const G4int)
Definition: G4INCLParticleTable.hh:219
G4INCL::ParticleTable::getChargeNumber
G4int getChargeNumber(const ParticleType t)
Get charge number from particle type.
Definition: G4INCLParticleTable.cc:885
G4INCL::ParticleTable::getProtonSeparationEnergy
G4double getProtonSeparationEnergy()
Getter for protonSeparationEnergy.
Definition: G4INCLParticleTable.cc:1137
G4INCL::ParticleTable::getNaturalIsotopicDistribution
IsotopicDistribution const & getNaturalIsotopicDistribution(const G4int Z)
Definition: G4INCLParticleTable.cc:1194
G4INCL::ParticleTable::getMomentumRMS
G4double getMomentumRMS(const G4int A, const G4int Z)
Return the RMS of the momentum distribution (light clusters)
Definition: G4INCLParticleTable.cc:1097
G4INCL::ParticleTable::getAntiKaonType
ParticleType getAntiKaonType(const G4int isosp)
Get the type of antikaon.
Definition: G4INCLParticleTable.cc:1303
G4INCL::ParticleTable::getSeparationEnergyRealForLight
G4double getSeparationEnergyRealForLight(const ParticleType t, const G4int A, const G4int Z)
Return the real separation energy only for light nuclei.
Definition: G4INCLParticleTable.cc:1129
G4INCL::ParticleTable::getNuclearRadius
G4double getNuclearRadius(const ParticleType t, const G4int A, const G4int Z)
Definition: G4INCLParticleTable.cc:965
G4INCL::ParticleTable::getWidth
G4double getWidth(const ParticleType t)
Get particle width (in s)
Definition: G4INCLParticleTable.cc:1313
G4INCL::ParticleTable::getShortName
std::string getShortName(const ParticleType t)
Get the short INCL name of the particle.
Definition: G4INCLParticleTable.cc:627
G4INCL::ParticleTable::getElementName
std::string getElementName(const G4int Z)
Get the name of the element from the atomic number.
Definition: G4INCLParticleTable.cc:1149
G4INCL::PhysicalConstants::Pf
const G4double Pf
Fermi momentum [MeV/c].
Definition: G4INCLGlobals.hh:57
G4INCL::MassDependentFermiMomentum
@ MassDependentFermiMomentum
Definition: G4INCLConfigEnums.hh:110
G4INCL::ConstantLightFermiMomentum
@ ConstantLightFermiMomentum
Definition: G4INCLConfigEnums.hh:109
G4INCL::RealForLightSeparationEnergy
@ RealForLightSeparationEnergy
Definition: G4INCLConfigEnums.hh:103
G4ThreadLocal
#define G4ThreadLocal
Definition: tls.hh:77