ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4INCLParticleTable.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4INCLParticleTable.cc
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 
54 namespace G4INCL {
55 
56  namespace ParticleTable {
57 
58  namespace {
59 
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 
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 
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 
317  G4ThreadLocal G4double constantFermiMomentum = 0.0;
318 
320  char iupacToInt(char c) {
321  return (char)(((G4int)'0')+elementIUPACDigits.find(c));
322  }
323 
325  char intToIUPAC(char n) { return elementIUPACDigits.at(n); }
326 
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 &s) {
532  if(s.theType==Composite)
533  return getShortName(s.theA,s.theZ);
534  else
535  return getShortName(s.theType);
536  }
537 
538  std::string getName(const ParticleSpecies &s) {
539  if(s.theType==Composite)
540  return getName(s.theA,s.theZ);
541  else
542  return getName(s.theType);
543  }
544 
545  std::string getName(const G4int A, const G4int Z) {
546  std::stringstream stream;
547  stream << getElementName(Z) << "-" << A;
548  return stream.str();
549  }
550 
551  std::string getShortName(const G4int A, const G4int Z) {
552  std::stringstream stream;
553  stream << getElementName(Z);
554  if(A>0)
555  stream << A;
556  return stream.str();
557  }
558 
559  std::string getName(const ParticleType p) {
560  if(p == G4INCL::Proton) {
561  return std::string("proton");
562  } else if(p == G4INCL::Neutron) {
563  return std::string("neutron");
564  } else if(p == G4INCL::DeltaPlusPlus) {
565  return std::string("delta++");
566  } else if(p == G4INCL::DeltaPlus) {
567  return std::string("delta+");
568  } else if(p == G4INCL::DeltaZero) {
569  return std::string("delta0");
570  } else if(p == G4INCL::DeltaMinus) {
571  return std::string("delta-");
572  } else if(p == G4INCL::PiPlus) {
573  return std::string("pi+");
574  } else if(p == G4INCL::PiZero) {
575  return std::string("pi0");
576  } else if(p == G4INCL::PiMinus) {
577  return std::string("pi-");
578  } else if(p == G4INCL::Lambda) {
579  return std::string("lambda");
580  } else if(p == G4INCL::SigmaPlus) {
581  return std::string("sigma+");
582  } else if(p == G4INCL::SigmaZero) {
583  return std::string("sigma0");
584  } else if(p == G4INCL::SigmaMinus) {
585  return std::string("sigma-");
586  } else if(p == G4INCL::KPlus) {
587  return std::string("kaon+");
588  } else if(p == G4INCL::KZero) {
589  return std::string("kaon0");
590  } else if(p == G4INCL::KZeroBar) {
591  return std::string("kaon0bar");
592  } else if(p == G4INCL::KMinus) {
593  return std::string("kaon-");
594  } else if(p == G4INCL::KShort) {
595  return std::string("kaonshort");
596  } else if(p == G4INCL::KLong) {
597  return std::string("kaonlong");
598  } else if(p == G4INCL::Composite) {
599  return std::string("composite");
600  } else if(p == G4INCL::Eta) {
601  return std::string("eta");
602  } else if(p == G4INCL::Omega) {
603  return std::string("omega");
604  } else if(p == G4INCL::EtaPrime) {
605  return std::string("etaprime");
606  } else if(p == G4INCL::Photon) {
607  return std::string("photon");
608  }
609  return std::string("unknown");
610  }
611 
612  std::string getShortName(const ParticleType p) {
613  if(p == G4INCL::Proton) {
614  return std::string("p");
615  } else if(p == G4INCL::Neutron) {
616  return std::string("n");
617  } else if(p == G4INCL::DeltaPlusPlus) {
618  return std::string("d++");
619  } else if(p == G4INCL::DeltaPlus) {
620  return std::string("d+");
621  } else if(p == G4INCL::DeltaZero) {
622  return std::string("d0");
623  } else if(p == G4INCL::DeltaMinus) {
624  return std::string("d-");
625  } else if(p == G4INCL::PiPlus) {
626  return std::string("pi+");
627  } else if(p == G4INCL::PiZero) {
628  return std::string("pi0");
629  } else if(p == G4INCL::PiMinus) {
630  return std::string("pi-");
631  } else if(p == G4INCL::Lambda) {
632  return std::string("l");
633  } else if(p == G4INCL::SigmaPlus) {
634  return std::string("s+");
635  } else if(p == G4INCL::SigmaZero) {
636  return std::string("s0");
637  } else if(p == G4INCL::SigmaMinus) {
638  return std::string("s-");
639  } else if(p == G4INCL::KPlus) {
640  return std::string("k+");
641  } else if(p == G4INCL::KZero) {
642  return std::string("k0");
643  } else if(p == G4INCL::KZeroBar) {
644  return std::string("k0b");
645  } else if(p == G4INCL::KMinus) {
646  return std::string("k-");
647  } else if(p == G4INCL::KShort) {
648  return std::string("ks");
649  } else if(p == G4INCL::KLong) {
650  return std::string("kl");
651  } else if(p == G4INCL::Composite) {
652  return std::string("comp");
653  } else if(p == G4INCL::Eta) {
654  return std::string("eta");
655  } else if(p == G4INCL::Omega) {
656  return std::string("omega");
657  } else if(p == G4INCL::EtaPrime) {
658  return std::string("etap");
659  } else if(p == G4INCL::Photon) {
660  return std::string("photon");
661  }
662  return std::string("unknown");
663  }
664 
665  G4double getINCLMass(const ParticleType pt) {
666  if(pt == Proton) {
667  return protonMass;
668  } else if(pt == Neutron) {
669  return neutronMass;
670  } else if(pt == PiPlus) {
671  return piPlusMass;
672  } else if(pt == PiMinus) {
673  return piMinusMass;
674  } else if(pt == PiZero) {
675  return piZeroMass;
676  } else if(pt == SigmaPlus) {
677  return SigmaPlusMass;
678  } else if(pt == SigmaMinus) {
679  return SigmaMinusMass;
680  } else if(pt == SigmaZero) {
681  return SigmaZeroMass;
682  } else if(pt == Lambda) {
683  return LambdaMass;
684  } else if(pt == KPlus) {
685  return KPlusMass;
686  } else if(pt == KZero) {
687  return KZeroMass;
688  } else if(pt == KZeroBar) {
689  return KZeroBarMass;
690  } else if(pt == KMinus) {
691  return KMinusMass;
692  } else if(pt == KShort) {
693  return KShortMass;
694  } else if(pt == KLong) {
695  return KLongMass;
696  } else if(pt == Eta) {
697  return etaMass;
698  } else if(pt == Omega) {
699  return omegaMass;
700  } else if(pt == EtaPrime) {
701  return etaPrimeMass;
702  } else if(pt == Photon) {
703  return photonMass;
704  } else {
705  INCL_ERROR("getMass : Unknown particle type." << '\n');
706  return 0.0;
707  }
708  }
709 
710  G4double getRealMass(const ParticleType t) {
711  switch(t) {
712  case Proton:
713  return theRealProtonMass;
714  break;
715  case Neutron:
716  return theRealNeutronMass;
717  break;
718  case PiPlus:
719  case PiMinus:
720  return theRealChargedPiMass;
721  break;
722  case PiZero:
723  return theRealPiZeroMass;
724  break;
725  case SigmaPlus:
726  return theRealSigmaPlusMass;
727  break;
728  case SigmaZero:
729  return theRealSigmaZeroMass;
730  break;
731  case SigmaMinus:
732  return theRealSigmaMinusMass;
733  break;
734  case Lambda:
735  return theRealLambdaMass;
736  break;
737  case KPlus:
738  case KMinus:
739  return theRealChargedKaonMass;
740  break;
741  case KZero:
742  case KZeroBar:
743  case KShort:
744  case KLong:
745  return theRealNeutralKaonMass;
746  break;
747  case Eta:
748  return theRealEtaMass;
749  break;
750  case Omega:
751  return theRealOmegaMass;
752  break;
753  case EtaPrime:
754  return theRealEtaPrimeMass;
755  break;
756  case Photon:
757  return theRealPhotonMass;
758  break;
759  default:
760  INCL_ERROR("Particle::getRealMass : Unknown particle type." << '\n');
761  return 0.0;
762  break;
763  }
764  }
765 
766  G4double getRealMass(const G4int A, const G4int Z, const G4int S) {
767 // assert(A>=0);
768  // For nuclei with Z<0 or Z>A, assume that the exotic charge state is due to pions
769  if(Z<0 && S<0)
770  return (A+S)*theRealNeutronMass - S*LambdaMass - Z*getRealMass(PiMinus);
771  else if(Z>A && S<0)
772  return (A+S)*theRealProtonMass - S*LambdaMass + (A+S-Z)*getRealMass(PiPlus);
773  if(Z<0)
774  return (A)*theRealNeutronMass - Z*getRealMass(PiMinus);
775  else if(Z>A)
776  return (A)*theRealProtonMass + (A-Z)*getRealMass(PiPlus);
777  else if(Z==0 && S==0)
778  return A*theRealNeutronMass;
779  else if(A==Z)
780  return A*theRealProtonMass;
781  else if(Z==0 && S<0)
782  return (A+S)*theRealNeutronMass-S*LambdaMass;
783  else if(A>1) {
784 #ifndef INCLXX_IN_GEANT4_MODE
785  return ::G4INCL::NuclearMassTable::getMass(A,Z,S);
786 #else
787  if(S<0) return theG4IonTable->GetNucleusMass(Z,A,std::abs(S)) / MeV;
788  else return theG4IonTable->GetNucleusMass(Z,A) / MeV;
789 #endif
790  } else
791  return 0.;
792  }
793 
794  G4double getINCLMass(const G4int A, const G4int Z, const G4int S) {
795 // assert(A>=0);
796  // For nuclei with Z<0 or Z>A, assume that the exotic charge state is due to pions
797  // Note that S<0 for lambda
798  if(Z<0 && S<0)
799  return (A+S)*neutronMass - S*LambdaMass - Z*getINCLMass(PiMinus);
800  else if(Z>A && S<0)
801  return (A+S)*protonMass - S*LambdaMass + (A+S-Z)*getINCLMass(PiPlus);
802  else if(Z<0)
803  return (A)*neutronMass - Z*getINCLMass(PiMinus);
804  else if(Z>A)
805  return (A)*protonMass + (A-Z)*getINCLMass(PiPlus);
806  else if(A>1 && S<0)
807  return Z*(protonMass - protonSeparationEnergy) + (A+S-Z)*(neutronMass - neutronSeparationEnergy) + std::abs(S)*(LambdaMass - lambdaSeparationEnergy);
808  else if(A>1)
809  return Z*(protonMass - protonSeparationEnergy) + (A-Z)*(neutronMass - neutronSeparationEnergy);
810  else if(A==1 && Z==0 && S==0)
811  return getINCLMass(Neutron);
812  else if(A==1 && Z==1 && S==0)
813  return getINCLMass(Proton);
814  else if(A==1 && Z==0 && S==-1)
815  return getINCLMass(Lambda);
816  else
817  return 0.;
818  }
819 
820  G4double getTableQValue(const G4int A1, const G4int Z1, const G4int S1, const G4int A2, const G4int Z2, const G4int S2) {
821  return getTableMass(A1,Z1,S1) + getTableMass(A2,Z2,S2) - getTableMass(A1+A2,Z1+Z2,S1+S2);
822  }
823 
824  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) {
825  return getTableMass(A1,Z1,S1) + getTableMass(A2,Z2,S2) - getTableMass(A3,Z3,S3) - getTableMass(A1+A2-A3,Z1+Z2-Z3,S1+S2-S3);
826  }
827 
828  G4double getTableSpeciesMass(const ParticleSpecies &p) {
829  if(p.theType == Composite)
830  return (*getTableMass)(p.theA, p.theZ, p.theS);
831  else
832  return (*getTableParticleMass)(p.theType);
833  }
834 
835  G4int getMassNumber(const ParticleType t) {
836  switch(t) {
837  case Proton:
838  case Neutron:
839  case DeltaPlusPlus:
840  case DeltaPlus:
841  case DeltaZero:
842  case DeltaMinus:
843  case SigmaPlus:
844  case SigmaZero:
845  case SigmaMinus:
846  case Lambda:
847  return 1;
848  break;
849  case PiPlus:
850  case PiMinus:
851  case PiZero:
852  case KPlus:
853  case KZero:
854  case KZeroBar:
855  case KShort:
856  case KLong:
857  case KMinus:
858  case Eta:
859  case Omega:
860  case EtaPrime:
861  case Photon:
862  return 0;
863  break;
864  default:
865  return 0;
866  break;
867  }
868  }
869 
870  G4int getChargeNumber(const ParticleType t) {
871  switch(t) {
872  case DeltaPlusPlus:
873  return 2;
874  break;
875  case Proton:
876  case DeltaPlus:
877  case PiPlus:
878  case SigmaPlus:
879  case KPlus:
880  return 1;
881  break;
882  case Neutron:
883  case DeltaZero:
884  case PiZero:
885  case SigmaZero:
886  case Lambda:
887  case KZero:
888  case KZeroBar:
889  case KShort:
890  case KLong:
891  case Eta:
892  case Omega:
893  case EtaPrime:
894  case Photon:
895  return 0;
896  break;
897  case DeltaMinus:
898  case PiMinus:
899  case SigmaMinus:
900  case KMinus:
901  return -1;
902  break;
903  default:
904  return 0;
905  break;
906  }
907  }
908 
909  G4int getStrangenessNumber(const ParticleType t) {
910  switch(t) {
911  case DeltaPlusPlus:
912  case DeltaPlus:
913  case DeltaZero:
914  case DeltaMinus:
915  case Proton:
916  case Neutron:
917  case PiPlus:
918  case PiZero:
919  case PiMinus:
920  case Eta:
921  case Omega:
922  case EtaPrime:
923  case Photon:
924  return 0;
925  break;
926  case Lambda:
927  case SigmaPlus:
928  case SigmaZero:
929  case SigmaMinus:
930  case KZeroBar:
931  case KMinus:
932  return -1;
933  break;
934  case KPlus:
935  case KZero:
936  return 1;
937  break;
938  case KShort:
939  return 0;
940  break;
941  case KLong:
942  return 0;
943  break;
944  default:
945  return 0;
946  break;
947  }
948  }
949 
950  G4double getNuclearRadius(const ParticleType t, const G4int A, const G4int Z) {
951 // assert(A>=0);
952  if(A > 19 || (A < 6 && A >= 2)) {
953  // For large (Woods-Saxon or Modified Harmonic Oscillator) or small
954  // (Gaussian) nuclei, the radius parameter is just the nuclear radius
955  return getRadiusParameter(t,A,Z);
956  } else if(A < clusterTableASize && Z>=0 && Z < clusterTableZSize && A >= 6) {
957  const G4double thisRMS = positionRMS[Z][A];
958  if(thisRMS>0.0)
959  return thisRMS;
960  else {
961  INCL_DEBUG("getNuclearRadius: Radius for nucleus A = " << A << " Z = " << Z << " is not available" << '\n'
962  << "returning radius for C12");
963  return positionRMS[6][12];
964  }
965  } else if(A <= 19) {
966  const G4double theRadiusParameter = getRadiusParameter(t, A, Z);
967  const G4double theDiffusenessParameter = getSurfaceDiffuseness(t, A, Z);
968  // The formula yields the nuclear RMS radius based on the parameters of
969  // the nuclear-density function
970  return 1.225*theDiffusenessParameter*
971  std::sqrt((2.+5.*theRadiusParameter)/(2.+3.*theRadiusParameter));
972  } else {
973  INCL_ERROR("getNuclearRadius: No radius for nucleus A = " << A << " Z = " << Z << '\n');
974  return 0.0;
975  }
976  }
977 
978  G4double getLargestNuclearRadius(const G4int A, const G4int Z) {
979  return Math::max(getNuclearRadius(Proton, A, Z), getNuclearRadius(Neutron, A, Z));
980  }
981 
982  G4double getRadiusParameter(const ParticleType t, const G4int A, const G4int Z) {
983 // assert(A>0);
984  if(A > 19) {
985  // radius fit for lambdas
986  if(t==Lambda){
987  G4double r0 = (1.128+0.439*std::pow(A,-2./3.)) * std::pow(A, 1.0/3.0);
988  return r0;
989  }
990  // phenomenological radius fit
991  G4double r0 = (2.745e-4 * A + 1.063) * std::pow(A, 1.0/3.0);
992  // HFB calculations
993  if(getRPCorrelationCoefficient(t)<1.){
994  G4double r0hfb = HFB::getRadiusParameterHFB(t,A,Z);
995  if(r0hfb>0.)r0 = r0hfb;
996  }
997  //
998  if(t==Neutron)
999  r0 += neutronSkin;
1000  return r0;
1001  } else if(A < 6 && A >= 2) {
1002  if(Z<clusterTableZSize && Z>=0) {
1003  const G4double thisRMS = positionRMS[Z][A];
1004  if(thisRMS>0.0)
1005  return thisRMS;
1006  else {
1007  INCL_DEBUG("getRadiusParameter: Radius for nucleus A = " << A << " Z = " << Z << " is not available" << '\n'
1008  << "returning radius for C12");
1009  return positionRMS[6][12];
1010  }
1011  } else {
1012  INCL_DEBUG("getRadiusParameter: Radius for nucleus A = " << A << " Z = " << Z << " is not available" << '\n'
1013  << "returning radius for C12");
1014  return positionRMS[6][12];
1015  }
1016  } else if(A <= 19 && A >= 6) {
1017  if(t==Lambda){
1018  G4double r0 = (1.128+0.439*std::pow(A,-2./3.)) * std::pow(A, 1.0/3.0);
1019  return r0;
1020  }
1021  // HFB calculations
1022  if(getRPCorrelationCoefficient(t)<1.){
1023  G4double r0hfb = HFB::getSurfaceDiffusenessHFB(t,A,Z);
1024  if(r0hfb>0.)return r0hfb;
1025  }
1026  return mediumRadius[A-1];
1027  // return 1.581*mediumDiffuseness[A-1]*(2.+5.*mediumRadius[A-1])/(2.+3.*mediumRadius[A-1]);
1028  } else {
1029  INCL_ERROR("getRadiusParameter: No radius for nucleus A = " << A << " Z = " << Z << '\n');
1030  return 0.0;
1031  }
1032  }
1033 
1034  G4double getMaximumNuclearRadius(const ParticleType t, const G4int A, const G4int Z) {
1035  const G4double XFOISA = 8.0;
1036  if(A > 19) {
1037  return getNuclearRadius(t,A,Z) + XFOISA * getSurfaceDiffuseness(t,A,Z);
1038  } else if(A <= 19 && A >= 6) {
1039  return 5.5 + 0.3 * (G4double(A) - 6.0)/12.0;
1040  } else if(A >= 2) {
1041  return getNuclearRadius(t, A, Z) + 4.5;
1042  } else {
1043  INCL_ERROR("getMaximumNuclearRadius : No maximum radius for nucleus A = " << A << " Z = " << Z << '\n');
1044  return 0.0;
1045  }
1046  }
1047 
1048  G4double getSurfaceDiffuseness(const ParticleType t, const G4int A, const G4int Z) {
1049  if(A > 19) {
1050  // phenomenological fit
1051  G4double a = 1.63e-4 * A + 0.510;
1052  // HFB calculations
1053  if(getRPCorrelationCoefficient(t)<1.){
1054  G4double ahfb = HFB::getSurfaceDiffusenessHFB(t,A,Z);
1055  if(ahfb>0.)a=ahfb;
1056  }
1057  //
1058  if(t==Lambda){
1059  // Like for neutrons
1060  G4double ahfb = HFB::getSurfaceDiffusenessHFB(Neutron,A,Z);
1061  if(ahfb>0.)a=ahfb;
1062  }
1063  if(t==Neutron)
1064  a += neutronHalo;
1065  return a;
1066  } else if(A <= 19 && A >= 6) {
1067  // HFB calculations
1068  if(getRPCorrelationCoefficient(t)<1.){
1069  G4double ahfb = HFB::getRadiusParameterHFB(t,A,Z);
1070  if(ahfb>0.)return ahfb;
1071  }
1072  return mediumDiffuseness[A-1];
1073  } else if(A < 6 && A >= 2) {
1074  INCL_ERROR("getSurfaceDiffuseness: was called for A = " << A << " Z = " << Z << '\n');
1075  return 0.0;
1076  } else {
1077  INCL_ERROR("getSurfaceDiffuseness: No diffuseness for nucleus A = " << A << " Z = " << Z << '\n');
1078  return 0.0;
1079  }
1080  }
1081 
1082  G4double getMomentumRMS(const G4int A, const G4int Z) {
1083 // assert(Z>=0 && A>=0 && Z<=A);
1084  return getFermiMomentum(A,Z) * Math::sqrtThreeFifths;
1085  }
1086 
1087  G4double getSeparationEnergyINCL(const ParticleType t, const G4int /*A*/, const G4int /*Z*/) {
1088  if(t==Proton)
1089  return theINCLProtonSeparationEnergy;
1090  else if(t==Neutron)
1091  return theINCLNeutronSeparationEnergy;
1092  else if(t==Lambda)
1093  return theINCLLambdaSeparationEnergy;
1094  else {
1095  INCL_ERROR("ParticleTable::getSeparationEnergyINCL : Unknown particle type." << '\n');
1096  return 0.0;
1097  }
1098  }
1099 
1100  G4double getSeparationEnergyReal(const ParticleType t, const G4int A, const G4int Z) {
1101  // Real separation energies for all nuclei
1102  if(t==Proton)
1103  return (*getTableParticleMass)(Proton) + (*getTableMass)(A-1,Z-1,0) - (*getTableMass)(A,Z,0);
1104  else if(t==Neutron)
1105  return (*getTableParticleMass)(Neutron) + (*getTableMass)(A-1,Z,0) - (*getTableMass)(A,Z,0);
1106  else if(t==Lambda)
1107  return (*getTableParticleMass)(Lambda) + (*getTableMass)(A-1,Z,0) - (*getTableMass)(A,Z,-1);
1108  else {
1109  INCL_ERROR("ParticleTable::getSeparationEnergyReal : Unknown particle type." << '\n');
1110  return 0.0;
1111  }
1112  }
1113 
1114  G4double getSeparationEnergyRealForLight(const ParticleType t, const G4int A, const G4int Z) {
1115  // Real separation energies for light nuclei, fixed values for heavy nuclei
1116  if(Z<clusterTableZSize && A<clusterTableASize)
1117  return getSeparationEnergyReal(t, A, Z);
1118  else
1119  return getSeparationEnergyINCL(t, A, Z);
1120  }
1121 
1122  G4double getProtonSeparationEnergy() { return protonSeparationEnergy; }
1123 
1124  G4double getNeutronSeparationEnergy() { return neutronSeparationEnergy; }
1125 
1126  G4double getLambdaSeparationEnergy() { return lambdaSeparationEnergy; }
1127 
1128  void setProtonSeparationEnergy(const G4double s) { protonSeparationEnergy = s; }
1129 
1130  void setNeutronSeparationEnergy(const G4double s) { neutronSeparationEnergy = s; }
1131 
1132  void setLambdaSeparationEnergy(const G4double s) { lambdaSeparationEnergy = s; }
1133 
1134  std::string getElementName(const G4int Z) {
1135  if(Z<1) {
1136  INCL_WARN("getElementName called with Z<1" << '\n');
1137  return elementTable[0];
1138  } else if(Z<elementTableSize)
1139  return elementTable[Z];
1140  else
1141  return getIUPACElementName(Z);
1142  }
1143 
1144  std::string getIUPACElementName(const G4int Z) {
1145  std::stringstream elementStream;
1146  elementStream << Z;
1147  std::string elementName = elementStream.str();
1148  std::transform(elementName.begin(), elementName.end(), elementName.begin(), intToIUPAC);
1149  elementName[0] = std::toupper(elementName.at(0));
1150  return elementName;
1151  }
1152 
1153  G4int parseElement(std::string pS) {
1154  // Normalize the element name
1155  std::transform(pS.begin(), pS.end(), pS.begin(), ::tolower);
1156  pS[0] = ::toupper(pS[0]);
1157 
1158  const std::string *iter = std::find(elementTable, elementTable+elementTableSize, pS);
1159  if(iter != elementTable+elementTableSize)
1160  return iter - elementTable;
1161  else
1162  return ParticleTable::parseIUPACElement(pS);
1163  }
1164 
1165  G4int parseIUPACElement(std::string const &s) {
1166  // Normalise to lower case
1167  std::string elementName(s);
1168  std::transform(elementName.begin(), elementName.end(), elementName.begin(), ::tolower);
1169  // Return 0 if the element name contains anything but IUPAC digits
1170  if(elementName.find_first_not_of(elementIUPACDigits)!=std::string::npos)
1171  return 0;
1172  std::transform(elementName.begin(), elementName.end(), elementName.begin(), iupacToInt);
1173  std::stringstream elementStream(elementName);
1174  G4int Z;
1175  elementStream >> Z;
1176  return Z;
1177  }
1178 
1179  IsotopicDistribution const &getNaturalIsotopicDistribution(const G4int Z) {
1180  return getNaturalIsotopicDistributions()->getIsotopicDistribution(Z);
1181  }
1182 
1183  G4int drawRandomNaturalIsotope(const G4int Z) {
1184  return getNaturalIsotopicDistributions()->drawRandomIsotope(Z);
1185  }
1186 
1187  G4double getFermiMomentumConstant(const G4int /*A*/, const G4int /*Z*/) {
1188  return constantFermiMomentum;
1189  }
1190 
1191  G4double getFermiMomentumConstantLight(const G4int A, const G4int Z) {
1192 // assert(Z>0 && A>0 && Z<=A);
1193  if(Z<clusterTableZSize && A<clusterTableASize) {
1194  const G4double rms = momentumRMS[Z][A];
1195  return ((rms>0.) ? rms : momentumRMS[6][12]) * Math::sqrtFiveThirds;
1196  } else
1197  return getFermiMomentumConstant(A,Z);
1198  }
1199 
1200  G4double getFermiMomentumMassDependent(const G4int A, const G4int /*Z*/) {
1201 // assert(A>0);
1202  static const G4double alphaParam = 259.416; // MeV/c
1203  static const G4double betaParam = 152.824; // MeV/c
1204  static const G4double gammaParam = 9.5157E-2;
1205  return alphaParam - betaParam*std::exp(-gammaParam*((G4double)A));
1206  }
1207 
1208  G4double getRPCorrelationCoefficient(const ParticleType t) {
1209 // assert(t==Proton || t==Neutron);
1210  return rpCorrelationCoefficient[t];
1211  }
1212 
1213  G4double getNeutronSkin() { return neutronSkin; }
1214 
1215  G4double getNeutronHalo() { return neutronHalo; }
1216 
1217  G4ThreadLocal G4double minDeltaMass = 0.;
1218  G4ThreadLocal G4double minDeltaMass2 = 0.;
1219  G4ThreadLocal G4double minDeltaMassRndm = 0.;
1220  G4ThreadLocal NuclearMassFn getTableMass = NULL;
1221  G4ThreadLocal ParticleMassFn getTableParticleMass = NULL;
1222  G4ThreadLocal SeparationEnergyFn getSeparationEnergy = NULL;
1223  G4ThreadLocal FermiMomentumFn getFermiMomentum = NULL;
1224 
1225  ParticleType getPionType(const G4int isosp) {
1226 // assert(isosp == -2 || isosp == 0 || isosp == 2);
1227  if (isosp == -2) {
1228  return PiMinus;
1229  }
1230  else if (isosp == 0) {
1231  return PiZero;
1232  }
1233  else {
1234  return PiPlus;
1235  }
1236  }
1237 
1238  ParticleType getNucleonType(const G4int isosp) {
1239 // assert(isosp == -1 || isosp == 1);
1240  if (isosp == -1) {
1241  return Neutron;
1242  }
1243  else {
1244  return Proton;
1245  }
1246  }
1247 
1248  ParticleType getDeltaType(const G4int isosp) {
1249 // assert(isosp == -3 || isosp == -1 || isosp == 1 || isosp == 3);
1250  if (isosp == -3) {
1251  return DeltaMinus;
1252  }
1253  else if (isosp == -1) {
1254  return DeltaZero;
1255  }
1256  else if (isosp == 1) {
1257  return DeltaPlus;
1258  }
1259  else {
1260  return DeltaPlusPlus;
1261  }
1262  }
1263 
1264 
1265  ParticleType getSigmaType(const G4int isosp) {
1266 // assert(isosp == -2 || isosp == 0 || isosp == 2);
1267  if (isosp == -2) {
1268  return SigmaMinus;
1269  }
1270  else if (isosp == 0) {
1271  return SigmaZero;
1272  }
1273  else {
1274  return SigmaPlus;
1275  }
1276  }
1277 
1278  ParticleType getKaonType(const G4int isosp) {
1279 // assert(isosp == -1 || isosp == 1);
1280  if (isosp == -1) {
1281  return KZero;
1282  }
1283  else {
1284  return KPlus;
1285  }
1286  }
1287 
1288  ParticleType getAntiKaonType(const G4int isosp) {
1289 // assert(isosp == -1 || isosp == 1);
1290  if (isosp == -1) {
1291  return KMinus;
1292  }
1293  else {
1294  return KZeroBar;
1295  }
1296  }
1297 
1298  G4double getWidth(const ParticleType pt) {
1299 // 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);
1300  if(pt == PiPlus) {
1301  return piPlusWidth;
1302  } else if(pt == PiMinus) {
1303  return piMinusWidth;
1304  } else if(pt == PiZero) {
1305  return piZeroWidth;
1306  } else if(pt == Eta) {
1307  return etaWidth;
1308  } else if(pt == Omega) {
1309  return omegaWidth;
1310  } else if(pt == EtaPrime) {
1311  return etaPrimeWidth;
1312  } else if(pt == SigmaPlus) {
1313  return SigmaPlusWidth;
1314  } else if(pt == SigmaZero) {
1315  return SigmaZeroWidth;
1316  } else if(pt == SigmaMinus) {
1317  return SigmaMinusWidth;
1318  } else if(pt == KPlus) {
1319  return KPlusWidth;
1320  } else if(pt == KMinus) {
1321  return KMinusWidth;
1322  } else if(pt == KShort) {
1323  return KShortWidth;
1324  } else if(pt == KLong) {
1325  return KLongWidth;
1326  } else {
1327  INCL_ERROR("getWidth : Unknown particle type." << '\n');
1328  return 0.0;
1329  }
1330  }
1331 
1332  } // namespace ParticleTable
1333 } // namespace G4INCL
1334